Go to Qiita Advent Calendar 2024 Top
LoginSignup
4
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@shi-works(shi works)

Cloud Optimized GeoTIFF（COG）をMapLibre GL JSで表示してみる

Last updated at Posted at 2023-06-15

はじめに

「Cloud Optimizedとは、特別なサーバ実装を必要とせず、大きな位置情報データの一部分を配信することを可能とするファイル形式の総称です。また、ラスターデータではすでにCloud Optimized形式が普及期であり、Cloud Optimized GeoTIFF（COG）がデファクトスタンダードになっています。」（引用：『現場のプロがわかりやすく教える 位置情報エンジニア養成講座』）本記事では、MapLibre GL JSを用いた、Cloud Optimized GeoTIFF（COG）の利用について説明します。

Cloud Optimized GeoTIFF（COG）の特徴

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

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

前提条件

  • OSGeo4Wがインストール済みであること。

GDALでCloud Optimized GeoTIFF（COG）を作成する

gdalinfo 08NF2340_AC.tif

出力結果

Driver: GTiff/GeoTIFF
Files: 08NF2340_AC.tif
       08NF2340_AC.tfw
Size is 2000, 1500
Origin = (52000.000000000000000,-97200.000000000000000)
Pixel Size = (0.200000000000000,-0.200000000000000)
Metadata:
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   52000.000,  -97200.000)
Lower Left  (   52000.000,  -97500.000)
Upper Right (   52400.000,  -97200.000)
Lower Right (   52400.000,  -97500.000)
Center      (   52200.000,  -97350.000)
Band 1 Block=2000x1 Type=Byte, ColorInterp=Red
Band 2 Block=2000x1 Type=Byte, ColorInterp=Green
Band 3 Block=2000x1 Type=Byte, ColorInterp=Blue
  • 下記のコマンドを実行してCloud Optimized GeoTIFF（COG）を作成します。
gdal_translate -a_srs EPSG:6676 -of COG 08NF2340_AC.tif 08NF2340_AC_cog.tif -co TILING_SCHEME=GoogleMapsCompatible -co COMPRESS=DEFLATE
  • 下記のコマンドを実行してCloud Optimized GeoTIFF（COG）のメタ情報を表示してみます。
gdalinfo 08NF2340_AC_cog.tif

出力結果

  • 上記で出力ファイルの座標参照系をEPSG:6676に設定しましたが、TILING_SCHEME=GoogleMapsCompatibleを追加することで、最終的な出力ファイルの座標参照系はWeb Mercator（EPSG:3857）になるようです。
Driver: GTiff/GeoTIFF
Files: 08NF2340_AC_cog.tif
Size is 1792, 1536
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 = (15481249.960766460746527,4180570.388001143932343)
Pixel Size = (0.298582141736655,-0.298582141738734)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  LAYOUT=COG
Tiling Scheme:
  NAME=GoogleMapsCompatible
  ZOOM_LEVEL=19
