rio-tilerを使ってLandsat-8のトゥルーカラー合成タイルを動的生成する
はじめに
Landsat-8という衛星のデータがさまざまなサイトで公開されています。
例えばAWS上ではシーンごとにindexが用意されており、データはCloud Optimized GeoTIFF （COG）の形式で入手可能です。
今回は衛星データをビジュアライズするときに定番のトゥルーカラー合成を、rio-tilerを使ってやってみます。
COGやrio-tilerがどんなものかはこちらの記事をどうぞ。
カラー合成
人工衛星は人間の目が赤、緑、青と感じる波長以外にも様々な波長で世界を観測しています。
そのため人が画像としてデータを確認するときはRGBにどの波長帯を割り当てるか決めないといけません。
代表的な組み合わせにはトゥルーカラー、フォルスカラー、ナチュラルカラーなど名前がついています。
トゥルーカラーは人間の目が赤、緑、青と感じる波長をそのままRGBに割り当てる合成方法です。
他の代表的な組み合わせはこちらなどを参照ください。
波長の組み合わせによって見えるものが変わってきます[1]。
Landsat-8のセンサは11バンドの観測波長を持っており、1バンドにつき1つのファイルで配布されています。
赤、緑、青はバンド4、３、２にあたり、どのファイルかはファイル名の末尾（BX.TIFのX部分）で判断します。
DN値タイル
ツールをインストールし、app.pyを用意して起動します[2]。
pipenv install fastapi uvicorn rio-tiler==2.0.0b19
import math
import os
import re
import requests
from enum import Enum
from typing import Any, Dict, List, Union, Optional
import numpy as np
import uvicorn
from fastapi import FastAPI, Path, Query, HTTPException
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 rio_tiler.tasks import multi_arrays
from rio_tiler.errors import TileOutsideBounds
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("/landsat8_dn_tile/{z}/{x}/{y}", **tile_routes_params)
def landsat8_dn_tile(
z: int,
x: int,
y: int,
base_url: str = Query(..., description="Landsat-8 Base URL."),
):
bands = [{
'url':'{}_B{}.TIF'.format(base_url, band)
} for band in [4,3,2]]
def _reader(band, *args: Any, **kwargs: Any):
with COGReader(band['url']) as cog:
try:
tile, mask = cog.tile(x, y, z, tilesize=256)
mask[tile[0] == 0] = 0
tile[tile < 0] = 0
return tile, mask
except TileOutsideBounds as e:
raise HTTPException(status_code=404, detail=str(e))
output = multi_arrays(bands, _reader)
tile = (255 * output[0] / 4095).astype(np.uint8) # 12bit to 8bit
mask = output[1]
driver = "PNG"
options = img_profiles.get(driver.lower(), {})
img = render(tile, mask, img_format=driver, **options)
return TileResponse(img, media_type="image/png")
pipenv shell
uvicorn app:app --reload
このコードは対象のシーンのファイル一式が置かれているパスをbase_urlで指定するように作っています。
AWSなら例えばこんなパスです。
https://landsat-pds.s3.amazonaws.com/c1/L8/139/045/LC08_L1TP_139045_20170304_20170316_01_T1/LC08_L1TP_139045_20170304_20170316_01_T1
そしてそのシーンに含まれている範囲のタイル座標を指定します。ここではz=7, x=94, y=56としてみましょう（表示まで数十秒かかるので辛抱強く待ちます）。
結果がこちら。ほとんど真っ黒で、目を細めるとモヤのようなものが見えるかもしれません。
というのも、Landsat-8のTIFFの各値は12bitのデジタルナンバーですが、ヒストグラムで確認すると多くの場合で左側に偏っています。左端右端には値がほとんどないケースも多いです。
そのため単純に0から4095でレンジを取ってしまうとこんな画像ができあがってしまいます。
人の目で見てわかる形にするには、画像毎にヒストグラム平滑化やダイナミックレンジの調整が必要です。
シグモイド関数フィルタ＋反射率タイル
濃度変換のやり方はいろいろありますが、今回はシグモイド関数フィルタ（S字トーンカーブによる変換）[3]を使うことにします[4]。
またこのときDN値をそのまま使うのでなく大気上端反射率に変換してからフィルタにかけることにします[5]。
大気上端反射率への変換は公式にある式をそのまま利用します。
@app.get("/landsat8_tile/{z}/{x}/{y}", **tile_routes_params)
def landsat8_tile(
z: int,
x: int,
y: int,
base_url: str = Query(..., description="Landsat-8 Base URL."),
a: Union[int, float] = Query(40, description="Coefficient 'a' Of Sigmoid Filter."),
b: Union[int, float] = Query(0.15, description="Coefficient 'b' Of Sigmoid Filter."),
):
mlt_path = '{}_MTL.txt'.format(base_url)
if base_url.startswith('http') or base_url.startswith('//:'):
r = requests.get(mlt_path)
mlt = r.text
else:
with open(mlt_path) as f:
mlt = f.read()
REFLECTANCE_MULT_BANDS = re.findall('REFLECTANCE_MULT_BAND_\d+\s+=\s+(.*)\n', mlt)
REFLECTANCE_ADD_BANDS = re.findall('REFLECTANCE_ADD_BAND_\d+\s+=\s+(.*)\n', mlt)
SUN_ELEVATION = re.findall('SUN_ELEVATION\s+=\s+(.*)\n', mlt)
denominator = math.sin(math.radians(float(SUN_ELEVATION[0])))
bands = [{
'url':'{}_B{}.TIF'.format(base_url, band),
'mult': float(REFLECTANCE_MULT_BANDS[band - 1]),
'add': float(REFLECTANCE_ADD_BANDS[band - 1]),
'denominator': denominator
} for band in [4,3,2]]
def _reader(band, *args: Any, **kwargs: Any):
with COGReader(band['url']) as cog:
try:
tile, mask = cog.tile(x, y, z, tilesize=256)
mask[tile[0] == 0] = 0
tile = (tile * band['mult'] + band['add'])/band['denominator']
tile[tile < 0] = 0
return tile, mask
except TileOutsideBounds as e:
raise HTTPException(status_code=404, detail=str(e))
output = multi_arrays(bands, _reader)
f0 = 1/(1 + np.exp(a*b))
f1 = 1/(1 + np.exp(a*(b - 1)))
tile = ((255/(1 + np.exp(a*(b - output[0]))) - f0) / (f1 - f0)).astype(np.uint8)
mask = output[1]
driver = "PNG"
options = img_profiles.get(driver.lower(), {})
img = render(tile, mask, img_format=driver, **options)
return TileResponse(img, media_type="image/png")
変換に必要な係数はmlt.txtに記述されているので、TIFFを取得する前にtxtファイルを取得する処理を追加しています。
tile = (tile * band['mult'] + band['add'])/denominatorの部分は、本当ならexpressionにして
tile, mask = cog.tile(x, y, z, expression, tilesize=256)で処理したいところですが、mask[tile[0] == 0] = 0でnodata部分をmaskに追加する必要があるためこのようにします。
これをやらないと、0の部分がmaskされず画像の端が黒くなってしまいます。
本来ならgeotiffにNodata値を設定しておくべきです。
また以降は必要なファイル（{シーンID}_mlt.txt, {シーンID}_B2.TIF, {シーンID}_B3.TIF, {シーンID}_B4.TIF）をすべてローカルに保存しているものとします(シーンIDのディレクトリを作って保存します。このときにTIFFの拡張子は.TIFとしてください)。
AWSに置いてあるCOGはOverviewsが用意されてないせいか、タイルを作るのに時間がかかってしまうためです。
もちろんAWSのパスのままでも基本は動きます。
シグモイド関数フィルタは次の式で実装しており、クエリでa、bを指定することができます。デフォルト値はa=40, b=0.15です。
パラメータを変えた例
これをOpenlayersを使って表示してやると、以下のようにキレイな絵を見ることができます。
openlayers.html
<!doctype html>
<html lang="ja">
<head>
<link rel="stylesheet"
href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.4.3/css/ol.css" type="text/css">
<style>
.map {
height: 480px;
width: 100%;
}
</style>
<script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.4.3/build/ol.js"></script>
<title>OpenLayers example</title>
</head>
<body>
<div id="map" class="map"></div>
<script type="text/javascript">
var sample = new ol.layer.Tile({
source: new ol.source.XYZ({
maxZoom: 12,
url: "http://127.0.0.1:8000/landsat8_tile/{z}/{x}/{y}?a=50&b=0.15&base_url=LC08_L1TP_139045_20170304_20170316_01_T1/LC08_L1TP_139045_20170304_20170316_01_T1",
crossOrigin: 'anonymous'
})
});
var map = new ol.Map({
target: 'map',
layers: [
new ol.layer.Tile({
source: new ol.source.OSM()
}),
sample,
new ol.layer.Graticule({
strokeStyle: new ol.style.Stroke({
color: 'rgba(255,120,0,0.9)',
width: 2,
lineDash: [0.5, 4],
}),
showLabels: true,
wrapX: false,
})
],
view: new ol.View({
center: ol.proj.fromLonLat([87, 21.5]),
zoom: 8
})
});
</script>
</body>
</html>
mosaic
この記事で利用しているmosaic機能(mosaic_reader)と組み合わせることもできます。
from rio_tiler.mosaic import mosaic_reader
assets = [
"LC08_L1TP_139045_20170304_20170316_01_T1/LC08_L1TP_139045_20170304_20170316_01_T1",
"LC08_L1TP_139046_20170304_20170316_01_T1/LC08_L1TP_139046_20170304_20170316_01_T1",
"LC08_L1TP_139044_20170304_20170316_01_T1/LC08_L1TP_139044_20170304_20170316_01_T1"
]
@app.get("/landsat8_mosaic/{z}/{x}/{y}", **tile_routes_params)
def landsat8_mosaic(
z: int,
x: int,
y: int,
a: Union[int, float] = Query(40, description="Coefficient 'a' Of Sigmoid Filter."),
b: Union[int, float] = Query(0.15, description="Coefficient 'b' Of Sigmoid Filter."),
):
def _reader(band, *args: Any, **kwargs: Any):
with COGReader(band['url']) as cog:
tile, mask = cog.tile(x, y, z, tilesize=256)
mask[tile[0] == 0] = 0
tile = (tile * band['mult'] + band['add'])/band['denominator']
tile[tile < 0] = 0
return tile, mask
def _assets_reader(base_url):
mlt_path = '{}_MTL.txt'.format(base_url)
if base_url.startswith('http') or base_url.startswith('//:'):
r = requests.get(mlt_path)
mlt = r.text
else:
with open(mlt_path) as f:
mlt = f.read()
REFLECTANCE_MULT_BANDS = re.findall('REFLECTANCE_MULT_BAND_\d+\s+=\s+(.*)\n', mlt)
REFLECTANCE_ADD_BANDS = re.findall('REFLECTANCE_ADD_BAND_\d+\s+=\s+(.*)\n', mlt)
SUN_ELEVATION = re.findall('SUN_ELEVATION\s+=\s+(.*)\n', mlt)
denominator = math.sin(math.radians(float(SUN_ELEVATION[0])))
bands = [{
'url':'{}_B{}.TIF'.format(base_url, band),
'mult': float(REFLECTANCE_MULT_BANDS[band - 1]),
'add': float(REFLECTANCE_ADD_BANDS[band - 1]),
'denominator': denominator
} for band in [4,3,2]]
return multi_arrays(bands, _reader)
output, _ = mosaic_reader(assets, _assets_reader)
if output[0] is None:
raise HTTPException(status_code=404, detail="Tile Outside Bounds.")
f0 = 1/(1 + np.exp(a*b))
f1 = 1/(1 + np.exp(a*(b - 1)))
tile = ((255/(1 + np.exp(a*(b - output[0]))) - f0) / (f1 - f0)).astype(np.uint8)
mask = output[1]
driver = "PNG"
options = img_profiles.get(driver.lower(), {})
img = render(tile, mask, img_format=driver, **options)
return TileResponse(img, media_type="image/png")
openlayers.htmlのスクリプト部分も書き換えます。
var mosaic = new ol.layer.Tile({
source: new ol.source.XYZ({
maxZoom: 12,
url: "http://127.0.0.1:8000/landsat8_mosaic/{z}/{x}/{y}?a=50&b=0.15",
crossOrigin: 'anonymous'
})
});
var map = new ol.Map({
target: 'map',
layers: [
new ol.layer.Tile({
source: new ol.source.OSM()
}),
mosaic,
new ol.layer.Graticule({
strokeStyle: new ol.style.Stroke({
color: 'rgba(255,120,0,0.9)',
width: 2,
lineDash: [0.5, 4],
}),
showLabels: true,
wrapX: false,
})
],
view: new ol.View({
center: ol.proj.fromLonLat([87, 21.5]),
zoom: 8
})
});
複数のシーンを組みあわせたレイヤーを表示することができました。
まとめ
rio-tilerのmulti_arraysを使うことでバンドが別ファイルに別れている衛星データのカラー合成をすることができました。
multi_arraysはmosaic_readerと組み合わせることも可能です。
しかもこのmulti_arraysとmosaic_readerは各ファイルを並列処理で呼び出してしてくれる優れもので、パフォーマンスを考えてよく作られています。v2.0が正式にリリースされるのが今から楽しみです[6]。
Discussion