Zenn
🙄

rio-tilerを使ってCloud Optimized GeoTIFFをブラウザで表示する

2020/11/03に公開
Python
#
gdal
#
cog
#
riotiler
tech

はじめに

GeoTIFFとはTIFFに地理座標を埋め込んだファイルフォーマットです。
さらにGeoTIFFをHTTP Rangeリクエストに対応させたものをCloud Optimized GeoTIFF （COG）と呼び、衛星データのフォーマットとして近頃利用されています(COGについて日本語でもう少し詳しく知りたい人はこちらをどうぞ)。

ファイル全体を利用したい人にとってCOGの恩恵はありません[1]が、ファイルのうち欲しいデータは一部分だけというようなケースではレスポンスが速くとてもメリットが大きいです。

gdallocationinfoを利用して1点の値を取得する例[2] 投影されているgeotiffであればこのコマンドは実行可能

gdallocationinfo /vsicurl/https://s3-us-west-2.amazonaws.com/planet-disaster-data/hurricane-harvey/SkySat_Freeport_s03_20170831T162740Z3.tif -wgs84 -95.321041 28.951376
> Report:
>   Location: (17817P,16000L)
>   Band 1:
>     Value: 147
>   Band 2:
>     Value: 144
>   Band 3:
>     Value: 142

そんなCOGをブラウザで表示するには、geotiff.jsを利用すると便利です。

しかしCOG単体で使うならこれでいいのですが、地理院タイルなど他のリソースと合わせて使おうとすると、そのままではうまくいかず色々と面倒です。

ウェブメルカトル投影された地図タイルを使うならOpenLayersを利用したいところですが、OpenLayersは現在(2020年時点)はCOGに対応していません。

従来であればGeoTIFFをブラウザで他のリソースと一緒に利用するのは、gdal2tilesでタイルに変換して大量のpngファイルを静的配信するのが常道でしたが、変換に時間はかかるし容量も食うので気軽にやりたいことではありませんでした。

ところが最近ではrio-tilerのDynamic Tiler機能を使うと、COGを事前に変換処理すること無く地図タイルとして呼び出すことができるようになりました。

本記事では実際にこれを試してみることにします。

サンプルを動かす

Dynamic Tilerのサンプルを動かしてみます。

ツールをインストールします。[3]

pipenv install fastapi uvicorn rio-tiler==2.0.0b19

app.pyを用意し起動します。サンプルのままではエラーが出たので108行目のrequest.url_forの引数を書き換えています。

"""rio-tiler tile server."""

import os
from enum import Enum
from typing import Any, Dict, List, Optional
from urllib.parse import urlencode

import uvicorn
from fastapi import FastAPI, Path, Query
from rasterio.crs import CRS
from starlette.background import BackgroundTask
from starlette.middleware.cors import CORSMiddleware
from starlette.middleware.gzip import GZipMiddleware
from starlette.requests import Request
from starlette.responses import Response

from rio_tiler.profiles import img_profiles
from rio_tiler.utils import render
from rio_tiler.io import COGReader

# From developmentseed/titiler
drivers = dict(jpg="JPEG", png="PNG")
mimetype = dict(png="image/png", jpg="image/jpg",)

class ImageType(str, Enum):
    """Image Type Enums."""

    png = "png"
    jpg = "jpg"



class TileResponse(Response):
    """Tiler's response."""

    def __init__(
        self,
        content: bytes,
        media_type: str,
        status_code: int = 200,
        headers: dict = {},
        background: BackgroundTask = None,
        ttl: int = 3600,
    ) -> None:
        """Init tiler response."""
        headers.update({"Content-Type": media_type})
        if ttl:
            headers.update({"Cache-Control": "max-age=3600"})
        self.body = self.render(content)
        self.status_code = 200
        self.media_type = media_type
        self.background = background
        self.init_headers(headers)


app = FastAPI(
    title="rio-tiler",
    description="A lightweight Cloud Optimized GeoTIFF tile server",
)
app.add_middleware(
    CORSMiddleware,
    allow_origins=["*"],
    allow_credentials=True,
    allow_methods=["GET"],
    allow_headers=["*"],
)
app.add_middleware(GZipMiddleware, minimum_size=0)

responses = {
    200: {
        "content": {"image/png": {}, "image/jpg": {}},
        "description": "Return an image.",
    }
}
tile_routes_params: Dict[str, Any] = dict(
    responses=responses, tags=["tiles"], response_class=TileResponse
)


