地球上の2点間の距離を知りたいことがある。
知りたいのは、3次元空間の直線距離ではなく、地球表面にそった、そのかぎりでの最短距離だとする。 地球の形は球で近似し、楕円や凹凸は考えない。すると、最短距離は大円の弧の長さになる。
つぎの球面三角の公式が利用できる。 ここで、a, b, c は 3つの辺の長さで、半径が 1 の球の大円の弧の長さとしてあらわされている。 A, B, C はそれぞれの辺にむかいあう内角で、ラジアンを単位としてあらわされている。
cos a = cos b cos c + sin b sin c cos A
参考: 国立天文台, 2021: 『理科年表 2022』 丸善出版, p. 1149 (附 15)。
北極を 点A とし、点B と 点C の大円距離をもとめることにすれば、公式をそのままつかえる。
地球上の大円距離は、a に地球の半径をかければよい。 地球の大円の周がわかっていれば、地球の半径は、それを 2π でわればよい。 ここでは、地球の半径として、6371 km をとることにする。 これは、赤道半径よりも小さく、極半径よりも大きい。
これをつかって、2地点の緯度・経度 を 「度」単位であたえられて、 距離を km 単位で返す Python の 関数 daien_kyori() をつくった。
それをつかうプログラム例として、 熊谷 (熊谷地方気象台の位置) から、 『理科年表 2022』に収録されている世界の気温・降水量の観測地点 それぞれまでの大円距離を計算して書きだすものをつくってみた。
陸上の地点の海からの距離を知りたいことがある。
そのためには、海岸線の位置の情報が必要である。 海岸線データセットには、海岸線を折れ線でひくことを想定して、 海岸線がとおる地点の緯度・経度がサンプリングされて収録されている。 ここではひとまず、陸上の地点から海までの距離を、 陸上の地点から海岸線データセットを構成する各地点までの大円距離の最小値で近似することにする。 これは、海岸線の位置を連続の線でつないで最短距離を計算したばあいにくらべて、 海岸線データセットのサンプリング間隔 (つぎにしめすデータのばあいは 1～2 km) と 同程度の距離の誤差をよけいにふくんでいる。 そのことを承知で、概算としてつかうことにする。
ここでは、Cartopy で海岸線をひくためにもつかっている 「Natural Earth」の世界の地理データのうち「1:10m Physical Vectors」 https://www.naturalearthdata.com/downloads/10m-physical-vectors/ から、「Coastline version 4.1.0」ne_10m_coastlines.zip をつかう。 (この「1:10m」は、地図の縮尺 1千万分の1を意味する。 Natural Earth の線データのうちでは比較的くわしいものである。) ただし、Cartopy の中で地理データとしてつかうのではなく、 数値データとしてつかうことにする。
上記のzipアーカイブをダウンロードして、展開すると、ne_10m_coastlines.shp がある。 これは、ESRI社がきめた Shapefile という形式になっている。 これを、Linux上で、GDAL というパッケージに属する ogr2ogr というプログラムで gmt 形式 (GMT = Generic Mapping Tools というソフトウェアの入力になりうるテキストファイル形式) に変換し、 さらに、折れ線のはじめ・おわりの情報をあらわす行をとりのぞいて、 1行に1地点の経度と緯度がスペースくぎりで書いてあるだけのテキストファイルにした。 (折れ線の構造の情報はうしなわれ、単なる点の列のデータになっている。) そして、[地形 (標高) データを Pythonで読み、作図する (試作版)] で地図を表示した範囲と同じ、 東経 138.5～141.5度、北緯 34.5～38.5度の範囲にふくまれる地点だけを選択した。
【ここまでの Linux上での 操作はつぎのとおり。
ogr2ogr -f gmt -lco GEOMETRY=AS_WKT work1.gmt downloads/ne_10m_coastline.shp grep -v \# work1.gmt > work2.gmt grep -v \> work2.gmt > work3.gmt awk '($1 >= 136.5) && ($1 <= 141.5) && ($2 >= 34.5) && ($2 <= 38.5)' work3.gmt > kaigan.txt
なお、ogr2ogr で 「-f csv」とすれば CSV への変換がおこなわれるが、それによってつくられる CSV ファイルは ここでの目的に適したものにならなかった。】
【[2022-11-27 補足] ogr2ogr コマンドをふくむ GDAL パッケージは MS Windows でも動くことが知られているが、 増田はまだ これを Windows にインストールすることができておらず、 Linux が動いているパソコンにデータを持って行って作業している。 なお、grep と awk は Windows 用の UnitUtils というオープンソースのパッケージを入れれば 「コマンドプロンプト」で動く (ただし awk のコマンド名は gawk)。】
陸上の地点の例として、気象庁のAMeDASをふくむ地上観測の地点のうち、 関東7都県と福島、新潟、長野、静岡、山梨県にある地点のリスト ([気象観測データによる風の空間分布の作図 (Python による試作)] の「第2回試作」のところでつくった amedas_stn_near_kanto.csv )をつかうことにする。 これにふくまれた地点のうちで、上記の緯度経度範囲の外のものは、 プログラム内で判断して、出力から除外することにした。 なお、上記の緯度経度範囲の境界に近い地点については、 結果が海からの距離の概算値としては不適切な値になっている可能性がある。
海からの距離の概算のための、東アジアの大陸本土の海岸線の経度・緯度のデータを作成した。
「Natural Earth」の世界の地理データのうち「1:50m Physical Vectors」 https://www.naturalearthdata.com/downloads/50m-physical-vectors/ から、「Coastline version 4.0.0」ne_50m_coastlines.zip をつかう。
zipアーカイブを展開し、MANDARA 10 に読みこんでマップエディタで表示させ、 ユーラシア・アフリカの大陸本土の海岸線が ne_50m_coastline.shp.1388 であることが わかったので、MANDARA 10の属性データ編集で、このオブジェクトだけ残してほかを消し、 ne_50m_coastline_eurasaf.shp (および .shx, .dbf, .prj) として保存した。 このシェープファイル から、上と同じ手順で、 東経 60～150度、北緯 0～60度の範囲にふくまれる地点だけを選択し、 結果を kaigan_eastasia.txt とした。
同じne_50m_coastlines.zip から、MANDARA 10 で 日本の主要4島の海岸線を取り出した。 オブジェクト名の最後の部分としてでてくる番号 ([2022-11-01] 「オブジェクト番号」とはくいちがっていることがわかった。])は、 109 (北海道)、110 (本州)、111 (四国)、112 (九州) である。 そこから上と同じ手順で、海岸線の経度・緯度だけのデータファイルにした。
[2022-11-02 暫定追加] 同じne_50m_coastlines.zip から、上記以外の日本の島のうち、 (当面つかいたい気象観測点のつごうで) 沖縄島 (オブジェクト名の最後の部分としてでてくる番号 707)、 佐渡島 (717)、対馬 (南部が727, 北部が728) を取り出し、 上記と同じ形のファイルにした。
[2022-11-07 暫定追加] 同じne_50m_coastlines.zip から、中国の島のうち、 (当面つかいたい気象観測点のつごうで) 海南 (Hainan) 島 (オブジェクト名の最後の部分としてでてくる番号 744) を取り出し、 上記と同じ形 (くぎり文字の行をふくまないほう) のファイルにした。