rio-tilerを使って標高地形図を動的に生成する
はじめに
rio-tilerを使うとCloud Optimized GeoTIFF(COG)から動的に地図タイルを生成できることは、こちらの記事で紹介しました。
本記事では実際にrio-tilerのmosaic_readerを使って標高地形図を生成し、Openlayersで表示してみます。
標高データを取得する
まずは標高データを取得します。
今回はJAXAが配布するALOS全球数値地表モデル(AW3D30)を利用します[1]。
ただし、AW3D30はDSM(Digital Surface Model)で地表面の高さ(樹木やビルの高さを含む)であり、厳密には標高とは異なることを知っておいてください。
データは 1deg * 1deg毎に1つのGeoTIFFとして、全球(陸域のすべて)配布されています。
ユーザー登録後、画面をポチポチしてzipをダウンロードします。
今回は全球はめんどうなので東経138 ~ 141deg, 北緯35 ~ 37degの範囲だけダウンロードしました。
ダウンロードしたGeoTIFFはGDALを使ってCOG化していきます[2]。ファイル分だけ繰り返してください。
gdal_translate N035E138/ALPSMLC30_N035E138_DSM.tif temp.tif -co TILED=YES
gdaladdo -r nearest temp.tif 2 4 8 16
gdal_translate temp.tif ALPSMLC30_N035E138_DSM2.tif -co TILED=YES -co COPY_SRC_OVERVIEWS=YES
できたCOGはapp.pyと同じディレクトリに、COGというディレクトリを作成しそこへ移動します[3]。
Dynamic Tiler
app.pyに次のコードを追加してサンプルを動かしてみます。
import numpy as np
from rio_tiler.colormap import cmap
from typing import Union
@app.get("/dsm_test/{z}/{x}/{y}", **tile_routes_params)
def dsm_test(
z: int,
x: int,
y: int,
url: str = Query(..., description="Cloud Optimized GeoTIFF URL."),
min: Union[int, float] = Query(0, description="Minimum Height."),
max: Union[int, float] = Query(8848, description="Maximum Height."),
):
if -9999 >= min:
min = 0
if min >= max:
max = 8848
with COGReader(url) as cog:
tile, mask = cog.tile(x, y, z, tilesize=256)
mask[tile[0] == -9999] = 0
tile[tile < min] = min
tile[tile > max] = max
tile = (tile/(max - min) * 255).astype(np.uint8)
format = ImageType.jpg if mask.all() else ImageType.png
driver = drivers[format.value]
options = img_profiles.get(driver.lower(), {})
colormap = cmap.get('jet')
img = render(tile, mask, img_format=driver, **options, colormap=colormap)
return TileResponse(img, media_type=mimetype[format.value])
プロダクト説明書によると、「無効画素は”-9999”を格納」とのことなのでmaskします。
min/maxをクエリで指定できるようにしましたが、未指定の場合はmin=0, max=8848(エベレストの高さ)となります。
タイル座標が指定したCOGの右下から少しはみ出しているため、透明の領域ができていますが標高で色付けしたタイルを得ることができました。
全体的に真っ青なのは0 ~ 8848mの範囲で色付けしているせいで高さの差がわかりにくくなっているためです。
クエリを変えてmax=3000にしてみます。
海岸線や山の稜線が先程よりもわかりやすくなりました(濃い赤のあたりは富士山山頂です)。
mosaic
いい感じにタイルを動的生成できましたが、このままではクエリでCOGを指定しないといけないため、範囲が変わるたびにパラメータのファイル名を変更しなければいけません。
Openlayersなどで表示する際に、それでは1つのCOGにつき1つのレイヤーを作成しなければならず不便です。
DSMとしてレイヤーが1つで完結するようにrio-tilerのmosaic機能(mosaic_reader)を使います。
mosaic_readerは事前に配列でCOGのパスを指定しておくと、それらをつなぎ合わせて1つの巨大なCOGとして扱うことができます。
GDALでGeoTIFFをつなぎ合わせたVRTファイルを作るのと同じような働きです[4]。
from rio_tiler.mosaic import mosaic_reader
# COGパス
assets = [
"COG/ALPSMLC30_N035E138_DSM.tif",
"COG/ALPSMLC30_N035E139_DSM.tif",
"COG/ALPSMLC30_N035E140_DSM.tif",
"COG/ALPSMLC30_N036E138_DSM.tif",
"COG/ALPSMLC30_N036E139_DSM.tif",
"COG/ALPSMLC30_N036E140_DSM.tif"]
def tiler(src_path: str, *args, **kwargs):
with COGReader(src_path) as cog:
return cog.tile(*args, **kwargs)
@app.get("/dsm/{z}/{x}/{y}", **tile_routes_params)
def dsm(
z: int,
x: int,
y: int,
min: Union[int, float] = Query(0, description="Minimum Height."),
max: Union[int, float] = Query(8848, description="Maximum Height."),
):
if -9999 >= min:
min = 0
if min >= max:
max = 8848
(tile, mask), _ = mosaic_reader(assets, tiler, x, y, z)
if tile is None:
raise HTTPException(status_code=404, detail="Tile Outside Bounds.")
mask[tile[0] == -9999] = 0
tile[tile < min] = min
tile[tile > max] = max
tile = (tile/(max - min) * 255).astype(np.uint8)
format = ImageType.png
driver = drivers[format.value]
options = img_profiles.get(driver.lower(), {})
colormap = cmap.get('jet')
img = render(tile, mask, img_format=driver, **options, colormap=colormap)
return TileResponse(img, media_type=mimetype[format.value])
浸水域を表示するエンドポイントも作ってみます。
water = np.array([
76 * np.ones((256, 256), dtype=np.uint8),
108 * np.ones((256, 256), dtype=np.uint8),
179 * np.ones((256, 256), dtype=np.uint8)])
@app.get("/flood/{z}/{x}/{y}", **tile_routes_params)
def flood(
z: int,
x: int,
y: int,
depth: Union[int, float] = Query(0, description="Water Depth."),
):
if -9999 >= depth:
depth = 0
(tile, mask), _ = mosaic_reader(assets, tiler, x, y, z)
if tile is None:
raise HTTPException(status_code=404, detail="Tile Outside Bounds.")
mask[tile[0] == -9999] = 0
mask[tile[0] > depth] = 0
format = ImageType.png
driver = drivers[format.value]
options = img_profiles.get(driver.lower(), {})
img = render(water, mask, img_format=driver, **options)
return TileResponse(img, media_type=mimetype[format.value])
depthで指定した高さ以下の領域を青色で塗りつぶします。
Openlayersで表示
Openlayersのサンプルを元に表示するHTML+JSを作成します。
<!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>
<dl>
<dt>DSM</dt>
<dd>
<div>
<label><input type="checkbox" name="dsm_display" id="dsm_display" value="show">表示</label>
</div>
<div>
不透明度 <input type="number" name="dsm_opacity" id="dsm_opacity" min="0" max="100" value="100">
</div>
</dd>
</dl>
<dl>
<dt>浸水域</dt>
<dd>
<div><label><input type="checkbox" name="flood_display" id="flood_display">表示</label>
</div>
<div>
不透明度 <input type="number" name="flood_opacity" id="flood_opacity" min="0" max="100" value="100">
</div>
</dd>
</dl>
<script type="text/javascript">
var queries = window.location.search.slice(1).split('&');
var params = {};
for (var i = 0; i < queries.length; i++) {
var kv = queries[i].split('=');
params[kv[0]] = kv[1];
}
var max = params.hasOwnProperty('max') ? parseFloat(params.max) : 3000;
var depth = params.hasOwnProperty('depth') ? parseFloat(params.depth) : 10;
var extent = ol.proj.transformExtent([138, 35, 141, 37], 'EPSG:4326', 'EPSG:3857');
var dsm = new ol.layer.Tile({
extent: extent,
visible: false,
source: new ol.source.XYZ({
attributions: "提供：JAXA",
maxZoom: 12,
url: "http://127.0.0.1:8000/dsm/{z}/{x}/{y}?max=" + max,
crossOrigin: 'anonymous'
})
});
var flood = new ol.layer.Tile({
extent: extent,
visible: false,
source: new ol.source.XYZ({
attributions: "提供：JAXA",
maxZoom: 12,
url: "http://127.0.0.1:8000/flood/{z}/{x}/{y}?depth=" + depth,
crossOrigin: 'anonymous'
})
});
var map = new ol.Map({
target: 'map',
layers: [
new ol.layer.Tile({
source: new ol.source.OSM()
}),
dsm,
flood,
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([139.5, 36]),
zoom: 8
})
});
document.getElementById('dsm_display').addEventListener('change', function (event) {
dsm.setVisible(event.target.checked)
});
document.getElementById('dsm_opacity').addEventListener('change', function (event) {
dsm.setOpacity(parseFloat(event.target.value) / 100)
});
document.getElementById('flood_display').addEventListener('change', function (event) {
flood.setVisible(event.target.checked)
});
document.getElementById('flood_opacity').addEventListener('change', function (event) {
flood.setOpacity(parseFloat(event.target.value) / 100)
});
</script>
</body>
</html>
openlayers.html?max=1500&depth=15のようにクエリで最大の標高と浸水域の高さを指定できるようにしています(デフォルトは3000mと10m)。
それぞれのレイヤーの表示/非表示や不透明度はinputで制御できます。
まとめ
rio-tilerのmosaic_readerと標高データのCOGを使うと標高地形図を動的に、しかもクエリでパラメータを指定して生成できます。
今回はダウンロードしたGeoTIFFを自分でCOG化してローカルに配置しましたが、標高データのCOGはとても便利なものなのでJAXAやどこかの機関が(zipでなく)公開してくれないものでしょうか。
Discussion