[memo] Google マイマップで円を描きたい

はじめに

Google マイマップで円とか四角とかを描きたかったのですが、フリーハンドではうまく描けなかったので、ちょこっと調べたことのメモです。
あまり精度は高くなくてもいいことが前提です。

Google マイマップの作り方

1. Google ドライブで「新規」→「その他」→「Google マイマップ」でマップを作成
2. レイヤに、ポイントを追加したり、図形を書き込んだり、ルートを作ったりできる。
3. レイヤをエクスポートすることもできる（KMZ ファイルと KML ファイル、CSV ファイルなど）
   • KMZ は ZIP 形式で圧縮されたファイル。KML や画像ファイルなどがまとめてあるみたい。
   • KML は XML 形式の 1 ファイル
4. 新規レイヤにインポートすることで、ポイントとか図形を取り込むことができる（KMZ/KML/CSV など）

とりあえず、適当なポイントや図形を描いて、KML 形式でエクスポートしてみました。
いろいろなスタイルなどが入っているけど、いらなそうなものを削っていって、図形一つでインポートができる形にしてみました。
（もしかしたら、さらに消せるものとかあるかもしれません。逆に仕様的には必須の物も Google マイマップでは動くからといって消しちゃっているかもしれません）

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document>
<name>テスト</name>
<Placemark>
<name>テスト</name>
<Polygon>
<outerBoundaryIs>
<LinearRing>
<tessellate>1</tessellate>
<coordinates>
 138.5000000,28.5000000,0
 141.5000000,28.5000000,0
 141.5000000,31.5000000,0
 138.5000000,31.5000000,0
 138.5000000,28.5000000,0
</coordinates>
</LinearRing>
</outerBoundaryIs>
</Polygon>
</Placemark>
</Document>
</kml>

上記は、適当な座標での四角を描いてみたものです。
coordinates 要素に、図形の頂点の緯度と経度を指定すればよいだけですね。
見たところ、最初と最後の点を同じにして、手動で線を閉じているしている感じですね。

また、複数の図形を KML ファイルに入れたい場合には、Placemark 要素（ <Placemark> ～ </Placemark> ）を複数並べる感じみたいです。

それと、おそらく正確な円を描くことはできないみたいで、正多角形を描いて円を近似することになるようです。

特定の緯度での距離を求める

あとは、図形の頂点を求めれば円でも四角でも好きに作れそうです。

問題は、特定の緯度での距離⇔度数の変換をする必要があることです。

少し調べてみると、かなり複雑な計算式があるようです。

【参考】 国土地理院の測量計算サイト の 経緯度を用いた2地点間の測地線長、方位角を求める計算

今回の内容的にそこまで精度が高くなくてもよいので、もうちょっと扱いやすい簡易的な式がないかと調べてみたら、「ヒュベニの公式」というのがあるみたいなので、これを使ってみました。

2 点間の距離を $D$ とすると

\begin{align}
D &= \sqrt{ ( D_y \cdot M )^2 + ( D_x \cdot N \cdot \cos{ P } )^2 }\\
M &= \frac{ R_x ( 1 - E^2) }{ W^3 }\\
N &= \frac{ R_x }{ W }\\
W &= \sqrt{ 1 - E^2 sin^2{P} }\\
E &= \sqrt{ \frac{ R_x^2 - R_y^2 }{ R_x2 } }\\
\end{align}

• $D_y$ : 2 点間の緯度の差（ラジアン単位）
• $D_x$ : 2 点間の経度の差（ラジアン単位）
• $P$ : 2 点の緯度の平均
• $M$ : 子午線曲率半径
• $N$ : 卯酉線曲率半径
• $E$ : 離心率
• $R_x$ : 長半径（赤道半径） WGS84 では 6 378 137.000 メートル
• $R_y$ : 短半径（極半径） WGS84 では 6 356 752.314 245 メートル

こんな感じのようです。$D$, $D_y$, $D_x$, $P$ 以外は定数ですね。

今回知りたいのは、2 点間の距離ではなく、1 点から特定の距離離れた緯経度が知りたいので、この式を変形してみます。

まずは子午線方向について。

ある点の緯度を $Lat$ 、そこから東西方向に $D$ メートル離れた地点の経度差 $D_x$ を求めたいとします。
同じ緯度の点なので $D_y = 0$ ですね。
ここでは $D_x$ が正のケースを考えてみます。

\begin{align}
D &= \sqrt{ ( D_y \cdot M )^2 + ( D_x \cdot N \cdot \cos{ P } )^2 }\\
  &= \sqrt{ ( D_x \cdot N \cdot \cos{ P } )^2 }\\
  &= D_x \cdot N \cdot \cos{ P }\\
D_x &= \frac{ D }{ N \cdot \cos{ P } }\\
\end{align}

次に、卯酉線方向について（ちなみに、十二支を方角に当てはめると、子（ねずみ）が北、午（うし）が南、卯（うさぎ）が東、酉（とり）が西になるので、卯酉線は東西方向のことです。当然、子午線は北極と南極をつなぐ南北の線）。

ある点の緯度を $Lat$ 、そこから南北方向に $( D / 2)$ メートルづつ離れた 2 つの地点の緯度差 $D_y$ を求めたいとします。
この緯度差 $D_y$ の半分が、ある点 $Lat$ との緯度差だと扱うことにします。
南北の2 点は同じ経度の点なので $D_x = 0$ ですね。
ここでも $D_y$ が正のケースを考えてみます。

\begin{align}
D &= \sqrt{ ( D_y \cdot M )^2 + ( D_x \cdot N \cdot \cos{ P } )^2 }\\
  &= \sqrt{ ( D_y \cdot M )^2 }\\
  &= D_y \cdot M\\
D_y &= \frac{ D }{ M }\\
\end{align}

ここで $D$ は、本来求めたい距離（ $r$ とします）の倍になるのですが（ある点から、例えば北に 1000 メートル、あるいは南に 1000 メートルの点の緯度を求める場合には $D = 2r = 2000$ となる）、最終的に求めた $D_y$ を $2$ で割るので、結局は、$ D_y \div 2 = \frac{ r * 2 }{ M } \div 2 = \frac{ r }{ M }$ となり、距離を $M$ で割ればよいだけですね。
本来、緯度によって同じ 1 度でも距離が変わるはずなんですが、ヒュベニの公式ではそのあたりは簡易的になっている感じなんですね。

緯度をラジアンではなく度で扱うようにしたら、以下のような感じでいいのかな（ Python の例。 lat は度、 d はメートル）。

def calcDxDy(lat: float, d: float) -> Tuple[float, float]:
    rx = 6378137.000
    ry = 6356752.314245
    e = math.sqrt( (rx**2 - ry**2) / rx**2 )
    w = math.sqrt( 1 - e**2 * math.sin(math.radians(lat))**2 )
    n = rx / w
    m = rx * (1 - e**2) / w**3

    dx = math.degrees(d / n / math.cos(math.radians(lat)))
    dy = math.degrees(d / m)
    return (dx, dy)

# 中心の緯度、経度 単位は度
(lat, lon) = ( 30.0000, 140.0000 )
# 円の半径（多角形だけど）単位はメートル
len = 2000
# 正 N 角形
N = 24

(dx, dy) = calcDxDy(lat, len)
coordinates = ""
for i in range(0, N+1):
    deg = i * 360 / N
    x = lon + math.sin(deg * math.pi / 180) * dx
    y = lat - math.cos(deg * math.pi / 180) * dy
    coordinates += " {x:.7f},{y:.7f},0\n".format(x=x, y=y)

参考にしたページ