Zenn
💭

ポリゴンをいい感じに地理座標系変換したい

2021/01/14に公開
GeoJSON
geo4326
tech

衛星データを扱っていると、EPSG:4326(地理座標系、経度緯度で座標を表す)以外のポリゴンをEPSG:4326に変換してGeoJSONを作りたくなることがままあります。

例えば人工衛星が撮影した画像の撮像範囲をGeoJSONにしたいときとか。

このとき、地理座標系では経度が180degを超えると-180degへワープするので、できればそれも考慮して変換してほしい。

そしてどうせならGeoJSONもRFCを満たす(Polygonが反時計回りだったり180degをまたぐ場合はMultiPolygonで表現したり)ようにしたい。

といったようにポリゴンをいい感じに地理座標系変換してGeoJSONを作ってくれるライブラリを探したのですが、見つけることができませんでした(あったら教えて下さい)[1]。

プロはだいたいArcGISという有償ソフトウェアを使っているらしい[2]。

無いものは仕方ないので自作してみましょう。

条件

  • 出力するGeoJSONはRFC7946を満たす。
  • 投影による辺の湾曲は辺を分割していい感じに近似する。
  • 180degをまたぐケースも正しく計算できることを目指す。
  • JavaScriptで実装できるようライブラリへの依存は最小限にする。
  • 両極点を含むポリゴンの変換は諦める(理由は後述)。

処理フロー

  1. 投影前の座標のLinearRing(先頭と末尾の点が等しい値となる[x, y]の配列、反時計回り)を作る
  2. 両極点を含む場合エラーとする
  3. 自己交差していたり潰れている場合エラーとする
  4. LinearRingの各辺をN分割する
  5. LinearRingをEPSG:4326に投影変換
  6. 投影変換したLinearRingが反時計回りか判定する

投影変換したLinearRingが反時計回りのとき

  1. 素直にbbox=[xの最小値, yの最小値, xの最大値, yの最大値]を計算し、GeoJSONを作成する。

投影変換したLinearRingが時計回りのとき(180degをまたいでいる)

  1. 経度方向の長さが180以上の辺は、180degをまたいでいると判定して180degとの交点を挿入する。
  2. 180degより左側と右側にポリゴンを分類する。
  3. bbox=[左側のxの最小値, yの最小値, 右側のxの最小値, yの最小値]とし、GeoJSONを作成する。

投影前の座標のLinearRing(先頭と末尾の点が等しい値となる[x, y]の配列、反時計回り)を作る

EPSG: 32645
Upper_Left: [382200.000, 2512500.000]
Lower_Left: [382200.000, 2279400.000]
Upper_Right: [610500.000, 2512500.000]
Lower_Right: [610500.000, 2279400.000]

というデータがあったとしたら、

反時計回りかつ先頭と末尾が同じ点となる配列を作ります。

[
  [382200.000, 2512500.000], 
  [382200.000, 2279400.000], 
  [610500.000, 2279400.000], 
  [610500.000, 2512500.000], 
  [382200.000, 2512500.000]
]

両極点を含む場合エラーとする

GeoJSONのPolygonは経度は-180 -> 180deg、緯度は-90 -> 90degの範囲であるべきである、とされています。
極点(北極点、南極点)を囲むようなポリゴンの場合、どこかで切らないことにはこの範囲に収まるポリゴンにすることができません。
一度polar sterographic projections(EPSG:3995かEPSG:3031)に投影して、切り開いてからEPSG:4326に投影することにします。[3]。

自己交差していたり潰れている場合エラーとする

衛星データは基本的には自己交差することはないし、潰れて線になってしまっていることもないはずなのでエラーとします。

LinearRingの各辺をN分割する

単純に四隅だけを指定したポリゴンを変換すると、本来曲線であるべき辺が直線になってしまいます。

投影変換前は四角形だったポリゴンも、地球は球なので変換すると各辺が歪んで曲線となる場合があります(表示サンプルは歪みが大きい極付近のポリゴン)。

before EPSG:3411

after EPSG:4326

しかし分割せずに四隅だけ変換してしまうと歪みまで考慮されません。

だめな例

辺を事前にN分割しておいて曲線を近似できるようにします。
分割数が多いほど元の曲線に近づきます。
ただ、あまり分割数を大きくしすぎると最終的に出力されるGeoJSONのサイズが大きくなるため、いい塩梅を探す必要があります。10以上がおすすめです(理由は後述)。

LinearRingをEPSG:4326に投影変換

変換にはproj4、Pythonならpyprojを使います。
あまりライブラリを使いたくないとはいえ、こればかりは仕方ないです。

投影変換したLinearRingが反時計回りか判定する

変換前のLinearRingを反時計回りとしたので、投影変換後も反時計回りとなるはずです[4]。
もし時計回りとなっていたら、それは投影されたポリゴンが180degをまたいでいると判断できます。

投影変換したLinearRingが反時計回りのとき

変換完了

bboxを[minlon, minlat, maxlon, maxlat]とし、geometryのcoordinatesに変換したLinearRingをセットすれば出力すべきGeoJSONとなります。

{
  "type": "Feature",
  "bbox": bbox,
  "geometry": {
  　　"type": "Polygon",
     "coordinates": [LinearRing]
  } 
}

投影変換したLinearRingが時計回りのとき

ポリゴンが180degをまたいでいると判断します。

経度方向の長さが180以上の辺は、180degをまたいでいると判定して180degとの交点を挿入する

