🌏

ラスターデータからCloud Optimized GeoTIFF（COG）を生成する方法

2023/08/24に公開

GDAL

はじめに

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

Cloud Optimized GeoTIFF（COG）の特徴

HTTP-Rang Requestsによって大きな位置情報データのバイト列の一部分をリクエストすることで、一部の領域だけのデータを取得する仕組み
ラスターデータの標準であるGeoTIFF形式をRang Requestsを実現するために拡張したファイル形式
GeoTIFFのピラミッド画像※という仕組みをベースに、部分データの配信を実現
JAXAがCOG形式で衛星画像を配信
※オリジナルの解像度のほか、低解像度の画像を内部的に保持する仕組み。

（引用元：『現場のプロがわかりやすく教える位置情報エンジニア養成講座』）

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/

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

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のファイルパスを設定します。
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

COGの生成

GeoTIFFからCloud Optimized GeoTIFF（COG）への変換にはOSGeo4Wのgdal_translateを使用します。
まず、使用するGeoTIFFがCOGではないことを確認します。
GeoTIFFがCOGかどうかを調べるpythonスクリプト（validate_cloud_optimized_geotiff.py）が下記にて提供されています。

OSGeo4W Shellを起動して、下記のコマンドを実行してGeoTIFFがCOGかどうかを調べます。

python validate_cloud_optimized_geotiff.py lp_orthophoto_add_srs.tif

COGではないので NOT a valid と表示されます。

The following warnings were found:
 - The file is greater than 512xH or Wx512, it is recommended to include internal overviews

lp_orthophoto_add_srs.tif is NOT a valid cloud optimized GeoTIFF.
The following errors were found:
 - The file is greater than 512xH or Wx512, but is not tiled

下記のコマンドを実行してCOGを作成します。

gdal_translate -a_srs EPSG:6676 -of COG lp_orthophoto_add_srs.tif lp_orthophoto_add_srs_cog.tif -co TILING_SCHEME=GoogleMapsCompatible -co COMPRESS=DEFLATE

再度、下記のコマンドを実行してGeoTIFFがCOGかどうかを調べます。

python validate_cloud_optimized_geotiff.py lp_orthophoto_add_srs_cog.tif

下記のように表示されればCOGです。

lp_orthophoto_add_srs_cog.tif is a valid cloud optimized GeoTIFF

The size of all IFD headers is 34246 bytes

下記のコマンドを実行してCOGのメタ情報を表示します。

gdalinfo lp_orthophoto_add_srs_cog.tif

上記のCOG作成時に出力ファイルの座標参照系をEPSG:6676に設定しましたが、TILING_SCHEME=GoogleMapsCompatibleを追加することで、最終的な出力ファイルの座標参照系はWeb Mercator（EPSG:3857）になるようです。

Driver: GTiff/GeoTIFF
Files: lp_orthophoto_add_srs_cog.tif
Size is 12032, 16384
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06ﾂｰS and 85.06ﾂｰN."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (15480179.842370461672544,4179882.454746577888727)
Pixel Size = (0.298582141738970,-0.298582141738970)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  LAYOUT=COG
Tiling Scheme:
  NAME=GoogleMapsCompatible
  ZOOM_LEVEL=19
Corner Coordinates:
Upper Left  (15480179.842, 4179882.455) (139d 3'38.96"E, 35d 7' 3.58"N)
Lower Left  (15480179.842, 4174990.485) (139d 3'38.96"E, 35d 4'54.15"N)
Upper Right (15483772.383, 4179882.455) (139d 5'35.14"E, 35d 7' 3.58"N)
Lower Right (15483772.383, 4174990.485) (139d 5'35.14"E, 35d 4'54.15"N)
Center      (15481976.113, 4177436.470) (139d 4'37.05"E, 35d 5'58.87"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 6016x8192, 3008x4096, 1504x2048, 752x1024, 376x512, 188x256

おまけ

上記の方法で生成したCOGをMapLibre GL JSで表示するデモサイトを下記のGitHubにて公開しています。