2025-01-20 海岸線データを作る試み

天気図の海岸線データを作りたいのだが、知財的にクリアされたデータでやりたい。

「coastline data」とかで雑にぐぐると Natural Earth のパブリックドメインデータがみつかる。
https://www.naturalearthdata.com/downloads/50m-physical-vectors/50m-coastline/

とにかくダウンロードしてみる。
https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/physical/ne_50m_coastline.zip

中はshapefileである。PRJファイルをみると普通に経緯度座標であるらしい。

ne_50m_coastline.prj

GEOGCS["GCS_WGS_1984",
 DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
 PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

shapefileを何らかのテキストファイル、たとえばGeoJSONに変換したいと思ってぐぐる。
https://zenn.dev/yuichiyazaki/articles/b70da4bb63ea02
GDALのogr2ogrを使うらしい。

# apt install gdal-bin
$ ogr2ogr -f GeoJSON ne_50m_coastline.geojson ne_50m_coastline.shp
$ wc ne_50m_coastline.geojson 
   1434  271663 2769871 ne_50m_coastline.geojson

中を見るとまっとうな GeoJSON ぽい。属性は min_zoom というのがあって、0.0〜4.0の範囲である。天気図としてはズームレベル6くらいを想定しているので、小さい島を捨てる処理は要らないかな。

とりあえずマトモに読めることと、構造を確認する。マトモに読めるとは、JSONであるからメモリを食いすぎてサーバが落ちる事態がありえるので先にあたりをつける趣旨。構造とは、GeoJSONには多様な図形種類がありえて、そのどれを使っているかである。

testgjlns.rb

#!/usr/bin/ruby

require 'json'
j=JSON.parse(File.read(ARGV.shift))
n=0
raise unless j['type']=='FeatureCollection'
j['features'].each{|f|
 raise f['type'] unless f['type']=='Feature'
 g = f['geometry']
 mz=f['properties']['min_zoom']
 case g['type']
 when 'LineString' then
   cs=g['coordinates'].size
 when 'MultiLineString' then
   cs=g['coordinates'].sum{|ls| ls.size}
 else
   raise g['type']
 end
 printf "%3u\t%3.1f\t%3u\n", n, mz, cs
 n+=1
}

最初はもうちょっとかんたんなプログラムで、主に各オブジェクトの type メンバを assertion（想定外なら raise で落とす）しながら落ちないようにプログラムを育てていった結果が上記。何故かは知らんが単純LineStringとMultiLineStringが混在していた。
出力。

$ ruby testgjlns.rb ne_50m_coastline.geojson | head
  0	 1.5	 56
  1	 4.0	  8
  2	 4.0	  8
  3	 3.0	  9
  4	 4.0	  7
  5	 4.0	  5
  6	  0.0	1154
  7	  0.5	 46
  8	 1.0	 65
  9	 1.0	101

CやFortranで読みやすいように整形出力する。

lnsgjtx.rb

#!/usr/bin/ruby

require 'json'
j=JSON.parse(File.read(ARGV.shift))

n=0
raise unless j['type']=='FeatureCollection'
j['features'].each{|f|
 raise f['type'] unless f['type']=='Feature'
 g = f['geometry']
 case g['type']
 when 'LineString' then
  n+=1
 when 'MultiLineString' then
  n+=g['coordinates'].size
 else
  raise g['type']
 end
}
printf("LINES %10u\n", n)

def show ls, i
  printf("%10u %10u\n", i, ls.size)
  ls.each{|x, y|
    printf("%+8.3f %+7.3f\n", x, y)
  }
end

i=0
j['features'].each{|f|
 g = f['geometry']
 case g['type']
 when 'LineString' then
  show(g['coordinates'], i)
  i+=1
 when 'MultiLineString' then
  g['coordinates'].each{|ls| show(ls, i); i+=1}
 end
}

次はこれを図化して確認する。プログラムを書こう。
https://github.com/etoyoda/grib2png/issues/48

素朴に全球プロットしてみるとまっとうに見える。

メルカトル図法に変更