rio-tilerを使って標高地形図を動的に生成する

10 min read

はじめに

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(エベレストの高さ)となります。

http://127.0.0.1:8000/dsm_test/9/453/202?url=COG/ALPSMLC30_N035E138_DSM.tif

タイル座標が指定したCOGの右下から少しはみ出しているため、透明の領域ができていますが標高で色付けしたタイルを得ることができました。
全体的に真っ青なのは0 ~ 8848mの範囲で色付けしているせいで高さの差がわかりにくくなっているためです。

クエリを変えてmax=3000にしてみます。

http://127.0.0.1:8000/dsm_test/9/453/202?url=COG/ALPSMLC30_N035E138_DSM.tif&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でなく)公開してくれないものでしょうか。

脚注
  1. フリーで利用できるDSMは、AW3D30以外にASTER GDEMというものもあります。 ↩︎

  2. GDALの環境が手元にない人はDockerを使うことをオススメします。 ↩︎

  3. 本来ならapp.pyとは別のサーバに置いてURLでアクセスするほうがデモとしてはいいのですが、サーバを建てるのがめんどうだったので今回は手を抜きます。 ↩︎

  4. 今回はapp.pyと同じサーバにCOGを置いているためmosaic_readerを使わずVRTでも事足りますが、COGが外部のURLで公開されているときmosaic_readerが真価を発揮します。 ↩︎

Discussion

ログインするとコメントできます