MySQL8.0.24の新しいGIS関数(2)~ST_LineInterpolatePoints()を試す

 ひとつ前のエントリ「MySQL8.0.24の新しいGIS関数(1)~LINESTRINGの経路点を求める~」で、LINESTRINGを 指定した割合ごとに区切ってPOINT集合を返す、ST_LineInterpolatePoints()関数を紹介しました。
 MySQL 8.0 のGIS機能(spatial機能)のウリは「測地系に対応したこと」ですので、この ST_LineInterpolatePoints()関数も測地系に対応しているのだと思います。平面座標と地球上の座標では「直線」の捉え方に差が出てきますので、このエントリではそのあたりを試してみたいと思います。

題材とするLINESTRING

 LINESTRING は2個以上の点(POINT)をつないでできる線です。2個以上なので、3個、4個、、、の点を順次つないだものを扱えるのですが、今回は2点間のシンプルなLINESTRINGで試してみることにします。扱う2点は、成田空港(日本)と シャルル・ド・ゴール空港(フランス)とします。測地系WGS84 (4326)です。

空港名 北緯 東経
成田空港 35.793537 度 140.386548 度
シャルルドゴール空港 48.997026 度 2.281060 度

単純にLINESTRINGを表示してみる

 MySQL Workbench を使って、まずこの2点を結ぶLINESTRINGを表示してみます。

SELECT ST_GeomFromText(
  'LINESTRING(35.793537 140.386548, 48.997026 2.281060)', 4326
)

f:id:sakaik:20210503114225p:plain
 少しわかりにくいのですが、東京(右のほう、北緯36度くらい)から、パリ(真ん中あたり、北緯49度くらい)に向けて直線が表示されています。実際の地球上の最短経路が表示されているわけではないようです。

経路点で区切って表示してみる

 そこで、経路点の集合を返してくれる ST_LineInterpolatePoints()関数を使って、このLINESTRINGを100分割したPOINTを表示してみましょう。 先ほどの ST_GeomFromText() を使って作成した SRID 4326 のLINESTRINGに対して、ST_LineInterpolatePoints()関数をカマしてみます。100分割なので第二引数には 0.01 を与えます。

SELECT ST_LineInterpolatePoints(
   ST_GeomFromText(
     'LINESTRING(35.793537 140.386548, 48.997026 2.281060 )', 4326
  ), 0.01
)

f:id:sakaik:20210503114630p:plain

 美しいですね。ST_LineInterpolatePoints() は(そしておそらくは ST_LineInterpolatePoint()も)、きちんと地球が丸いことを知った上で経路上の点を求めてくれていることがわかりました。


ほんとに測地系を考慮した結果なのか?

 こういった実験は、予測した結果が「得られた」ケースとともに「得られなかった」ケースについても確認することで、より内容に確信が持てるようになるものです。ということで、平面座標である SRID 0 について、同じ2点間の座標の Interpolate を求めてみることにします。
 SQLは以下の通り。測地系の指定がなくなったことと、SRID 0 は Axis が SRID 4326 の場合とは異なるので x, y を入れ替えている点が先ほどとの違いです。
 

SELECT ST_LineInterpolatePoints(
   ST_GeomFromText(
     'LINESTRING(140.386548 35.793537, 2.281060 48.997026)'
  ), 0.01
)

 結果は・・・・
f:id:sakaik:20210503115443p:plain
 直線になりました! SRID 4326 でのあの曲線は、ちゃんと測地系を考慮してくれた結果だということが確認できました。

余談:日付変更線を越える場合

 MySQLの中での(SRIDを指定しての)地球は「丸い」ので、東経 180度と 西経180度はつながっています。それを確認してみましょう。サンフランシスコ国際空港(北緯 37.621438 度、西経 -122.376194 度)と成田空港との経路点を表示させてみます。

SELECT ST_LineInterpolatePoints(
   ST_GeomFromText(
     'LINESTRING(35.793537 140.386548, 37.621438 -122.376194)', 4326
  ), 0.01
)

 MySQL Workbench の Spatial View 機能がちょっとイケてなくて、全体を一画面中に表示することができない(縮尺変更が期待通りに動作しない)のですが、以下のように 東京(本当は千葉)を出発して、(緯線に対してやや北よりの進路を取り)日付変更線を越えて、サンフランシスコに到着していることが確認できます。
f:id:sakaik:20210503120043p:plain
f:id:sakaik:20210503120105p:plain




 得られた MULTIPOINT を CAST()関数を使って LINESTRING に戻してあげれば、元のLINESTRINGを刻んだLINESTRINGを生成することができ、例えば フレシェ値(あるいはハウスドルフ値)のより精度の高い結果を得るのに役に立つのかもしれません。
CAST()の使用例については私もまだ試すことができていないので、おいおい挑戦してみたいと思います。


追記

 折角なので CAST() で MULTIPOINT を LINESTRING に変換するのもやってみました。
分割数が多すぎると「きれいに」見えてしまうので、敢えて成田からシャゴールまでを 0.2 ごとに分割する例で試しました。

SELECT CAST(ST_LineInterpolatePoints(
   ST_GeomFromText(
     'LINESTRING(35.793537 140.386548, 48.997026 2.281060 )', 4326
  ), 0.2)
  AS LINESTRING )

f:id:sakaik:20210503122328p:plain

SQL入力エリアに赤バッテンがついていますが、実行自体は正しくできます(MySQL Workbench が、まだこの構文、または LINESTRING という型を認識していないということでしょうか)。
 結果は、カクッカクとした点を経由したLINESTRINGとなったことが見てわかります。
ただ、ひとつ重大な欠点に気づいてしまいまして、、、、、このLINESTRINGは、というかその元となったMULTIPOINTは、「始点を含まない」ということです。0.2ずつに区切る場合は、始点から 0.2 の位置、0.4の位置、、、、を MULTIPOINT として返しますから、0.0 の位置(始点)がここには含まれないのです。
 この用途で使用するときには一工夫が必要となりそうです。(元のLINESTRINGの始点を、InterpolatePoints で得られた MULTIPOINT の先頭にくっつけた上で CASTする等)