Corner Coordinates:
Upper Left  (15481249.961, 4180570.388) (139d 4'13.56"E, 35d 7'21.78"N)
Lower Left  (15481249.961, 4180111.766) (139d 4'13.56"E, 35d 7' 9.65"N)
Upper Right (15481785.020, 4180570.388) (139d 4'30.87"E, 35d 7'21.78"N)
Lower Right (15481785.020, 4180111.766) (139d 4'30.87"E, 35d 7' 9.65"N)
Center      (15481517.490, 4180341.077) (139d 4'22.22"E, 35d 7'15.71"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 896x768, 448x384, 224x192
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 896x768, 448x384, 224x192
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 896x768, 448x384, 224x192
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 896x768, 448x384, 224x192
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 896x768, 448x384, 224x192
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 896x768, 448x384, 224x192
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 896x768, 448x384, 224x192

Cloud Optimized GeoTIFF（COG）をMapLibre GL JSで表示する

  • 下記の記事のコードを参考にさせていただきました。
  • Cloud Optimized GeoTIFF（COG）のJavaScript実装はgeotiff.js一択のようです。
  • MapLibre GL JSでは画像ファイルとしてバイナリデータを渡してやる必要があるとのことなので、fast-pngを用いてPNG形式へエンコードしています。

コード

<html>

<head>
    <title>COG on MapLibre GL JS</title>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <script src="https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.js"></script>
    <link href="https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.css" rel="stylesheet">
    <script src="https://cdn.jsdelivr.net/npm/geotiff"></script>
    <style>
        body {
            margin: 0;
            padding: 0;
        }

        #map {
            position: absolute;
            top: 0;
            bottom: 0;
            width: 100%;
        }
    </style>
</head>

<body>
    <div id="map"></div>
    <script type="module">
        // 'GeoTIFF'パッケージから 'fromUrl' と 'Pool' 関数をインポート
        const { fromUrl, Pool } = GeoTIFF;

        // 'fast-png'パッケージから'encode'関数をインポート。これは画像データをPNG形式にエンコードするために使用。
        import { encode as fastPngEncode } from 'https://cdn.jsdelivr.net/npm/fast-png@6.1.0/+esm';

        // タイルのx/y/z座標をWebメルカトルの境界ボックスに変換する関数を定義
        const merc = (x, y, z) => {
            const GEO_R = 6378137;
            const orgX = -1 * ((2 * GEO_R * Math.PI) / 2);
            const orgY = (2 * GEO_R * Math.PI) / 2;
            const unit = (2 * GEO_R * Math.PI) / Math.pow(2, z);
            const minx = orgX + x * unit;
            const maxx = orgX + (x + 1) * unit;
            const miny = orgY - (y + 1) * unit;
            const maxy = orgY - y * unit;
            return [minx, miny, maxx, maxy];
        };

        // Cloud Optimized GeoTIFF(COG)ソースを生成する非同期関数を定義
        const generateCogSource = async (url) => {
            // GeoTIFFから指定されたURLの画像を読み込み
            const tiff = await fromUrl(url);
            const pool = new Pool();
            maplibregl.addProtocol('cog', (params, callback) => {
                const segments = params.url.split('/');
                console.log('segments:' + segments);
                // タイルインデックスの取得
                const [z, x, y] = segments.slice(segments.length - 3).map(Number);
                // webmercator-bboxを取得
                const bbox = merc(x, y, z);
                // 画像のサイズを指定
                const size = 256;
                // GeoTIFFから特定の境界ボックス(bbox)内のラスターデータを読み込み
                tiff.readRasters({
                    bbox,
                    samples: [0, 1, 2, 3], // 4つのバンド（赤、緑、青、アルファ）を取得する
                    width: size,  // 読み込む画像のサイズを指定
                    height: size, // 読み込む画像のサイズを指定
                    interleave: true, // バンドを交互に読み込む
                    pool, // 並列処理のためのワーカープールを指定
                }).then((data) => {
                    const img = new ImageData(
                        //@ts-ignore
                        new Uint8ClampedArray(data),
                        size,
                        size,
                    );
                    // imgデータをPNG形式にエンコード
                    const png = fastPngEncode(img);
                    callback(null, png, null, null);
                });
                return { cancel: () => { } };
            });

            // RasterSourceSpecificationに対応するソースオブジェクトを生成
            const source = {
                type: 'raster',
                tiles: [`cog://${url.split('://')[1]}/{z}/{x}/{y}`],
                tileSize: 256,
                minzoom: 2,
                maxzoom: 18,
                attribution: '<a href="https://www.geospatial.jp/ckan/dataset/aac-2021710-atami-lp">G空間情報センター 2021年7月10日静岡県熱海市土石流災害航空レーザ計測データ他(CC By 4.0)</a>'
            };

            return { source };
        };

        // マップを初期化
        const map = new maplibregl.Map({
            container: 'map',
            style: {
                version: 8,
                sources: {
                    osm: {
                        type: 'raster',
                        tiles: ['https://tile.openstreetmap.org/{z}/{x}/{y}.png'],
                        attribution:
                            '<a href="http://osm.org/copyright">OpenStreetMap contributors</a>',
                    },
                },
                layers: [
                    {
                        id: 'osm',
                        type: 'raster',
                        source: 'osm',
                    },
                ],
            },
            hash: true,
            zoom: 17,
            center: [139.073147, 35.12128],
            pitch: 0,
            bearing: 0,
            attributionControl: false,
        });


        // ズーム・回転
        map.addControl(new maplibregl.NavigationControl());

        // フルスクリーンモードのオンオフ
        map.addControl(new maplibregl.FullscreenControl());

        // 現在位置表示
        map.addControl(new maplibregl.GeolocateControl({
            positionOptions: {
                enableHighAccuracy: false
            },
            fitBoundsOptions: { maxZoom: 18 },
            trackUserLocation: true,
            showUserLocation: true
        }));

        // スケール表示
        map.addControl(new maplibregl.ScaleControl({
            maxWidth: 200,
            unit: 'metric'
        }));

        // Attributionを折りたたみ表示
        map.addControl(new maplibregl.AttributionControl({
            compact: true,
            customAttribution: '（<a href="https://twitter.com/syanseto" target="_blank">Twitter</a> | <a href="https://github.com/shi-works/cog-on-maplibre-gl-js" target="_blank">Github</a>） '
        }));

        // マップのロードが完了したら、COGソースを生成しレイヤーに追加
        map.on('load', async () => {
            const { source } = await generateCogSource(
                './08NF2340_AC_cog.tif'
            );
            map.addSource('cogSample', source);
            map.addLayer({
                id: 'cogLayer',
                type: 'raster',
                source: 'cogSample',
            });
        });

        map.showTileBoundaries = true;

    </script>
</body>

</html>

デモサイト

image.png

他の実装例

  • Sentinel HubのEO Browserから取得した、Sentinel-1の衛星画像（GeoTIFF）をCloud Optimized GeoTIFF（COG）に変換して、MapLibre GL JSで表示した事例になります。

image.png

参考文献

4
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?