MapLibre GL JSとCOGによるサーバーレスなラスタータイル配信・表示 #TypeScript

はじめに - Cloud-Optimizedという用語について

昨今のFOSS4G界隈のキーワードのひとつはCloud-Optimizedです。これはいわゆる「サーバーレス」を志向するファイルフォーマットを表す形容詞です。具体的には、HTTP-Range-Requestを活用し、巨大なデータから必要最小限の部分のみを取り出せるようにした形式と言えます。言うなればリモートファイルへのランダムアクセスで、この仕組みを使えば、ウェブサーバーさえあればGISサーバーが必要なくなる、つまり「サーバーレス」が実現します。

Range-Requestでは、ブラウザからサーバーへ「AというファイルのBからCの範囲のバイト列を返してくれ」というリクエストを出す事になります。しかし、リモートサーバーにあるファイルのバイト列のうちブラウザで必要な範囲がどこなのかは、何もアテがなければ特定することが出来ません。このアテを用意するのがCloud-Optimizedなファイル形式と言えます。

Cloud-Optimizedなファイル形式では、バイト列の先頭に、内部に保存されているデータの情報を格納するための一定のバイト長が確保されています。ブラウザではまず、その一定の長さのバイト列をRange-Requestで読み出し、データがどのように格納されているかを把握します。これをアテにして、抽出したいデータのバイト列の範囲をブラウザ側で計算し改めてRange-Requestを送ることで、必要最小限のデータのみを取得することが出来ます。下記で紹介する各種フォーマットのExampleでブラウザの通信を確認すれば、どのようにしてRange-Requestが発出されているか確認することができるでしょう。

ベクター・ラスターでそれぞれで、Cloud-Optimizedな形式が実装されています。ベクターはFlatGeobufが最も成熟しており、GeoParquetが追いかける形です。ラスターではCloud-Optimized GeoTIFF(COG)がデファクトスタンダードです。変わり種ではPMTilesというタイル志向の規格も実装されています。

私の観測範囲ではCOGはすでに普及期のようです（従来のGeoTIFFの機能拡張なので普及しやすいのかも）。

本記事では、COGを用いてMapLibre GL JSのラスタータイルレイヤーをサーバーレスで実現出来るか実験した結果を紹介します。なお実現はできて、条件付きで実用性もあるだろうという結論に至りました。

データについて

本記事で用いるデータは下記で取得出来るオルソ画像を加工したものとします。
航空レーザー測量データ（平成28年　地方創生推進交付金）【当麻町】

なお加工したデータのファイル名はinput.tifとします。

先行事例

このあたりの記事を読めば、COGを用いたラスタータイルサーバーを構築することができるでしょう。そのうえで本記事は、COGの仕組みがあるならば、そもそもサーバーを介さずにラスタータイル実装を実現出来るはずという仮説を検証することを目的とします。

Cloud-Optimized GeoTIFF(COG)を作成する

COGはGDALで簡単に作る事ができます。ラスタータイルとして使いたいので、ここではTILING_SCHEMEオプションを設定し、ウェブメルカトルに適したCOGを作成します。

gdal_translate input.tif -of COG cog.tif -co TILING_SCHEME=GoogleMapsCompatible -co COMPRESS=DEFLATE

EPSG:3857に再投影されているのでデータ自体は変わっているものの、QGISなどで表示する限りは変換前とデータにほぼ差がない事がわかります。普通のGeoTIFFとしてそのまま使えるのもCOGの利点と言えるでしょう。

ここで、投影法のほかブロックサイズも調整の余地があるようですが、実際にどのようにパフォーマンスに影響するかは深掘りしていません。

下記のようなGeoTIFFファイルとなりました。

