南極のウェブ地図を試作する

[hfu]

データソース

https://x.com/gsi_chiriin/status/1780870248318837120?s=46&t=f8vcjeJOLjOap3CynrX6Cw →https://www.gsi.go.jp/antarctic/download_index_03.html

• 提供されているのは、印刷図画像
• 59MB 程度の PDF ファイル
• 扱いは、国土地理院技術資料（通常、国土地理院コンテンツ利用規約適用 - CC BY）
• 見た感じ、巨大な地図画像である

ツール

Dezoomify https://ourdigitalworld.net/2022/03/30/broken-tiles-a-retroconversion-project/?amp=1

普通に GDAL で捌く。

[hfu]

gdalinfo でチェック

Driver: PDF/Geospatial PDF
Files: R5_1000万分1南極大陸.pdf
Size is 9933, 7016
Metadata:
  CREATION_DATE=D:20240322170412+10'00'
  CREATOR=Adobe Illustrator 27.8 (Windows)
  DPI=300
  PRODUCER=Adobe Photoshop for Windows -- Image Conversion Plug-in
  TITLE=R5_1000万分1南極大陸
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 7016.0)
Upper Right ( 9933.0,    0.0)
Lower Right ( 9933.0, 7016.0)
Center      ( 4966.5, 3508.0)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  Overviews: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  Overviews: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Blue
  Overviews: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110
Band 4 Block=1024x1024 Type=Byte, ColorInterp=Alpha
  Overviews: 4966x3508, 2483x1754, 1242x877, 621x438, 310x219, 155x110

[hfu]

GCP を設定する

ChatGPT

ウェブメルカトルのタイルの z=0 の範囲は、世界全体をカバーする単一のタイルです。このタイルの範囲をEPSG:3857で表現するためには、その範囲の左下隅と右上隅の座標を求める必要があります。

ウェブメルカトル投影では、地球を球体ではなく楕円体であるWGS84楕円体に投影します。この楕円体上でのz=0のタイルは、緯度経度座標で左下隅が (-180, -85.05112877980659)、右上隅が (180, 85.05112877980659) になります。これらの座標をEPSG:3857に変換することで、ウェブメルカトルの z=0 のタイルの extent を示すことができます。

座標の変換には、GDAL（Geospatial Data Abstraction Library）のコマンドラインツールである gdaltransform を使用します。以下は、左下隅と右上隅の座標を変換するコマンドです。

左下隅の座標を変換する:

echo "-180 -85.05112877980659" | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857

右上隅の座標を変換する:

echo "180 85.05112877980659" | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857

これらのコマンドを実行すると、左下隅の座標が (-20037508.342789 -20037508.342789) に、右上隅の座標が (20037508.342789 20037508.342789) になります。

したがって、EPSG:3857で表現されるウェブメルカトルの z=0 のタイルの extent は、左下隅が (-20037508.342789, -20037508.342789) であり、右上隅が (20037508.342789, 20037508.342789) です。

[hfu]

紙の上で初等的な計算をして、このような形で変換を試しています。 上記の 20037508.342789 を M と呼び、余裕を持って M=20030000 としています。

gdal_translate -of MBTiles -co "ZOOM_LEVEL_STRATEGY=UPPER" -a_srs EPSG:3857 -gcp 0 0 -20030000 14147839 -gcp 9933 7016 20030000 -14147839 -gcp 9933 0 20030000 14147839 -gcp 0 7016 -20030000 -14147839 R5_1000万分1南極大陸.pdf output.mbtiles

[hfu]

[hfu]

First quality check

Not very bad.

[hfu]

修正

gdal_translate -of MBTiles -co BLOCKSIZE=512 -co TILE_FORMAT=WEBP -co ZOOM_LEVEL_STRATEGY=UPPER -co RESAMPLING=LANCZOS -a_srs EPSG:3857 -gcp 0 0 -20030000 14147839 -gcp 9933 7016 20030000 -14147839 -gcp 9933 0 20030000 14147839 -gcp 0 7016 -20030000 -14147839 ~/Desktop/R5_1000万分1南極大陸.pdf output.mbtiles; gdaladdo -r lanczos output.mbtiles 2 4 8 16 

minzoom = 1, maxzoom = 5 のタイルセットができる。

[hfu]

ファイルシステム展開

./mb-util output.mbtiles zxy

1028タイル、33.6MB

[hfu]

リポジトリ作成

https://github.com/optgeo/tma

[hfu]

z=0 追加

gdal_translate -of MBTiles -co BLOCKSIZE=512 -co TILE_FORMAT=WEBP -co ZOOM_LEVEL_STRATEGY=UPPER -co RESAMPLING=LANCZOS -a_srs EPSG:3857 -gcp 0 0 -20030000 14147839 -gcp 9933 7016 20030000 -14147839 -gcp 9933 0 20030000 14147839 -gcp 0 7016 -20030000 -14147839 ~/Desktop/R5_1000万分1南極大陸.pdf output.mbtiles; gdaladdo -r lanczos output.mbtiles 2 4 8 16 32

mb-util --image_format=webp output.mbtiles zxy

1029タイル、33.7MB

[hfu]

Observable サイト作成

https://observablehq.com/d/d980baa9701dd18c

[hfu]

完成したのでクローズします。