@app.get("/{z}/{x}/{y}", **tile_routes_params)
def tile(
    z: int,
    x: int,
    y: int,
    url: str = Query(..., description="Cloud Optimized GeoTIFF URL."),
):
    """Handle tiles requests."""
    with COGReader(url) as cog:
        tile, mask = cog.tile(x, y, z, tilesize=256)

    format = ImageType.jpg if mask.all() else ImageType.png

    driver = drivers[format.value]
    options = img_profiles.get(driver.lower(), {})
    img = render(tile, mask, img_format=driver, **options)

    return TileResponse(img, media_type=mimetype[format.value])


@app.get("/tilejson.json", responses={200: {"description": "Return a tilejson"}})
def tilejson(
    request: Request,
    url: str = Query(..., description="Cloud Optimized GeoTIFF URL."),
    minzoom: Optional[int] = Query(None, description="Overwrite default minzoom."),
    maxzoom: Optional[int] = Query(None, description="Overwrite default maxzoom."),
):
    """Return TileJSON document for a COG."""
    tile_url = request.url_for("tile", z="{z}", x="{x}", y="{y}").replace("\\", "")

    kwargs = dict(request.query_params)
    kwargs.pop("tile_format", None)
    kwargs.pop("tile_scale", None)
    kwargs.pop("minzoom", None)
    kwargs.pop("maxzoom", None)

    qs = urlencode(list(kwargs.items()))
    tile_url = f"{tile_url}?{qs}"

    with COGReader(url) as cog:
        center = list(cog.center)
        if minzoom:
            center[-1] = minzoom
        tjson = {
            "bounds": cog.bounds,
            "center": tuple(center),
            "minzoom": minzoom or cog.minzoom,
            "maxzoom": maxzoom or cog.maxzoom,
            "name": os.path.basename(url),
            "tiles": [tile_url],
        }

    return tjson

pipenv shell
uvicorn app:app --reload

あとはcurlでもブラウザでも良いので、
http://127.0.0.1:8000/tilejson.json?url={cogファイルのURL}を呼び出します。

tilejson
http://127.0.0.1:8000/tilejson.json?url=https://s3-us-west-2.amazonaws.com/planet-disaster-data/hurricane-harvey/SkySat_Freeport_s03_20170831T162740Z3.tif

{
   "bounds":[
      -95.46993601221669,
      28.869054414503687,
      -95.23861524309007,
      29.068191707419608
   ],
   "center":[
      -95.35427562765338,
      28.968623060961647,
      10
   ],
   "minzoom":10,
   "maxzoom":17,
   "name":"SkySat_Freeport_s03_20170831T162740Z3.tif",
   "tiles":[
      "http://127.0.0.1:8000/{z}/{x}/{y}?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif"
   ]
}

tilejson.jsonはCOGのメタ情報から座標やタイルにしたときのzoomの範囲、地図タイルのエンドポイント(tiles)を返してくれます。
centerの値から https://gsj-seamless.jp/labs/tools/tileCoord.html でタイル座標を調べると、z=10, x=240, y=425近辺なのでzを変化させて数タイル呼び出してみます(tilesの{z}/{x}/{y}を置き換える)。

tile
http://127.0.0.1:8000/10/240/425?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif

http://127.0.0.1:8000/12/962/1702?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif

http://127.0.0.1:8000/16/15401/27241?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif

ファイルの範囲外の座標を指定すると、rio_tiler.errors.TileOutsideBoundsで500エラーとなります。

Landsat-8を試す

Landsat-8という衛星のデータもCOGで配布されているのでrio-tilerを試してみます。
https://landsat-pds.s3.amazonaws.com/c1/L8/139/045/LC08_L1TP_139045_20170304_20170316_01_T1/index.html

ここまで試してきたSkySatのサンプルは1ファイルの中に3バンド(RGB)すべてが含まれていましたが、Landsat-8のデータはバンドごとにファイルが分かれています。

とりあえず1ファイル、モノクロで表示するところを目標とします。

tilejson.json
http://127.0.0.1:8000/tilejson.json?url=https://landsat-pds.s3.amazonaws.com/c1/L8/139/045/LC08_L1TP_139045_20170304_20170316_01_T1/LC08_L1TP_139045_20170304_20170316_01_T1_B1.TIF

{
   "bounds":[
      85.85282002580081,
      20.609905339496255,
      88.07610888865462,
      22.719911220251923
   ],
   "center":[
      86.9644644572277,
      21.66490827987409,
      7
   ],
   "minzoom":7,
   "maxzoom":12,
   "name":"LC08_L1TP_139045_20170304_20170316_01_T1_B1.TIF",
   "tiles":[
      "http://127.0.0.1:8000/{z}/{x}/{y}?url=https%3A%2F%2Flandsat-pds.s3.amazonaws.com%2Fc1%2FL8%2F139%2F045%2FLC08_L1TP_139045_20170304_20170316_01_T1%2FLC08_L1TP_139045_20170304_20170316_01_T1_B1.TIF"
   ]
}

