untitled .engineer

技術系のブログ(仮)

MySQL8.0でGIS機能を試す No.4 - 投影座標系のデータを使ってみる


目次


本エントリの概要

  • 今まで地理座標系で登録したデータをいじくっていましたが、投影座標系を使ってみるテストをします。
  • このエントリもわかっていないことが多いままです。どこか勘違いしている自信があります。

前回まで

前回試したのは以下の構文。

INSERT INTO portals
(SELECT
  `jiscode`,
  `building`,
  `address`,
  `lat`,
  `lon`,
  ST_GeomFromText(CONCAT('POINT(',`lat`,' ',`lon`,')'), 4612)
FROM local_governments);

地理座標系なのでどのエリアもSRID=4612で統一で問題ありません。
後で計算をする項目についてSRID=4612でいくつかサンプルを取ります。

北海道庁沖縄県庁の距離をST_Distance()で測ってみます。

SET @geom1 = (SELECT geom FROM portals WHERE `id` = 1000);
SET @geom2 = (SELECT geom FROM portals WHERE `id` = 47000);
SELECT ST_Distance(@geom1, @geom2);
+-----------------------------+
| ST_Distance(@geom1, @geom2) |
+-----------------------------+
|          2243864.5690358933 |
+-----------------------------+
1 row in set (0.00 sec)

2243km、妥当な距離です。

北海道庁東京都庁の距離をST_Distance()で測ってみます。

SET @geom1 = (SELECT geom FROM portals WHERE `id` = 1000);
SET @geom2 = (SELECT geom FROM portals WHERE `id` = 13000);
SELECT ST_Distance(@geom1, @geom2);
+-----------------------------+
| ST_Distance(@geom1, @geom2) |
+-----------------------------+
|           831044.6635025819 |
+-----------------------------+
1 row in set (0.00 sec)

831km、妥当な距離です。

投影座標系

さて、前回のエントリで、ST_Union()が地理座標系のGEOMETRY型を引数に取れなかったので、投影座標系を試してみます。
投影座標系の場合エリアによってSRIDが違います。
これを正しく?投入するために、条件分岐してINSERTする方法を考えます。

JGD2000を利用する場合、UTM座標系と平面直角座標系を選択できます。
平面直角座標系の場合は地方自治境で分割されるのでそれで条件分岐すればいいのですが、 2453~2461が離島部で「北海道の一部」とか「東京都の一部」とか「沖縄県の一部」でここの処理がちょっと面倒なのに対して、UTM座標系は東経だけで条件分岐できるので楽です。

JGD2000のUTM座標系は以下のようになっています。

ゾーン SRID 区域
51 3097 東経120-126
52 3098 東経126-132
53 3099 東経132-138
54 3100 東経138-144
55 3101 東経144-150

SRIDが違う=地図が違うということだと思うので、同一のテーブル内に複数のSRIDのレコードがあるのは意味がないように思いますが、今回は試し。
経度線直上の場合にどちらに属するのかわかりませんが、今回は一旦左に寄せて以下のように雑に取り込みます。

DELETE FROM portals;
INSERT INTO portals
(SELECT
  `jiscode`,
  `building`,
  `address`,
  `lat`,
  `lon`,
  ST_GeomFromText(CONCAT('POINT(',`lat`,' ',`lon`,')'), 
    CASE WHEN (`lon` >= 120 AND `lon` < 126) THEN 3097 ELSE
    CASE WHEN (`lon` >= 126 AND `lon` < 132) THEN 3098 ELSE
    CASE WHEN (`lon` >= 132 AND `lon` < 138) THEN 3099 ELSE
    CASE WHEN (`lon` >= 138 AND `lon` < 144) THEN 3100 ELSE
    CASE WHEN (`lon` >= 144 AND `lon` < 150) THEN 3101 ELSE
    3100 END END END END END
  )
FROM local_governments);

データを見てみます。
都道府県庁のjiscodeがn*1000であることを利用して抜粋します。
東経144-150には都道府県庁がありませんので、念のため網走市を足します。