GeoJSONは3.1.9. Antimeridian Cuttingで、180degをまたいでいるポリゴンは分割することを推奨(SHOULD)しています。
180degとの交点を見つけて複数のLinearRingに分割しましょう。

ただ、辺が180degと交わっているかは、端点からだけでは判定することができません。
[170, -10], [-170, -10]の2点は、長さが20([170, -10]->[-170, -10])で180degと交わっているかもしれないし、長さが340([-170, -10]->[170, -10])で裏側を回っているかもしれないためです。

そこで、この時点でLinearRingの各辺の経度方向の距離がワープしていない状態で180未満である、という条件を与えます。

例えば、
[170, -10]->[-170, -10]の辺であればワープしていない状態では[170, -10]->[190, -10]となり、経度方向の距離が180より小さいためOKです。
[-170, -10]->[170, -10]の辺であれば、経度方向の距離が180より大きいため、辺を分割して出直してこい、となります。

幸いにも曲線を近似するため変換前に辺を分割しているので、条件をクリアしている、はずです[5]。
よって、2点の経度方向の距離が180より大きい辺は180degをまたいでいる、と判定しましょう。

180degより左側と右側にポリゴンを分類する。

bboxの計算のためにポリゴンを選り分けます。

変換完了

GeoJSONのbboxは5.2. The Antimeridianによると、図形が180degをまたいでいるときは[177.0, -20.0, -178.0, -16.0]のようになります。
bbxoの定義が[minlon, minlat, maxlon, maxlat]でなく[leftlon, northlat, rightlon, southlat]ということです。

そんなわけで180degより左側と右側に選り分けたポリゴンたちから次のようにbboxを計算します。
bbox=[左側のxの最小値, yの最小値, 右側のxの最小値, yの最小値]

そしてgeometryはMultiPolygonとしましょう。

{
  "type": "Feature",
  "bbox": bbox,
  "geometry": {
  　　"type": "MultiPolygon",
     "coordinates": [左側のLinearRings, 右側のLinearRings]
  } 
}

出力例

ここまでのロジックをJavaScriptで実装したのがgeo4326です。

EPSG: 32645
Upper_Left: [382200.000, 2512500.000]
Lower_Left: [382200.000, 2279400.000]
Upper_Right: [610500.000, 2512500.000]
Lower_Right: [610500.000, 2279400.000]

というデータから、分割数10で作ったのがこちら。

{
  "type":"Feature",
  "bbox":[85.85296718933647,20.610041795245515,88.07596179098907,22.719775713801845],
  "geometry":{
    "type":"Polygon",
    "coordinates":[[[85.85296718933647,22.71566587084141],[85.85471144948144,22.505128271679137],[85.85643508166952,22.29458504609756],[85.85813822451726,22.08403623614523],[85.8598210145163,21.87348188415423],[85.86148358606195,21.662922032738383],[85.86312607148135,21.452356724791283],[85.86474860106108,21.24178600348436],[85.86635130307411,21.031209912264952],[85.86793430380634,20.820628494854294],[85.8694977275825,20.610041795245515],[86.08856914891048,20.611335735650073],[86.30765561712383,20.612352097028005],[86.52675351639061,20.613090808785984],[86.7458592298015,20.613551819609576],[86.96496913971066,20.613735097472322],[87.18407962807663,20.613640629641363],[87.40318707680342,20.6132684226797],[87.6222878680816,20.612618502445013],[87.84137838472925,20.611690914085163],[88.06045501053292,20.610485722030123],[88.06192160852325,20.82107736606614],[88.0634065712253,21.031663741215127],[88.06491001527168,21.242244803670587],[88.06643205913487,21.45282050989972],[88.06797282315154,21.663390816645528],[88.06953242954742,21.873955680928745],[88.0711110024627,22.084515060049903],[88.07270866797775,22.295068911591237],[88.07432555413966,22.5056171934187],[88.07596179098907,22.71615986368381],[87.8536827240295,22.717500978966182],[87.63138882771706,22.718533183697783],[87.4090839625065,22.71925640487522],[87.18677198995947,22.71967059134279],[86.96445677235506,22.719775713801845],[86.74214217229994,22.719571764815985],[86.51983205233836,22.719058758812512],[86.29753027456248,22.71823673207974],[86.07524070022231,22.717105742760463],[85.85296718933647,22.71566587084141]]]
    }
}

Landsat-8のMLT.TXTを見ると、計算したbboxとCORNER_UL_LAT_PRODUCT~CORNER_LR_LON_PRODUCTがだいたい同じ値になってます。

見つからなかった+勉強しようのモチベーションで作って見ましたが、これよりもシンプルな解法や既存のライブラリがあればそちらに鞍替えしたいのでぜひ教えて下さい。

脚注

  1. 衛星データを扱う際にデファクトになりつつあるrasterioは180degをまたぐケースの計算がテキトウで今回は使えません。 ↩︎

  2. 無償でも一部機能は使えるが、欲しい機能が無償かはわからなかったので深追いはせず。 ↩︎

  3. 両極を含んでいるような図形(正距方位図法で北極南極を含むような大きなポリゴンなど)はさすがにどうしようもなかったので諦めました。 ↩︎

  4. 変換しても回転方向が保持される、と思っていますが、すべての条件を試したわけではないので、もしかしたら例外があるかも。 ↩︎

  5. はずとしか言えないのが分割数によっては条件を満たせていないかもしれないからです。エイヤで10分割もすれば、ほとんどのケースで大丈夫だとは思うのですが、元の投影法によっては満たしていると必ずしも言い切れない......。上手くいかなかったら分割数を増やしてみてください。 ↩︎

Discussion