地球の球体を考慮した２点間の最短ライン

[キタ]

QGISで世界地図を表示し、２点を結ぶラインを作成したとき（例えば東京とパリ）、地球の球体を考慮したラインの軌跡（弧になると思います）を描くにはどうしたら良いでしょうか？
座標系はWEBメルカトルのEPSG3857で描きたいと思っています。
よろしくお願いします。

キタ

[ありた]

こんにちは

地球表面上の最短ルート（いわゆる大圏コース）を求める方法として、ぱっと思い浮かぶのは

始点を中心とした正距正方位図法上で終点の値を得て、ルートを作成する方法です。

ただし始点と終点のみのラインデータだと、たとえ正距正方位図法のジオメトリデータでも

（真球メルカトル等）他の投影法で投影した際は、始点と終点の情報から直線で表示されるため

正距正方位図法上のルートで、適宜折れ点を作成しておく必要があります。

ためしに QGIS のフィールド計算機の関数エディタで great_circle_route 関数を作ってみました。

（折れ点は約100km間隔で入るようにしています。）

from qgis.core import *

from qgis.gui import *

import pyproj

@qgsfunction(args='auto', group='Custom')

def great_circle_route(lat1, lon1, lat2, lon2, feature, parent):

wgs84 = pyproj.Proj('+init=EPSG:4326')

azeq = pyproj.Proj('+proj=aeqd +lat_0=%f +lon_0=%f +x_0=0 +y_0=0 '

'+datum=WGS84 +units=m +no_defs' % (lat1, lon1))

x2, y2 = pyproj.transform(wgs84, azeq, lon2, lat2)

dist = (x2 ** 2 + y2 ** 2) ** 0.5

nodes = []

edges = int(dist // 100000) + 1

for i in range(edges):

lon, lat = pyproj.transform(azeq, wgs84, i * x2 / edges, i * y2 / edges)

nodes.append(QgsPoint(lon, lat))

nodes.append(QgsPoint(lon2, lat2))

return QgsGeometry.fromPolyline(nodes)

たとえば EPSG:4326 のポリラインのレイヤで、適当な地物をひとつ作成して保存します。

フィールド計算機で既存フィールドの更新でジオメトリを選択し、次のような式で更新します。

great_circle_route(35.7653, 140.3856, 49.0097, 2.5486)

もちろん属性値が整備されていれば、座標は直接入力ではなく次のようにもできます。

great_circle_route("lat1", "lon1", "lat2", "lon2")

ジオメトリデータとしては EPSG:4326 に戻していますが、データとして EPSG:3857 で

ほしい場合には、フィールド計算機の transform() をかませてください。

作ってみたものの、プラグイン化していないので QGIS では使いにくいです。

ただ、始点終点ペアが与えられていて、ローカルで Python (+PyProj) が使える場合や、

WebGIS で始点終点を選んだら、大圏コースをプロットするなどは比較的容易そうですね。

（ Proj や PyProj は最近大きなアップデートがあったため、インストール済みバージョンに注意）

[adachi]

キタ様

考え方としてはありた様と同じものかと思いますが、なんとかGUIでの操作メインでやる方法がないか探してみました。

以下を参考にしています。

https://gis.stackexchange.com/questions/133026/great-circles-in-qgis-and-export-in-3857-webmap

1. 任意の点を中心とした正距方位図のカスタム投影法を作成しておく（今回は東京駅付近の座標としています）

下記式で、緯度経度を中心点に応じて変更してください。

+proj=aeqd +lat_0=中心点緯度 +lon_0=中心点経度 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

2. プロジェクト座標系をEPSG:4326として、東京とパリのポイントを追加（見た目はかなりずれてますがちゃんと乗っかっています）

3. プロジェクト座標系を1で作成したカスタム座標系に変更し、２点間を結ぶラインデータを追加

4. できたラインに対して「頂点の高密度化（個数ベース）」を使用し、+99個程度で処理（ここの個数は任意で良いと思います）

5. プロジェクト座標系をEPSG:3857に変更

[adachi]

追記：

説明ではプロジェクト座標系しか書きませんでしたが、レイヤの座標系は

ポイントデータ：4326

ラインデータ：カスタム座標系

で作成しています。

[福岡]

キタ　様

おはようございます。

以前「測地線は取り扱えるの？」という質問があり、Shape Tools というプラグインが使えるように思うと投稿したことがあります。

測地線は取り扱えるの？

https://groups.google.com/forum/#!topic/qgisshitumon01/ZN_eb0h11js

確認してみたところQGIS3.xでもプラグインにありましたので、こちらを使うという方法もあるかと思います。

[キタ]

ありたさん
adachiさん
ありがとうございます。

福岡さん
確か同じような質問見たことあるなあと思ってはいたんですが、言葉がわからなくて探せませんでした。
この掲示板もおかげさまで質問がたくさん来ているのですが、探すのも大変になってきてますね。
細かな説明、ありがとうございました。
プラグインも含め、試してみようと思います。

キタ

[キタ]

adachiさんの方法でできました。

ありがとうございました。

キタ