SELECT id, name, ST_SRID(geom) FROM portals WHERE id LIKE '%000' OR name LIKE '網走市役所';
+-------+-----------------+---------------+
| id    | name            | ST_SRID(geom) |
+-------+-----------------+---------------+
|  1000 | 北海道庁        |          3100 |
|  1211 | 網走市役所      |          3101 |
|  2000 | 青森県庁        |          3100 |
|  3000 | 岩手県庁        |          3100 |
|  4000 | 宮城県庁        |          3100 |
|  5000 | 秋田県庁        |          3100 |
|  6000 | 山形県庁        |          3100 |
|  7000 | 福島県庁        |          3100 |
|  8000 | 茨城県庁        |          3100 |
|  9000 | 栃木県庁        |          3100 |
| 10000 | 群馬県庁        |          3100 |
| 11000 | 埼玉県庁        |          3100 |
| 12000 | 千葉県庁        |          3100 |
| 13000 | 東京都庁        |          3100 |
| 14000 | 神奈川県庁      |          3100 |
| 15000 | 新潟県庁        |          3100 |
| 16000 | 富山県庁        |          3099 |
| 17000 | 石川県庁        |          3099 |
| 18000 | 福井県庁        |          3099 |
| 19000 | 山梨県庁        |          3100 |
| 20000 | 長野県庁        |          3100 |
| 21000 | 岐阜県庁        |          3099 |
| 22000 | 静岡県庁        |          3100 |
| 23000 | 愛知県庁        |          3099 |
| 24000 | 三重県庁        |          3099 |
| 25000 | 滋賀県庁        |          3099 |
| 26000 | 京都府庁        |          3099 |
| 27000 | 大阪府庁        |          3099 |
| 28000 | 兵庫県庁        |          3099 |
| 29000 | 奈良県庁        |          3099 |
| 30000 | 和歌山県庁      |          3099 |
| 31000 | 鳥取県庁        |          3099 |
| 32000 | 島根県庁        |          3099 |
| 33000 | 岡山県庁        |          3099 |
| 34000 | 広島県庁        |          3099 |
| 35000 | 山口県庁        |          3098 |
| 36000 | 徳島県庁        |          3099 |
| 37000 | 香川県庁        |          3099 |
| 38000 | 愛媛県庁        |          3099 |
| 39000 | 高知県庁        |          3099 |
| 40000 | 福岡県庁        |          3098 |
| 41000 | 佐賀県庁        |          3098 |
| 42000 | 長崎県庁        |          3098 |
| 43000 | 熊本県庁        |          3098 |
| 44000 | 大分県庁        |          3098 |
| 45000 | 宮崎県庁        |          3098 |
| 46000 | 鹿児島県庁      |          3098 |
| 47000 | 沖縄県庁        |          3098 |
+-------+-----------------+---------------+
48 rows in set (0.02 sec)

想定した通り取り込めています。

ST_Distance()で距離を測る

ゾーンをまたぐ測定

参考までにゾーンが違う北海道庁沖縄県庁の距離をST_Distance()で測ってみます。

SET @geom1 = (SELECT geom FROM portals WHERE `id` = 1000);
SET @geom2 = (SELECT geom FROM portals WHERE `id` = 47000);
SELECT ST_Distance(@geom1, @geom2);
ERROR 3033 (HY000): Binary geometry function st_distance given two geometries of different srids: 3100 and 3098, which should have been identical.

やはりエラーです。同じグループの投影座標系であってもSRIDが違うとダメらしいです。
ST_Transform()も実装されていないので変換もできず、ゾーンをまたいでの距離計測はできそうにありません。

同一ゾーン内の測定

では、同じゾーンの中で北海道庁東京都庁の距離を測ってみます。