gdalinfo input.tif 
Driver: GTiff/GeoTIFF
Files: input.tif
Size is 37632, 21248
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        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 = (15862823.605966057628393,5450800.924044203013182)
Pixel Size = (0.298582141739031,-0.298582141738910)
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  (15862823.606, 5450800.924) (142d29'53.41"E, 43d54'19.13"N)
Lower Left  (15862823.606, 5444456.651) (142d29'53.41"E, 43d51'51.26"N)
Upper Right (15874059.849, 5450800.924) (142d35'56.78"E, 43d54'19.13"N)
Lower Right (15874059.849, 5444456.651) (142d35'56.78"E, 43d51'51.26"N)
Center      (15868441.728, 5447628.787) (142d32'55.09"E, 43d53' 5.21"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 18816x10624, 9408x5312, 4704x2656, 2352x1328, 1176x664, 588x332, 294x166, 147x83

一般的なウェブサーバーならどこでもいいのですが、今回はAWS S3にアップロードしておくことにします。
以降、URLはhttps://<BUCKET_NAME>.s3.ap-northeast-1.amazonaws.com/cog.tifとします。

MapLibre GL JSとgeotiff.jsでCOGタイルを実装する

ブラウザでの地図ライブラリはMapLibre GL JSを選択しました。個人的に最もよく使うライブラリであることと、最近実装された仕組みを試したかったからです。ラスタータイルはふつうの地図ライブラリでは必ず実装があるものなので、ひとつのライブラリで検証が成功すれば横展開は難しくないはずです。

また、COGのJavaScript実装はgeotiff.js一択となります。以下、全コードと実装のポイントを示します。

実装のポイント

MapLibre GL JSのaddProtocol()を用いて、COGへのタイルリクエストの振る舞いを定義する（ちなみにMapbox GL JSにはこの実装はありません）
geotiff.jsのGeoTIFFBase.readRasters()を用いて、リモートサーバーのCOGに対してRange-Requestを送る
（これは今回の本筋ではありませんが補足）readRasters()で取得した値はUint8Arrayとなります（0-255の範囲でのRGBA値）。一方でMapLibre GL JSでは画像ファイルとしてのバイナリデータを渡してやる必要があるので、fast-pngを用いてPNG形式へエンコードしています。canvasでやるより100倍程度高速でした。

コード

検証なのであまり綺麗な実装ではありません。

HTML

<div id="map" style="height:90vh;"></div>

TypeScript

import maplibregl, { Map, RasterSourceSpecification } from 'maplibre-gl';
import 'maplibre-gl/dist/maplibre-gl.css';
import { Pool, fromUrl } from 'geotiff';
import { encode } from 'fast-png';

/**
 * transform x/y/z to webmercator-bbox
 * @param x
 * @param y
 * @param z
 * @returns {number[]} [minx, miny, maxx, maxy]
 */
const merc = (x: number, y: number, z: number): number[] => {
    // 参考: https://qiita.com/MALORGIS/items/1a9114dd090e5b891bf7
    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];
};

const generateCogSource = async (
    url: string,
): Promise<{ source: RasterSourceSpecification }> => {
    const tiff = await fromUrl(url);
    const pool = new Pool();
    maplibregl.addProtocol('cog', (params, callback) => {
        const segments = params.url.split('/');
        const [z, x, y] = segments.slice(segments.length - 3).map(Number);
        const bbox = merc(x, y, z);
        const size = 256;
        tiff.readRasters({
            bbox,
            samples: [0, 1, 2, 3], // 取得するバンドを指定
            width: size,
            height: size,
            interleave: true,
            pool,
        }).then((data) => {
            const img = new ImageData(
                //@ts-ignore
                new Uint8ClampedArray(data),
                size,
                size,
            );
            const png = encode(img);
            callback(null, png, null, null);
        });
        return { cancel: () => {} };
    });
    const source: RasterSourceSpecification = {
        type: 'raster',
        tiles: [`cog://${url.split('://')[1]}/{z}/{x}/{y}`],
        tileSize: 256,
        minzoom: 8,
        maxzoom: 16, // 今回用いたCOGは19までの解像度を持っているようだが制限しておく
    };
    return { source };
};

const map = new 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',
            },
        ],
    },
});

map.on('load', async () => {
    const { source } = await generateCogSource(
        'https://<BUCKET_NAME>.s3.ap-northeast-1.amazonaws.com/cog.tif',
    );
    map.addSource('cogSample', source);
    map.addLayer({
        id: 'cogLayer',
        type: 'raster',
        source: 'cogSample',
    });
});

結果と結論

上記の実装で、サーバーを用いずCOGをラスタータイルとして表示することができました。

なお下記のような注意点があります。以下が許容・迂回できるなら実用できると思います。

注意点

ちょっと重たい

画像の読み込み時にややカクツキがあります。読み込み速度自体もタイルサーバーを用いるときよりは若干遅いです（これは当然か）。

一度に大量のリクエストが発生すると、geotiff.jsがエラーを出すことがある

現状未解決のIssueのようです。なるべくリクエストが少なくなるようにした方がよさそうです（MapLibre GL JSで言えば、tileSizeを大きくする・maxzoomをやや低くしておくなど）。

また、エラーが発生した場合はそのタイルの読み込みがストップしてしまいます。