ラスター画素値演算の例 - PNG標高タイル方式による基盤地図情報DEMの表現

2016年6月7日に開催された「第5回地理院地図パートナーネットワーク会議」（国土地理院主催）では、地理院タイルの利用に関するさまざまな事例発表が行われましたが、「新シームレス地質図3Dと標高データ」（産業技術総合研究所）で紹介のあった「PNG標高タイル」は、ウェブブラウザで利用するための軽量なDEMラスター仕様として大変興味深く拝聴しました。発表資料は地理院地図パートナーネットワーク会議ホームページで公開されているので、ご参照ください。

これまでも何度かラスターの画素値演算を含むワークスペース例を掲載したことがありますが、ここでは画素値演算に焦点を絞り、基盤地図情報数値標高モデルデータを「PNG標高タイル」の画素値の求め方に準じたPNGラスターに変換するワークスペース例を掲げます。

FME 2016.1.0.1 build 16494

FMEワークスペース例

[JP_FGD_DEM] リーダー: 基盤地図情報DEM (XML) をラスターに変換して読み込む。
RasterBandNodateRemover: ラスターから Nodata の定義を削除。
RasterExpressionEvaluator: 画素値 (m単位標高, Real64) を100倍したラスター (cm単位標高, Int32) に変換。
RasterExpressionEvaluator_2: 「PNG標高タイル」に準じた RGBA32 ラスターに変換。
[PNGRASTER] ライター: PNG形式のファイルに出力。

[JP_FGD_DEM] リーダーは、FME Hub (旧 FME Store) で公開されているカスタムフォーマットで、基盤地図情報数値標高モデル (XML) をラスターとして読み込みます。海域などの標高データがない場所の画素値は -9999 とし、Nodata 値として定義しますが、演算を行う前に RasterBandNodataRemover によって Nodata の定義を削除しました。

ひとつ目の RasterExpressionEvaluator では各画素値を100倍するとともにバンドのデータ型をInt32 (整数) に変更し、ふたつ目の RasterExpressionEvaluator で「PNG標高タイル」の方式で RGBA32 (1バイト x 4バンド) ラスターに変換しました。各バンドの演算内容は以下のとおりです。

R: @if(A[0]<0, 255, @floor(A[0]/(256*256))
cm単位標高が負ならば255、0以上ならば (256 x 256) で除して小数部を切り捨てた値。

G: @if(A[0]<0, 255, @int32(@floor(A[0]/256))%256)
cm単位標高が負ならば255、0以上ならば 256 で除して小数部を切り捨てた値を 256 で除した余り。

B: @if(A[0]<0, 255, A[0]%256)
cm単位標高が負ならば255、0以上ならば 256 で除した余り。

A: @if(A[0]==-999900, 0, 255)
cm単位標高が -999900 (元の Nodata 値 x 100) に等しければ 0 (透過)、それ以外は 255 (不透過)。

各式における A[0] は、A ポートから入力されたラスターの第1バンドの画素値を示します。

基盤地図情報10mメッシュDEM, 富士山周辺二次メッシュ22区画の変換結果
注: 掲載の都合上、データサイズを小さくするために1/10に縮小して画素値をリサンプルしたため、この画像の画素値から元の標高値を正確に再現することはできません。

要は、cm単位の標高値(整数)を256進数に基数変換した3桁で表現し、最上位桁から順番に R、G、B バンドに割り当てるとともに、Nodata だった画素の Alpha バンドに 0 を割り当てることにより透明にしたということです。m単位の標高値は次の式によって求められます。

h = (R x 256 x 256 + G x 256 + B) / 100
ただし A = 0 の場合は h = -9999 (Nodata)

2016-06-09 補足: 前述の画素値演算の式とこの式は、負の標高値を考慮していない簡略版です。産業技術総合研究所の「PNG標高タイル」の仕様は、RGB3バンドを24ビット符号付き整数として解釈することにより、負の標高値も表現できるようになっています。「PNG標高タイル」の仕様については次のページを参照してください。
産業技術総合研究所 地質情報研究部門 シームレス地質情報研究グループ: シームレス標高サービス実験公開 > PNG標高タイル ( PNG Elevation tile )

元のラスターが倍精度浮動小数点数1バンド (1画素あたり64ビット) であったのに対して、「PNG標高タイル」の方式では1バイト符号なし整数4バンド（1画素あたり32ビット）、Alpha バンドを不要とすれば3バンド（1画素あたり24ビット）で表現できることになり、広い意味でのデータ圧縮とも言えます。

この方式では、RGB3バンドを24ビット符号付き整数として解釈するならば -2²³ (-8,388,608) 以上 2²³ (8,388,608) 未満、24ビット符号なし整数として解釈するならば 0 以上 2²⁴ (16,777,216) 未満の値が表現できるので、数値標高モデルだけでなく、さまざまなデータの表現に応用できるかも知れません。

なお、式を分かり易くするためワークスペース例では RasterExpressionEvaluator をふたつ使いましたが、次のパラメーター設定によってひとつにまとめることもできます

==========
2016-06-10 追記: 産業技術総合研究所「PNG標高タイル」仕様にしたがい、RGB3バンドを24ビット符号付き整数として解釈するとともに、Nodata 画素値を 0,0,0,0 (黒, 透過) とする場合の式は次のとおりです（標高分解能=0.01m）。

RasterExpressionEvaluator: 1バンドInt32ラスターに変換
Int32: @if(A[0]<0, @if(A[0]==-9999, -1, @int32(A[0]*100)+16777216), @int32(A[0]*100))
標高値が負の場合、Nodata (-9999) ならば -1 (負)、そうでなければ 16,777,216 (= 2²⁴) を加算することにより24ビットでの2の補数 (負値の表現) とする。

RasterExpressionEvaluator_2: RGBA32ラスターに変換
R: @if(A[0]<0, 0, @floor(A[0]/(256*256)))
G: @if(A[0]<0, 0, @int32(@floor(A[0]/256))%256)
B: @if(A[0]<0, 0, A[0]%256)
A: @if(A[0]<0, 0, 255)

このRGBA32ラスターの画素値から元の標高値(m)を求める式は次のようになります。

x = R * 256 * 256 + G * 256 + B
R < 128 の場合: h = x * 0.01
R >= 128 の場合: h = (x - 2²⁴) * 0.01
ただし A = 0 の場合は h = -9999 (Nodata)