SET @geom1 = (SELECT geom FROM portals WHERE `id` = 1000);
SET @geom2 = (SELECT geom FROM portals WHERE `id` = 13000);
SELECT ST_Distance(@geom1, @geom2);
+-----------------------------+
| ST_Distance(@geom1, @geom2) |
+-----------------------------+
|           7.558147767060658 |
+-----------------------------+
1 row in set (0.02 sec)

7.5m。とっても近いですね。(違

マニュアルを見ると

Returns the distance between g1 and g2, measured in the length unit of the spatial reference system (SRS).

Google翻訳

空間参照系(SRS)の長さ単位で測定されたg1とg2の間の距離を返します。

length unit→長さ単位になっているわけですが、これを探し出してみます。

@sakaik さんに教えていただいた「st_spatial_reference_systems」を見てみます。

SELECT DEFINITION FROM INFORMATION_SCHEMA.st_spatial_reference_systems WHERE srs_id = 3100\G

出力をパースすると以下の通りです。

PROJCS["JGD2000 / UTM zone 54N",
  GEOGCS["JGD2000",
    DATUM["Japanese Geodetic Datum 2000",
      SPHEROID["GRS 1980",6378137,298.257222101,
        AUTHORITY["EPSG","7019"]
      ],
      TOWGS84[0,0,0,0,0,0,0],
      AUTHORITY["EPSG","6612"]
    ],
    PRIMEM["Greenwich",0,
      AUTHORITY["EPSG","8901"]
    ],
    UNIT["degree",0.017453292519943278,
      AUTHORITY["EPSG","9122"]
    ],
    AXIS["Lat",NORTH],
    AXIS["Lon",EAST],
    AUTHORITY["EPSG","4612"]
  ],
  PROJECTION["Transverse Mercator",
    AUTHORITY["EPSG","9807"]
  ],
  PARAMETER["Latitude of natural origin",0,
    AUTHORITY["EPSG","8801"]
  ],
  PARAMETER["Longitude of natural origin",141,
    AUTHORITY["EPSG","8802"]
  ],
  PARAMETER["Scale factor at natural origin",0.9996,
    AUTHORITY["EPSG","8805"]
  ],
  PARAMETER["False easting",500000,
    AUTHORITY["EPSG","8806"]
  ],
  PARAMETER["False northing",0,
    AUTHORITY["EPSG","8807"]
  ],
  UNIT["metre",1,
    AUTHORITY["EPSG","9001"]
  ],
  AXIS["E",EAST],
  AXIS["N",NORTH],
  AUTHORITY["EPSG","3100"]
]

UNIT["degree",0.017453292519943278 っていうのと UNIT["metre",1, っていうのがあります。
どっちなんでしょうね。。。。

仮に前者のdegreeだとすると「度(°)」です。
なのでこの場合 約 7.55° ということになります。
「地球は球体じゃない」と何度も口酸っぱく言われたばかりなのですが、7.55°がどのくらいの距離かを仮にうっかり地球が球体だとした場合の計算で求めると、

  • 地球の外周は40,075km
  • 40075/360≒111.3 で1°≒111.3km
  • 7.55*111.3 ≒ 840.3km

SRID=4612で求めた831kmに少し近い数値となりましたので、今回のST_Distance()の戻り値の単位は"degree"であろうと思われます。

ちなみにSRID=4612の単位はどうなってるのでしょう。

SELECT DEFINITION FROM INFORMATION_SCHEMA.st_spatial_reference_systems WHERE srs_id = 4612\G
GEOGCS["JGD2000",
  DATUM["Japanese Geodetic Datum 2000",
    SPHEROID["GRS 1980",6378137,298.257222101,
      AUTHORITY["EPSG","7019"]
    ],
    TOWGS84[0,0,0,0,0,0,0],
    AUTHORITY["EPSG","6612"]
  ],
  PRIMEM["Greenwich",0,
    AUTHORITY["EPSG","8901"]
  ],
  UNIT["degree",0.017453292519943278,
    AUTHORITY["EPSG","9122"]
  ],
  AXIS["Lat",NORTH],
  AXIS["Lon",EAST],
  AUTHORITY["EPSG","4612"]
]

こっちはUNITが"degree"しかありません。なのにST_Distance()の戻り値はどう見てもメートルです。

これはマニュアルからは読み取れなかったのですが、たまたま見ていたPostGISのST_Distance()のドキュメントで、
PostGIS 2.4.5dev Manual 8.9. Spatial Relationships and Measurements https://postgis.net/docs/ST_Distance.html
原文

ST_Distance - For geometry type Returns the 2D Cartesian distance between two geometries in projected units (based on spatial ref). For geography type defaults to return minimum geodesic distance between two geographies in meters.

Google翻訳

ST_Distance - ジオメトリ・タイプの場合2つのジオメトリ間の2次元デカルト距離を投影単位で返します(空間参照に基づく)。地理タイプのデフォルトでは、2つの地域間の最小測地線距離がメートル単位で返されます。

とありますので、MySQL8.0の場合も地理座標系の場合はメートル単位になるということなのではないかと推測します。OpenGISの仕様なのかどうかまでは調べてません。

ST_Union()で結合を試す

前回地理座標系のSRID=4612のデータではできなかったST_Union()を試してみようと思います。

SET @geom1 = (SELECT geom FROM portals WHERE `id` = 13214);
SET @geom2 = (SELECT geom FROM portals WHERE `id` = 13218);
SELECT ST_AsText(ST_Union(@geom1, @geom2));
+-----------------------------------------------------------+
| ST_AsText(ST_Union(@geom1, @geom2))                       |
+-----------------------------------------------------------+
| MULTIPOINT((35.710948 139.462236),(35.738692 139.326691)) |
+-----------------------------------------------------------+
1 row in set (0.00 sec)

あっさり成功です。ついでにST_ConvexHull()も試します。

SELECT ST_AsText(ST_ConvexHull(ST_Union(@geom1, @geom2)));
+-------------------------------------------------------+
| ST_AsText(ST_ConvexHull(ST_Union(@geom1, @geom2)))    |
+-------------------------------------------------------+
| LINESTRING(35.710948 139.462236,35.738692 139.326691) |
+-------------------------------------------------------+
1 row in set (0.00 sec)

問題なさそうです。
SRIDはどうでしょうか。

SELECT ST_SRID(ST_ConvexHull(ST_Union(@geom1, @geom2)));
+--------------------------------------------------+
| ST_SRID(ST_ConvexHull(ST_Union(@geom1, @geom2))) |
+--------------------------------------------------+
|                                             3100 |
+--------------------------------------------------+

こちらも大丈夫です。

さらに「Ingressリンク可能なリンク先の選別」を試してみます。

MySQL8.0でGIS機能を試す No.2 - Ingressのリンク可能なリンク先の選別 - untitled .engineer

これも同じゾーン内の比較をするため、SIRD=3100以外のデータをportalsから削除しておきます。

DELETE FROM portals WHERE ST_SRID(geom) != 3100;
Query OK, 1050 rows affected (0.32 sec)

linksテーブルにはSRID=3100でデータを入れなおします。

DELETE FROM links;
SET @srid=3100;
INSERT INTO gis.links VALUES (null,13104,13103,ST_GeomFromText(CONCAT((SELECT CONCAT('LINESTRING(',lat,' ',lon,',') FROM portals WHERE `id` = 13104),(SELECT CONCAT(lat,' ',lon,')') FROM portals WHERE `id` = 13103)),3100));
INSERT INTO gis.links VALUES (null,13103,13106,ST_GeomFromText(CONCAT((SELECT CONCAT('LINESTRING(',lat,' ',lon,',') FROM portals WHERE `id` = 13103),(SELECT CONCAT(lat,' ',lon,')') FROM portals WHERE `id` = 13106)),3100));
INSERT INTO gis.links VALUES (null,13106,13120,ST_GeomFromText(CONCAT((SELECT CONCAT('LINESTRING(',lat,' ',lon,',') FROM portals WHERE `id` = 13106),(SELECT CONCAT(lat,' ',lon,')') FROM portals WHERE `id` = 13120)),3100));
INSERT INTO gis.links VALUES (null,13112,13120,ST_GeomFromText(CONCAT((SELECT CONCAT('LINESTRING(',lat,' ',lon,',') FROM portals WHERE `id` = 13112),(SELECT CONCAT(lat,' ',lon,')') FROM portals WHERE `id` = 13120)),3100));
INSERT INTO gis.links VALUES (null,13113,13112,ST_GeomFromText(CONCAT((SELECT CONCAT('LINESTRING(',lat,' ',lon,',') FROM portals WHERE `id` = 13113),(SELECT CONCAT(lat,' ',lon,')') FROM portals WHERE `id` = 13112)),3100));

距離のメートル単位の指定ができないので、10°未満で見てみます。

SET @srid=3100;
SELECT
  PT.name,
  ST_Distance(PF.geom, PT.geom) AS DIST
FROM portals PT, portals PF
WHERE PF.id = 13104
AND PT.id != 13104
AND ST_Distance(PF.geom, PT.geom) < 10
AND NOT EXISTS(
  SELECT `id`
  FROM links L
  WHERE 
    ST_Crosses(
      ST_GeomFromText(
        CONCAT(
          'LINESTRING(',PF.lat,' ',PF.lon,',',PT.lat,' ',PT.lon,')'
        )
      ,@srid),
      L.geom
    ) = 1)
AND NOT EXISTS(
  SELECT `id`
  FROM links L
  WHERE
    ST_Equals(
      ST_GeomFromText(
        CONCAT(
          'LINESTRING(',PF.lat,' ',PF.lon,',',PT.lat,' ',PT.lon,')'
        )
      ,@srid),
      L.geom
    ) = 1)
ORDER BY ST_Distance(PF.geom, PT.geom);
+--------------------+---------------------+
| name               | DIST                |
+--------------------+---------------------+
| 東京都庁           | 0.01260299456476942 |
| 渋谷区役所         | 0.03199157615685517 |
| 中野区役所         | 0.04182679796492804 |
| 千代田区役所       | 0.05011300360186995 |
| 文京区役所         | 0.05070790677004672 |
| 練馬区役所         | 0.06606636471004149 |
| 世田谷区役所       | 0.06941099015285372 |
| 台東区役所         | 0.07873210499535442 |
| 品川区役所         | 0.08910147312474771 |
| 大田区役所         | 0.13318854164304347 |
| 袖ケ浦市役所       | 0.36437785659394795 |
| 木更津市役所       |  0.3831966879188805 |
| 君津市役所         |  0.4142333975731584 |
| 富津市役所         |  0.4189912983022505 |
| 鋸南町役場         |  0.5977174078943974 |
| 南房総市役所       |  0.6649814985651841 |
| 大多喜町役場       |  0.6790321828845963 |
| 鴨川市役所         |  0.7018729573683606 |
| 館山市役所         |  0.7169419424876777 |
| 御宿町役場         |  0.8172220519938307 |
| 勝浦市役所         |  0.8212857028616766 |
| 八丈町役場         |  2.5827249671749817 |
| 青ヶ島村役場       |  3.2276579056202657 |
| 小笠原村役場       |   8.952218327483864 |
+--------------------+---------------------+
24 rows in set (0.85 sec)

小笠原まで届きました。
距離の妥当性はすぐにはわかりませんが、候補とその順序は同等です。
パフォーマンスも少し改善しています。(SRID=4612で距離1000km未満のとき:1.24 secでした)

(投影座標系をいじってみた)今回のまとめ

  • ST_Union()はちゃんと動く
  • ST_Distance()の結果は使いづらい
  • ゾーン(SRID)をまたいでの比較や計測などはできない
  • やっぱりGoogle翻訳機械翻訳のレベルを超えている

次回はもう一度Ingressのルールに基づく検証に戻ります。