Zenn
🌏

ラスターデータからPMTilesを生成する方法

2023/06/26に公開
GIS
GDAL
geotiff
PMTiles
mbtiles
tech

はじめに

GIS界隈ではCloud Optimizedが重要なキーワードになっています。Cloud Optimizedとは、特別なサーバ実装を必要とせず、大きな位置情報データの一部分を配信することを可能とするファイル形式の総称です。本記事は、Cloud Optimized形式のうち、PMTiles（ラスタータイル）の生成方法に焦点を当てた記事になります。具体には、静岡県（VIRTUAL SHIZUOKA）がG空間情報センターを通じてオープンデータとして公開している、オルソ画像データ（GeoTIFF）からPMTilesを生成する方法について説明します。

PMTilesの特徴

  • 単一ファイル形式のタイルデータ
  • ラスタータイル、ベクトルタイルともに対応
  • 地球規模のタイルマップの保存
  • S3 のようなストレージとHTTP Range Requestsのみを使用

https://protomaps.com/docs/pmtiles

PC環境

OS：Windows 10 Pro
CPU：Intel Core i7-9700
メモリ：32GB

前提条件

  • OSGeo4W（GDAL）がインストール済みであること。
  • GDALのバージョンは、GDAL 3.6.3（released 2023/03/07）です。

https://trac.osgeo.org/osgeo4w/

  • rio-mbtilesがインストール済みであること。
    （GeoTIFFからMBTilesを生成するツールです。）

https://github.com/mapbox/rio-mbtiles

  • pipでインストールできます。
pip install rio-mbtiles
  • go-pmtilesをダウンロード済みであること。
  • go-pmtiles_1.8.0_Windows_x86_64.zip
    （MBTilesからPMTilesを生成するツールです。）

https://github.com/protomaps/go-pmtiles/releases

  • G空間情報センターからオルソ画像データを取得済みであること。
  • 本記事では、LPデータ オルソ画像データを使用します。

https://www.geospatial.jp/ckan/dataset/shizuoka-2019-pointcloud

GeoTIFFのマージ

  • GeoTIFFのマージには、OSGeo4Wのgdal_merge.pyを使用します。
  • GeoTIFFは、下記のとおり、lp_orthophotoフォルダ内に格納されているとします。
lp_orthophoto
│  08NF2268.tif
│  08NF2269.tif
│  08NF2278.tif
│  08NF2279.tif
│  08NF2288.tif
│  08NF2289.tif
│  08NF2298.tif
│  08NF2299.tif
（省略）
  • 下記のbatファイルにGeoTIFFが存在するディレクトリや出力するGeoTIFFのファイルパスを設定します。
  • OSGeo4W Shellを開いて、その中でバッチファイル（gdal_merge.bat）を実行します。
gdal_merge.bat
@echo off
setlocal enabledelayedexpansion

REM gdal_merge.pyへのパスを指定します。
set GDAL_MERGE_PATH="C:\OSGeo4W\apps\Python39\Scripts\gdal_merge.py"

REM 出力するGeoTIFFの名前を設定します。
set OUTPUT_FILE=D:\GIS\geospatial\VIRTUAL_SHIZUOKA\lp_orthophoto.tif

REM 処理するファイルの拡張子を設定します。
set EXTENSION=*.tif

REM GeoTIFFファイルが存在するディレクトリを設定します。
set DIRECTORY=D:\GIS\geospatial\VIRTUAL_SHIZUOKA\lp_orthophoto

REM 指定したディレクトリ内の全ての.tifファイルを一覧表示します。
dir %DIRECTORY%\%EXTENSION%

REM 指定したディレクトリ内の全ての.tifファイルを集めます。
for /R %DIRECTORY% %%f in (%EXTENSION%) do (
    set "files=!files! %%f"
)

REM gdal_merge.pyを使用して、全ての.tifファイルを一つの.tifファイルにマージします。
python %GDAL_MERGE_PATH% -o %OUTPUT_FILE% !files!

