rio-tilerを使ってCloud Optimized GeoTIFFをブラウザで表示する
はじめに
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
{
"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
ファイルの範囲外の座標を指定すると、rio_tiler.errors.TileOutsideBounds
で500エラーとなります。
Landsat-8を試す
Landsat-8という衛星のデータもCOGで配布されているのでrio-tilerを試してみます。
ここまで試してきたSkySatのサンプルは1ファイルの中に3バンド(RGB)すべてが含まれていましたが、Landsat-8のデータはバンドごとにファイルが分かれています。
とりあえず1ファイル、モノクロで表示するところを目標とします。
tilejson.json
{
"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
のあたりです。
しかし、呼び出してみても真っ黒な画像が表示されるだけのはずです。
というのも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])
画像の隅っこですが、インドのカタック付近、山や川の筋がぼんやり表示されるはずです[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の他の機能も調べてみると面白いでしょう。
-
COGにすることで全体のファイルサイズが大きくなり、ファイル全体の利用しか想定しない場合はむしろデメリットが勝ります。 ↩︎
-
URLが無効となってしまったときは別のサンプルファイルを探してみてください。 ↩︎
-
本記事ではpipenvを使いますがpipのままいきたい人は適宜読み替えてください。 ↩︎
-
https://www.usgs.gov/core-science-systems/nli/landsat/using-usgs-landsat-level-1-data-product ↩︎
-
AWSのLandsat-8はサムネイルのピラミッドを持たないため、zoomが小さいとレスポンスがかなり遅いです。 ↩︎
-
expressionを使うともっとキレイに書けます。 ↩︎
-
2.0.0rc1からサンプルから消えてしまいました。代わりに
rio_tiler.tasks.multi_arrays
を使えとのこと。 ↩︎
Discussion