PNG標高タイルからDEMラスターへの変換

「ラスター画素値演算の例 - PNG標高タイル方式による基盤地図情報DEMの表現」でRGB画像ラスターによるDEM - 数値標高モデルの表現方法を取り上げましたが、それと同じ方式で作成された基盤地図情報数値標高モデルのPNG標高タイルが公開されています。

国土地理院 | 地理院タイル一覧 | 標高タイル（基盤地図情報数値標高モデル）

前回は、DEMラスターをPNG標高タイル方式のRGB画像ラスターに変換しましたが、今回はその逆に、地理院地図のPNG標高タイルをDEMラスターに変換してみます。どちらも RasterExpressionEvaluator による画素値演算の分かりやすい例だと思います。

FME 2017.0.0.1 build 17271

FMEワークスペース例

JpGsiTileFetcher (カスタムトランスフォーマー): PNG標高タイルを取得
RasterExpressionEvaluator: PNG標高タイル (RGB24) をUInt64ラスターに変換
RasterExpressionEvaluator_2: UInt64ラスターをDEMラスター (Real64) に変換
RasterBandNodataSetter: Nodata値 (-9999) を設定
RasterMosaicker: 複数のタイルラスターを結合

JpGsiTileFetcher を使うと、ワークスペースの実行時に、入力フィーチャーのジオメトリをカバーする範囲について、パラメーターで指定したズームレベル、種類の地理院タイルデータを取得することができます。このカスタムトランスフォーマーは、次のページで公開しています。

FMEサポート | 国内データ変換のサポート

このカスタムトランスフォーマーのパラメーターを次のように設定すれば、入力フィーチャーのジオメトリをカバーする範囲のズームレベル 14 のPNG標高タイルが取得できます。もちろん、実行時にはインターネットに接続していなければならず、また、入力フィーチャーには座標系が適切に設定されている必要があります。

• GSI Tile Data ID: dem_png ※基盤地図情報DEM10B PNG標高タイルのデータID
• Format Extension: png
• Zoom Level: 14

三宅島の区域を例として、上記ワークスペースによってPNG標高タイルを取得、DEMラスターに変換した結果は次のとおりです。

左: 入力フィーチャーのジオメトリ (三宅島の区域), 中: PNG標高タイル, 右: DEMラスター
掲載用に画像サイズを縮小したので、解像度は実際のPNG標高タイルやDEMラスターとは異なります。

地理院地図 | 標高タイルの詳細仕様によれば、PNG標高タイルの画素値 (RGB値) から標高値 h (m) への換算式は次のとおりです。

x = 2¹⁶R + 2⁸G + B ... (1)
x < 2²³ の場合　h = xu ... (2)
x = 2²³ の場合　h = NA ... (3)
x > 2²³ の場合　h = (x-2²⁴)u ... (4)
u: 標高分解能 (0.01m)

RasterExpressoinEvaluator パラメーター

Interpretation	Expression
UInt64	@pow(2,16)*A[0] + @pow(2,8)*A[1] + A[2]

最初の RasterExpressionEvaluator では、PNG標高タイル = 8ビット符号なし整数 x 3バンド (R, G, B) ラスター を、(1) 式によって求めた x を画素値とする UInt64 (64ビット符号なし整数) x 1 バンドラスターに変換しました。

RasterExpressionEvaluator_2 パラメーター

Interpretation	Expression
Real64	@if(A[0]==@pow(2,23), -9999.0, @if(A[0]<@pow(2,23), A[0]*0.01, (A[0]-@pow(2,24))*0.01))

2番目の RasterExpressionEvaluator では、上記変換後の UInt64 ラスター (画素値 = x) に基づき、(2), (3), (4) 式によって標高値を求め、Real64 (64ビット浮動小数点数) x 1 バンドラスターに変換しました。(3) 式の NA は値がないことを意味していますが、ラスターの画素にはなんらかの値を設定する必要があるので、標高としてあり得ない値: -9999 を設定しました。この値は、後で RasterBandNodataSetter によって Nodata値として取り扱うように設定します。

式 (Expression) 中の A[n] は、入力ラスターの n 番目 (n は 0 から始まる連番) のバンドの画素値を示します。最初の RasterExpressionEvaluator の入力ラスターはRGB3バンドであり、式中の A[0], A[1], A[2] が R, G, B の値を示しています。2番目の RasterExpressionEvaluator の入力ラスターは (1) 式で求めた x の値を格納している 1 バンドラスターなので、A[0] (x の値) のみを使用しています。

@pow(a, b) 関数は a の b 乗を返し、@if(a, b, c) 関数は a が真のときは b、偽のときは c を返します。換算式と対比し易いように、2 の累乗を行うのに @pow 関数を使いましたが、2ⁿ の値を定数として式の中に記述しても構いません。例えば、@pow(2,16) は 65536、@pow(2,8) は 256 に置き換えられます。

RasterBandNodateSetter で -9999 を Nodata値として設定し、RasterMosaicker でタイルごとのラスターを結合してひとつのラスターに変換すればできあがりです。

ワークスペース例には示していませんが、DEMラスターを格納するのに適したフォーマット (GeoTIFF, ERDAS IMAGINE など) のライターを追加し、RasterMosaicker の後に接続すれば、結果をファイルに出力することができます。