このシーンの場合、座標はだいたいz=7, x=94, y=56のあたりです。

http://localhost:8000/7/94/56?url=https%3A%2F%2Flandsat-pds.s3.amazonaws.com%2Fc1%2FL8%2F139%2F045%2FLC08_L1TP_139045_20170304_20170316_01_T1%2FLC08_L1TP_139045_20170304_20170316_01_T1_B3.TIF

しかし、呼び出してみても真っ黒な画像が表示されるだけのはずです。
というのもLandsat-8のデータはDN（デジタルナンバー）といってそのままでは画像として表示するには不適当な値となっており、人が見てそれらしい画像にするためには大気上端反射率(TOA Reflectance)などになおしてあげる必要があります[4][5]。

補正に必要な値は{シーンID}_MTL.txtで知ることができます。

本当はシーンIDからMTL.txtを取得し、必要な値を読み取り大気上端反射率に変換したタイルを返してくれるエンドポイントを作れるといいのですが、今回はあくまでデモということで、シーンIDはLC08_L1TP_139045_20170304_20170316_01_T1固定でタイルを返すエンドポイントを作ります。

@app.get("/test_landsat8/{z}/{x}/{y}", **tile_routes_params)
def test_landsat8(
    z: int,
    x: int,
    y: int,
):
    url = "https://landsat-pds.s3.amazonaws.com/c1/L8/139/045/LC08_L1TP_139045_20170304_20170316_01_T1/LC08_L1TP_139045_20170304_20170316_01_T1_B1.TIF"

    # from LC08_L1TP_139045_20170304_20170316_01_T1_MTL.txt
    REFLECTANCE_MULT_BAND_1 = 2.0000E-05
    REFLECTANCE_ADD_BAND_1 = -0.100000
    SUN_ELEVATION = 51.69540297

    with COGReader(url) as cog:
        tile, mask = cog.tile(x, y, z, tilesize=256)

    tile = (tile * REFLECTANCE_MULT_BAND_1 + REFLECTANCE_ADD_BAND_1) \
            / math.sin(math.radians(SUN_ELEVATION))
    tile[tile<0] = 0
    tile=(tile*255).astype(np.uint8)

    format = ImageType.jpg if mask.all() else ImageType.png

    driver = drivers[format.value]
    options = img_profiles.get(driver.lower(), {})
    img = render(tile, mask, img_format=driver, **options)

    return TileResponse(img, media_type=mimetype[format.value])

http://localhost:8000/test_landsat8/7/94/56

画像の隅っこですが、インドのカタック付近、山や川の筋がぼんやり表示されるはずです[6]。

tile, mask = cog.tile(x, y, z, tilesize=256)のtileはnumpyの配列(256 * 256)なため演算も簡単です[7]。

あとはmulti_assets[8]を参考に各バンドを並列で呼び出してやれば、任意のバンドを組み合わせてカラー合成ができるはずです。

まとめ

rio-tilerを使うと簡単に地図タイルをCOGから呼び出すことができるのをサンプルを使って紹介しました。

Landsat-8の例のように、標高データのCOGを使うと、クエリによって標高Xm以下は塗り潰すなど条件を変えたタイルを動的に作ることもできそうです。

rio-tilerの他の機能も調べてみると面白いでしょう。

脚注

  1. COGにすることで全体のファイルサイズが大きくなり、ファイル全体の利用しか想定しない場合はむしろデメリットが勝ります。 ↩︎

  2. URLが無効となってしまったときは別のサンプルファイルを探してみてください。 ↩︎

  3. 本記事ではpipenvを使いますがpipのままいきたい人は適宜読み替えてください。 ↩︎

  4. http://rs.aoyaman.com/seminar/about5.html ↩︎

  5. https://www.usgs.gov/core-science-systems/nli/landsat/using-usgs-landsat-level-1-data-product ↩︎

  6. AWSのLandsat-8はサムネイルのピラミッドを持たないため、zoomが小さいとレスポンスがかなり遅いです。 ↩︎

  7. expressionを使うともっとキレイに書けます。 ↩︎

  8. 2.0.0rc1からサンプルから消えてしまいました。代わりにrio_tiler.tasks.multi_arraysを使えとのこと。 ↩︎

Discussion