座標参照系の付与

  • G空間情報センターから取得した、オルソ画像データには座標参照系が設定されていないため、GDALのgdalwarpを使用して、座標参照系を付与します。
  • なお、付与する座標参照系は、日本測地系2011／平面直角座標系第8系（EPSG：6676）になります。
  • OSGeo4W Shellを起動して、下記のコマンドを実行して座標参照系を付与します。
gdalwarp -t_srs EPSG:6676 lp_orthophoto.tif lp_orthophoto_add_srs.tif
  • また、下記のコマンドを実行して、座標参照系を付与されているかGeoTIFFのメタ情報を確認します。
gdalinfo lp_orthophoto_add_srs.tif
  • メタ情報の出力結果は下記のとおりになります。
Driver: GTiff/GeoTIFF
Files: lp_orthophoto_add_srs.tif
Size is 14000, 19500
Coordinate System is:
PROJCRS["JGD2011 / Japan Plane Rectangular CS VIII",
    BASEGEOGCRS["JGD2011",
        DATUM["Japanese Geodetic Datum 2011",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",6668]],
    CONVERSION["Japan Plane Rectangular CS zone VIII",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",36,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",138.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
        AREA["Japan - onshore - Honshu between approximately 137ﾂｰ45'E and 139ﾂｰE - Niigata-ken; Nagano-ken; Yamanashi-ken; Shizuoka-ken."],
        BBOX[34.54,137.32,38.58,139.91]],
    ID["EPSG",6676]]
Data axis to CRS axis mapping: 2,1
Origin = (51200.000000000000000,-97800.000000000000000)
Pixel Size = (0.200000000000000,-0.200000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   51200.000,  -97800.000) (139d 3'42.18"E, 35d 7' 1.70"N)
Lower Left  (   51200.000, -101700.000) (139d 3'41.31"E, 35d 4'55.14"N)
Upper Right (   54000.000,  -97800.000) (139d 5'32.76"E, 35d 7' 1.17"N)
Lower Right (   54000.000, -101700.000) (139d 5'31.84"E, 35d 4'54.62"N)
Center      (   52600.000,  -99750.000) (139d 4'37.02"E, 35d 5'58.16"N)
Band 1 Block=14000x1 Type=Byte, ColorInterp=Red
Band 2 Block=14000x1 Type=Byte, ColorInterp=Green
Band 3 Block=14000x1 Type=Byte, ColorInterp=Blue

MBTilesの生成

  • GeoTIFFからMBTilesの生成にはrio-mbtilesを使用します。
  • コマンドプロンプトを起動して、以下のコマンドを実行します。
  • 出力されるタイルのフォーマットは、PNG形式を指定しています。
  • 地図のズームレベル0から18までの地図タイルを作成します。
  • タイルサイズは、512x512ピクセルのタイルを指定しています。
rio mbtiles lp_orthophoto_add_srs.tif lp_orthophoto_add_srs.mbtiles --format PNG --zoom-levels 0..18 --tile-size 512 --resampling bilinear

PMTilesの生成

  • MBTilesからPMTilesの生成にはgo-pmtilesを使用します。
  • ダウンロード済みのgo-pmtiles_1.8.0_Windows_x86_64.zipの中からpmtiles.exeをコピーして、MBTilesと同じフォルダに格納します。
  • コマンドプロンプトを起動して、以下のコマンドを実行します。
pmtiles convert lp_orthophoto_add_srs.mbtiles lp_orthophoto_add_srs.pmtiles

生成したPMTilesの確認

  • 生成したPMTilesの確認には、PMTiles Viewerが便利です。
  • PMTilesをドラッグ&ドロップすると、画像が地図上に表示されます。

https://protomaps.github.io/PMTiles/

おまけ

  • 上記の方法で生成したPMTilesとMapLibre GL JSで作成したデモサイトを下記のGitHubにて公開しています。

https://github.com/shi-works/Raster-PMTiles-MapLibre-GL-JS

Discussion