概要
地球は楕円体なので、地図に真円を描画しても中心から等しい距離の軌跡でありません(数学感)。先人達の記事を見てもよくわからなかったので自分でやってみました。
注意
ちゃちゃっと作ったものなのでミスがあるかもしれません。コードはCopilotで移植したものを使えるよう調整したものです。命名が雑だったりおかしかったりするのは気にしないでください。
メルカトル図法とか計算が面倒でよくわからないので、正距円筒図法(座標をそのまま投影するもの)を使用しています。別の方法でもたぶん大丈夫です。
もうちょっとちゃんと作ったら追記します。
内容
どんな感じになるのか？
地図と中心座標を用意して、1pxごとに距離を計算して色を塗りました。距離の計算は経度緯度の情報から「測地線航海算法」で距離を計算するを参考にしました。
描画部分を抜き取ったものは下の通りです。
var dist = GeodesicDistance.Dist([centerLat, centerLon], [imgLat, imgLon]) / 1000d;//距離(km)
if (double.IsNaN(dist))//緯度経度が一致するとNaNになる？
continue;
var level = dist / cor;//色段階
var co = level > 765 ? Color.FromArgb(a, 128, 128, 128) : level > 510 ? Color.FromArgb(a, 0, 0, (int)level - 510) : level > 255 ? Color.FromArgb(a, 0, (int)level - 255, 0) : Color.FromArgb(a, (int)level, 0, 0);
g.DrawLine(new Pen(new SolidBrush(co)), x, y, x + 1, y);//クソコード
(点の描画は幅1px長さ2pxの線を1pxぶん上書きしていく適当な実装)(aは透明度、corは距離を合わせるための値です。(255*cor)kmで色が切り替わります(たぶん))
無茶苦茶な計算量ですが、私の環境では1080x1080で2秒弱で描画できました(地図描画除く)。
なんか地震界隈の某有名ソフトは到達円内にグラデーションがありますがこんな感じですね
それはさておき当然ながらきれいな楕円でもないです。
こんな感じで複数の楕円でもうまくいかなそうです。
Vincenty法の利用
調べているとどうやらVincenty法の順解法を使うことでできるようです。[Python]始点の緯度経度と方位角と距離から、終点の緯度経度と方位角を求める
を参考にしました。
var psApprox = new List<Point>();
for (double d = 0; d <= 360; d += deltaD)
{
var pt = Vincenty.VincentyDirect(centerLat, centerLon, d, newDist, 1);
psApprox.Add(new Point((int)((pt.Value.lon - config.LonSta) * zoomW), (int)((config.LatEnd - pt.Value.lat) * zoomH)));
}
終点の緯度経度が分かればいいですね。deltaDを小さくすれば円に近くできます。
できました！
どちらも同じ楕円体モデルで高精度のため、うまくできています。
PS波描画への応用(2024/10/01追加)
走時表を使えば深さ・震央距離からPS波それぞれの到達時間が分かります。しかし深さと時間からPS波それぞれの震央距離が欲しいです。これは線形補間を使って近似的に求めることができることは先人達の記事で知っていたので利用させていただきます。
未到達だと(処理方法によっては)距離がマイナスになって未到達なのに到達していたり変なことになるので0kmとかにしときましょう。
JMA2001走時表(tjma2001.zipの方)を使用しました。震源は40.3°N,145.2°E,D=170kmとしました。
ばらつきはありますがだいたい20msで描画できます(どれくらい内部で最適化されているかわかりませんが、最初でも100ms行かないくらいです)。楕円だと違和感があるのでできればメルカトル図法を使いたいですね。
まとめ
Vincentyの式を利用することで簡単かつ正確に、高速に描画できました。また、これをもとに走時表を用いてPS波を描画できました。
距離単位がmだったりkmだったりするのでよく確認すること。(2敗)
汚いので使い物にならない気がしますが、ソースコードはここです。