HTMLレポートページの作成 - 東京都知事選挙投票結果

FMEの文字列処理用のトランスフォーマーは充実しており、「e-Stat 統計表情報ウェブページの自動作成 (FME Server)」でも紹介したように、従来から、工夫次第でHTML文書を作成することが可能でしたが、FME 2016.1 ではHTMLに関していくつかの機能が追加され、図表を中心とする定型のHTMLレポートページが簡単に作成できるようになりました。

次の Safe Software ブログ記事も参照してください。
FME 2016.1 New Functionality: HTML and HTML Transformers

ここでは、2016年7月31日に投開票が行われた東京都知事選挙の投票結果について、開票区別の投票率を示す棒グラフと開票区別の投票結果一覧表を含むHTMLレポートページを作成するワークスペース例を掲げます。

【ソースデータ】
東京都選挙管理委員会事務局ウェブサイトにおいて公開されている東京都知事選挙投票結果Excelファイル
東京都選挙管理委員会事務局 > 都知事選挙（平成28年7月31日執行）投開票結果
投票結果: h28tochiji_touhyo.xls [Sheet1]

【変換内容】
開票区別の投票率（男女平均）を示す棒グラフ、および、開票区別の有権者数（男、女、計）、投票者数（男、女、計）、投票率（男、女、平均）の一覧表をウェブブラウザで閲覧できるHTML文書に変換する。

FME 2016.1.1.0 build 16609

FME ワークスペース例

[XLSXR] (Excel) リーダー: Excelスプレッドシートから投票結果を読み込む。
AttributeRenamer: 属性名（列名）を分かり易い名称に変更する。
AttributeTrimmer: 開票区名の先頭にある余分な全角空白を削除する。
Tester: 集計行（開票区名が「***計」、「***支庁」であるもの）を除く。
StringFormatter x 2: 投票者数、投票率の数値書式調整
HTMLReportGenerator: HTML図表作成
HTMLlayouter: HTMLページのレイアウト調整
[HTML] ライター: HTMLファイル出力

Excelリーダーをワークスペースに追加する際のパラメーター設定によって、列名が記述されている行番号を指定し、それらを属性名とすることができますが、その行に同一の列名があった場合、2番目以降の列名には連番 (00, 01, 02...) が自動的に付加されます。この例では、6行目を列名としたので、2番目以降の「男」は、「男00」、「男01」、「男02」となります。ただし、3番目の「男01」（棄権者数）は使用しないので、フィーチャータイププロパティ設定画面の User Attributes (ユーザー属性) タブで非表示に設定したため、上の図には現れていません。「女」等についても同様です。

Excelに数値として記録されている値には整数と小数点数の区別はなく、Excel画面上に表示される小数点以下の桁数は、セルの書式として設定されています。FME はセルの書式のいかんに関わらず、Excelの数値データを全て小数点数として読み込むので、HTML文書を作成する前に StringFormatter によって数値の書式を整えました（有権者数、投票者数は小数部なしの整数、投票率は小数点以下2桁の小数点数）。

以上の前処理の後、FME 2016.1 で導入された HTMLReportGenerator, HTMLLayouter, [HTML] ライターによって、次のようなHTML文書=ウェブページに変換しました（実際のページはこちら > HTMLReportGenerator DEMO）。

HTMLReportGenerator では、入力データを次の要素を含むHTML文書に変換することができます。

Chart (Bar, Line, Pie)	チャート (棒グラフ, 折れ線グラフ, 円グラフ)
Header	ヘッダー（H*要素文字列)
List	リスト
Image	画像
Map	マップ (注)
Separator	セパレーター (hr要素)
Table	表
Custom HTML (FME 2017.0 以降)	任意のHTML要素を含むHTML文書のパーツ

注: マップ要素は、Esri Leaflet, Google, または Mapbox Leaflet のAPIを利用して地物を表示するためのスクリプトをHTML文書に組み込むものです。利用にあたっては各APIの利用規約等にしたがってください。

HTML文書に組み込む要素の種類、数、順序は任意です。この例では、ヘッダーx3, セパレーター, チャート (棒グラフ), ヘッダー, 表の順で使用しました。

HTMLReportGenerator パラメーター設定画面

HTMLLayouter では、HTMLReportGenerator によって作成したひとつ以上のHTML文書を1ページのHTML文書に統合し、それらのレイアウトを整えることができます。この例では HTMLReportGenerator でひとつの文書しか作成していませんが、 Layout Type パラメーターを Vertical (垂直) とすることにより、左右のマージンなどが整えられます。

HTMLLayouter パラメーター設定画面

[HTML] ライターは、HTMLReportGenerator, HTMLLayouter で作成されたHTML文書（html_content属性に格納された文字列）をそのままファイルに出力します。

[HTML] ライターの代わりに FeatureWriter トランスフォーマーによってファイル出力すれば、引き続き FTPCaller トランスフォーマー等によってウェブサーバーにアップロードするところまでをひとつのワークスペースで自動化することもできます。

HTMLReportGenerator, HTMLLayouter は、特に、静的ではあるが、しばしばデータの更新があるようなウェブページついて、データ更新に伴うHTML文書の更新を正確かつ迅速に行うために活用できそうです。

なお、現時点 (FME 2016.1.1) の HTMLReportGenerator はまだ機能が限られており、思い通りのウェブページを作成するのが難しいケースもあります。例えば:

• 表の各欄の右揃え、センタリング等の設定をサポートしておらず、全て左揃えになる。
• 図の凡例文字列中のスペースや非ASCII文字の前で改行され、縦書きのように表示される。
• 入力文字列中の <, > 等の特殊文字が自動的に実体参照に置き換えられるため、任意のHTML要素（例えば他のページへのリンクなど）を埋め込むことができない (FME 2017.0 では "Custom HTML" によって、任意のHTML要素を含むHTML文書のパーツを組み込むことができるようになりました)。

これらについてはすでに多数のユーザーから改良、機能強化の要望があがっているので、今後のアップグレードで改良されていくことと思われます。

GISポリゴンデータに基づく面積計算 (_AZMEA_の利用)

GISポリゴンデータに基づいて地物の面積を求めることはよくあると思います。GISソフトウェア等の機能として面積計算はありふれたものですが、データの精度は元の地図の縮尺やデジタイズの方法等によって異なるとともに、同じデータであっても投影先の座標系によって地物の形状は異なるため、使用するデータや投影方法によって、同一の地物について求められる面積が異なることもあるということは認識しておく必要があるでしょう。

既存のGISデータから求めた面積を研究論文・報文や調査業務報告書などで用いるのならば、データの出典を明記するのはもちろんのこと、元のデータの座標系が緯度経度である場合には、どの座標系に投影した形状の面積なのかも説明するべきであると考えます。

緯度経度の座標系で作成された日本国内のGISデータに基づいて地物の面積を求める場合、平面直角座標系またはUTM座標系に変換することが多いかも知れません。

しかし、平面直角座標系の複数の系、UTM座標系の複数のゾーンにまたがった広い範囲に分布する地物の面積を求める場合には、個々の地物が所在する位置に応じて変換先座標系の系番号、ゾーン番号を選択するといったプロセスが必要になります。そのようなプロセスを含むワークスペースを作成することもできますが、合理的な説明ができれば面積計算にあたっての投影法は任意で良いという場合には、_AZMEA_ (Dynamic Reprojection Equal Area) という名前の座標系を利用するのが便利です。

_AZMEA_ は、任意の座標を基準点とする「ランベルト正積方位図法」 (Lambert azimuthal equal-area projection) によって定義される座標系 (単位: メートル) であり、座標系変換を行う Reprojector または CsmapReprojector トランスフォーマーで変換先の座標系として _AZMEA_ を指定すると、ワークスペースの実行時に、個々の地物について元の座標系におけるバウンディングボックス中心を基準点とするランベルト正積方位図法の座標系が定義され、それに変換されます。これにより、系番号やゾーン番号を考慮しなくても、広い範囲に分布する地物の面積を統一的な手法で求めることができます。

ここでは、国土数値情報「行政区域」データに基づいて、_AZMEA_ を利用して全国の市区町村の面積を求めるワークスペース例を掲げます。

【ソースデータ】
国土数値情報ダウンロードサービスサイトからダウンロードした「行政区域」データ（平成28年1月1日時点, 全国, Esri Shapefile 形式ポリゴンデータ）

FME 2016.1.1.0 build 16609

FMEワークスペース例 1

[ESRISHAPE] リーダー: 国土数値情報行政区域（全国）データを読み込む。
Tester: 市区町村コード（N03_007属性値）を持っているフィーチャーとそうでないものに分岐する。
Aggregator: 同じ市区町村コードを持っているフィーチャーを集約する。
Reprojector: フィーチャーごとに _AZMEA_ に座標系を変換する。
AreaCalculator: フィーチャーごとに面積を計算する。
Inspector: 変換結果を FME Data Inspector で確認する。

国土数値情報の行政区域データ（ポリゴン）は、飛び地や離島などがある市区町村について、それらが個別のレコードとなっている、いわゆるシングルパートポリゴンです。また、埋立地などでどの市区町村に属するか決定していない「所属未定地」のレコードには、市区町村コードが記録されていません。

そのため、このワークスペース例ではまず、Tester によって市区町村コードを持っているフィーチャーと持っていないフィーチャー（所属未定地）に分け、市区町村コードを持っているフィーチャーについて Aggregator によって市区町村ごとに集約しました。これによって、飛び地や離島がある市区町村についても、市区町村単位で全域の面積を求めることができます。

Reprojector の Destination Coordinate System （変換先の座標系）パラメーターに _AZMEA_ を設定することにより、前述のとおり、フィーチャー（市区町村または所属未定地）ごとに、それぞれバウンディングボックス中心を基準点とする「ランベルト正積方位図法」による座標系に変換されます。

Reprojector パラメーター設定画面

AreaCalculator は、入力フィーチャーの面積を計算し、その値を格納する新たな属性を付加します。ここで、Multiplier パラメーターによって、求められた値に乗じる値を指定することができます。_AZMEA_ によって変換された座標値の単位はメートルなので、Multiplier パラメーターに 1e-6 (= 0.000001) を設定すれば、平方キロメートル単位の値に換算されます。


このワークスペースを実行して Data Inspector に変換後のポリゴンを表示させると分かりますが、_AZMEA_ は、個々のフィーチャーのバウンディングボックス中心を基準点、つまり原点 (0, 0) とする座標系なので、変換後のポリゴンは原点周辺で重なり合うことになります。地図としては用を成しません。

面積を求めることだけが目的ならばこのままでも構いません。しかし、元の座標系（この場合は JGD2011 緯度経度）におけるジオメトリとしての処理も必要な場合には、元の座標系における位置、形状を復元する必要があります。

もうひとつ Reprojector を追加して元の座標系に再変換するのは簡単ですが、座標系変換にはどうしても計算誤差が伴い、完全に復元できるとは限りません。正確に元に戻すには、あらかじめ元の座標系名とジオメトリを属性データとして抽出しておき、面積を求めた後で、それらの属性データに基づいて元のジオメトリと座標系の定義を復元するという方法があります。

上記ワークスペースにそれらの手順を追加すると、次のようになります。

FME ワークスペース例 2

(座標系変換前）
CoordinateSystemExtractor: フィーチャーに設定されている座標系名を属性 (_coordsys) として抽出する。
GeometryExtractor: ジオメトリ（形状のデータ）を属性 (_geometry) として抽出する。

（面積計算後）
GeometryReplecer: 属性 (_geometry) に基づいてジオメトリを復元する。
CoordinateSystemSetter: 属性 (_coordsys) に基づいて元の座標系の定義を復元する。

GeometryExtractor/Replacer では、Geometry Encoding パラメーターでいくつかの属性値（文字列）としてのジオメトリの表現方法が選択できますが、そのうち FME Binary は、ジオメトリの抽出/復元に伴う誤差は全くなく、かつ、最も効率が良いので、この例のようにジオメトリを復元することを目的とする場合には最も適しています。

GeometryReplacer は、フィーチャーのジオメトリ（形状のデータ）だけを置き換えるものであり、座標系の定義は復元しないことに注意してください。そのため、このワークスペース例では、CoordinateSystemExtractor によって元の座標系名をあらかじめ抽出しておき、CoordinateSystemSetter によってその座標系の定義を復元しています。

なお、GeometryExtractor/Replacer は、FME Binary のほか、GeoJSON, GML, OGC Well Known Binary, OGC Well Known Text などのメジャーなジオメトリの文字列表現をサポートしており、用途に応じて使い分けることができます。

MapnikRasterizer による画像作成 - 基盤地図情報等高線と標高点

複数のベクター、ラスターレイヤを重ね合わせ、ひとつの地図として視覚的に表現することは、GISが得意とするところです。FMEはGISではありませんが、MapnikRasterizer トランスフォーマーによって、GISにも匹敵する高度な地図画像（ラスター）を作成することも可能です。

名前からも分かるように、このトランスフォーマーは Mapnik というオープンソースの地図レンダリング用ツールキットを内部で使用することにより、非常に多彩な表現をサポートしています。

ここでは、基盤地図情報基本項目「等高線」および「標高点」のレンダリングを例として、MapnikRasterizer の基本的な使い方を示すワークスペース例を掲げます。

【ソースデータ】
三宅島周辺の二次メッシュ4区画（513903, 513904, 513913, 513914）について、基盤地図情報ダウンロードサービスサイトからダウンロードした基盤地図情報基本項目の等高線および標高点データ（GML形式）

【変換内容】
等高線をライン、標高点の位置をマーカー、標高点の標高値をテキストとしてレンダリングする。
等高線は「100m間隔」、「50m間隔」、「その他」に分類し、それぞれ異なる色で表現する。
変換結果は、三次メッシュ区画単位でPNG形式のファイルに出力する。

FME 2016.1.1.0 build 16609

FME ワークスペース例

[GML] リーダー: 基盤地図情報等高線、標高点データ（GML形式）を読み込む。
JpStdGridAccumulator: 入力フィーチャーをカバーする範囲の三次メッシュ区画ポリゴンを作成する。
BoundsExtractor: 各メッシュ区画の東西端、南北端の座標を抽出する。
Clipper: メッシュ区画ポリゴンで等高線、標高点をクリップする。
GeometryFilter: ラインフィーチャー（等高線）とポイントフィーチャー（標高点）に振り分ける。
TestFilter: 等高線を「100m間隔」、「50m間隔」、「その他」の3区分に振り分ける。
MapnikRasterizer: メッシュ区画ごとに等高線（ライン）、標高点（マーカー、テキスト）をレンダリングする。
[PNGRASTER] ライター: メッシュ区画ごとにラスターをPNG形式のファイルに出力する。

基盤地図情報基本項目データはGML形式で提供されており、FMEは基盤地図情報応用スキーマをバンドルしているので、標準のGMLリーダーで直接読み込むことができます。

JpStdGridAccumulator は FME Hub で公開しているカスタムトランスフォーマーで、パラメーターの設定に応じて、入力フィーチャーをカバーする範囲の一次～1/10細分メッシュ区画ポリゴンまたはポイントを作成するとともに、それらにメッシュコードを属性として付加します。この例では、三次メッシュ区画ポリゴンを作成しました。

BoundsExtractor は、各フィーチャーのバウンディングボックス境界の座標を抽出します。この例では各三次メッシュ区画の境界の座標が得られ、それらによって後述の MapnikRasterizer で画像の作成範囲を指定しています。

Clipper によって、等高線と標高点をメッシュ区画ポリゴンでクリップしました。このとき、Clipper のパラメーター設定画面で Merge Attributes（属性を結合する）オプションをオンにすることにより、メッシュコードやバウンディングボックス境界座標を格納したメッシュ区画ポリゴンの属性が、区画ごとにその内部となる全ての等高線、標高点に結合されます。

GeometryFilter でフィーチャーの流れをライン（等高線）とポイント（標高点）に分岐した後、等高線についてはさらに、TestFilter によって、各フィーチャーが持っている標高値（alti属性の値）に応じて「100m間隔」、「50m間隔」、「その他」の3区分に振り分けました。

MapnikRasterizer には、必要に応じて任意の数の入力ポートを追加することができ、ポートごとに、そこから入力されるフィーチャーに対するひとつ以上のレンダリングルール（視覚的な表現方法）を設定することができます。

レンダリングルールの設定内容は多岐にわたります。主なものだけ挙げると、ラインについては線幅、色、ダッシュスタイル（実線部の長さ、間隔の任意の繰り返しパターン）など、マーカーについてはサイズ、色、あるいは外部画像ファイルなど、テキストについてはフォント、サイズ、色、ポイントに対する相対的な位置などです。

MapnikRasterizer パラメーター設定画面

ソースデータ（二次メッシュ4区画）のうち、等高線、標高点が存在する三次メッシュ区画数は70でした。そのうち、雄山（おやま）火口周辺の次の9区画（下記）の画像をモザイク（結合）した結果を示します。
5139-14-11, 5139-14-12
5139-14-01, 5139-14-02
5139-04-91, 5139-04-92
この画像の解像度は、上記 MapnikRasterizer での設定と同じ（三次メッシュ1区画あたり1200 x 800ピクセル）です。ただし、ファイルサイズを小さくするため、パレット（エントリー数=16）つきの1バンドラスターに変換しました。ピクセル等倍でご覧になるには、ダウンロードして適当なビューアーで開いてください。

一見、うまくいっているようですが、上記ワークスペース例には、メッシュ区画境界付近で標高点テキストが描画されないことがあるという欠陥があります。上の図では、画像下部、三次メッシュ区画 5139-04-91 と 5139-04-92 の境界付近の 5139-04-91 側、つまり、西側の区画の東端付近にある標高点について、標高値テキストが描画されていません。これは、MapnikRasterizer では、指定した画像作成範囲をはみ出すようなテキストはレンダリングされないためです。

これは、ワークスペースに次のような改良を加えることによって解消することができます。
・MapnikRasterizerによる画像作成範囲を三次メッシュ区画ぴったりの範囲でなく、隣接区画の標高点テキスト描画範囲を含むのに十分な範囲まで拡大する。
・作成するラスターの解像度（縦横のピクセル数）をその範囲にあわせて修正する。
・MapnikRasterizerで作成された画像を、対応する三次メッシュ区画の範囲によってクリップする。

この改良を行うために追加したり置き換えたりするトランスフォーマーは、Scaler, SpatialRelator, ListExploder, JpMeshCodeReplacer (FME Hub) などです。

==========
2016-09-08 追記: クリップではなく、RasterSubsetter によって三次メッシュ区画に相当する範囲のラスターに変換する方法もあり、場合によってはその方が効率が良いかも知れません。
==========

改良後のワークスペースで作成した三次メッシュ区画 5139-04-91, 5139-04-92 の画像をモザイクした結果は次のとおりです。改良前には描画できなかった中央付近の標高値テキスト「723」も描画されました。

MapnikRasterizer は、複数のベクター、ラスターデータをまとめてひとつの地図画像を作成するのに非常に便利なトランスフォーマーですが、現時点では、非ASCII文字を含むテキストのレンダリングについて、MS UI Gothic (Normal) 以外のフォントがサポートされていないという制約があります。

その他のフォント（MS ゴシック、MS 明朝など）での日本語文字の描画が必要な場合、TextAdder 等によってテキストフィーチャーを作成し、さらに TextStroker によってフォントを指定してポリゴンに変換したうえでレンダリングするという代替策はありますが、サイズや位置の調整が少し面倒になるかも知れません。近い将来、MS UI Gothic 以外の非ASCII文字フォントのサポートも追加されることを期待しています。

ポイントクラスによるLAS点群の分割

ASPRS LiDAR Data Exchange Format (LAS) フォーマットの点群データは、点の座標を示す x, y, z のほか、通常、レーザーパルスのリターン強度を示す intensity、地物等の種類 (ポイントクラス) を示す classification  などのコンポーネントを持っています。

ASPRS LAS バージョン 1.4 仕様書によれば classification コンポーネントの値は 0～255 の整数であり、各値の意味は次のように定義されています。なお、LAS バージョン 1.1, 1.2, 1.3 では 0～31 の値が使われ、また、その範囲の値でも一部、下表と意味が異なるものがあります。

Classicifation Value	Meaning	意味 (仮訳)
0	Created, never classified	分類不能
1	Unclassified	未分類
2	Ground	地表
3	Low Vegetation	植生(低)
4	Medium Vegetation	植生(中)
5	High Vegetation	植生(高)
6	Building	建築物
7	Low Point (noise)	低ポイント (ノイズ)
8	Reserved	予約
9	Water	水部
10	Rail	鉄道
11	Road Surface	路面
12	Reserved	予約
13	Wire - Guard (Shield)	電線 - 保護線
14	Wire - Conductor (Phase)	電線 - 導線
15	Transmission Tower	送電塔
16	Wire-structure Connector (e.g. Insulator)	電線接続構造物 (例: 碍子)
17	Bridge Deck	橋梁床版
18	High Noise	高ノイズ
19-63	Reserved	予約
64-255	User definable	ユーザー定義

注: "Classification Value" および "Meaning" は、ASPRS ウェブサイト掲載の "LAS SPECIFICATION VERSION 1.4 – R13 15 July 2013" - Table 17: ASPRS Standard LIDAR Point Classes (Point Data Record Formats 6-10) を引用したものです。「意味 (仮訳)」は、国内関連業界においてオーソライズされた訳語・用語ではなく、不適切なものがあるかも知れません。実務においては仕様書原文を参照するとともに、実データの内容に応じて適切に解釈してください。

FMEでは、PointCloudSplitter トランスフォーマーによって、任意のコンポーネントの値に基づいて点群データに含まれる個別の点を分類し、その分類区分ごとの点群に分割することができます。

ここでは、このトランスフォーマーを使い、単一の LAS 点群データを classification コンポーネントの値に基づいて「地表」、「植生(低)」、「植生(中)」、「植生(高)」、「建築物」の点群に分割し、それぞれ LAS ファイルに出力するワークスペース例を掲げます。

サンプルデータは、Safe Software 社 FME Server デモサイトの次のページで作成、ダウンロードしました。

FME Server Demos | Point Cloud Data Distribution

このページでは、カナダ・バンクーバー市の一部 (黒枠で示された4タイルの範囲) について、地図画面上で描画したポリゴンに含まれる点群データを作成、ダウンロードすることができます。以下のワークスペース例で使用したデータは、地図画面左下のタイルについて次の条件でダウンロードしたものです。
・Output Format: ASPRS Lidar Data Exchange Format (LAS)
・Output Coordinate System: NAD83 datum, UTM Zone 10
・Return Tiled Output: YES (タイル境界で分割する)
・Split Point Cloud: NO (コンポーネントの値による分割はしない)
・Use Point Cloud Thinning: NO (点の間引きはしない)

本題ではありませんが、 同サイトではこのほかにもいくつかの FME Server のデモページが公開されていますので、関心のある方はご覧ください。> FME Server Demos

FME 2016.1.0.1 build 16494


FMEワークスペース例

[LAS] リーダー: LAS形式点群データを読み込む。
PointCloudSplitter: classification コンポーネントの値 (クラス) に基づいて点群を分割する。
AttributeValueMapper: 各クラスに対応する意味 (文字列) 属性を作成する。
[LAS] ライター: 分割後の点群データをクラス別のLASファイルに出力する。

変換結果 FME Data Inspector による表示 (合成)
左: 分割前 (15,935,683 ポイント)
右: 分割後 上から 建築物 (3,435,083), 植生(高) (1,505,395), 植生(低) (2,839,508), 地表 (1,439,096)
注1: 分割前の点群には「未分類」等のポイントも含まれるため、分割後のポイント数の合計と分割前のポイント数は一致しません。
注2: サンプルデータには「植生(中)」 (classification = 4) に該当する点は含まれていませんでした。低、中、高の3区分で分類したならば「植生(中)」が全く現れないということは考えにくいので、このデータの植生の分類は低と高の2区分で行われていると考えるのが妥当と思われます。

PointCloudSplitter によって、任意のコンポーネントの値に基づいて点群を分割することができます。この例では、次のようにパラメーターを設定することにより、classification コンポーネントの値に応じて「地表」、「植生(低)」、「植生(中)」、「植生(高)」、「建築物」クラスの点群に分割しました (その他のクラスの点は破棄)。

Split By: classification (classification コンポーネントの値で分割する)
Split Type: Unique (個別の値ごとに分割する)
Output Attribute: _split_value (コンポーネントの値を _split_value 属性に格納して分割後の点群に付加する)
Unique Values to Keep: 2, 3, 4, 5, 6 (これらの値に該当する点のみを保持する)

PointCloudSplitter パラメーター設定画面

PointCloudSplitter だけでコンポーネントの値による点群の分割はできますが、classification コンポーネントの値 (整数) だけでは意味が分かりにくいので、次の AttributeValueMapper によって値に対応する意味 (文字列) を格納する _class 属性を作成しました。後述するように、この文字列を出力先のファイル名とすることを意図しています。

AttributeValueMapper では、入力フィーチャーの属性 (Source Attribute パラメーターで選択) の値に対して 1:1 または N:1 で対応する値を格納する新たな属性 (Destination Attribute パラメーターで属性名を指定) を作成することができます。値の対応関係は、Value Map パラメーター (リスト) で定義します。

AttributeValueMapper パラメーター設定画面

AttributeValueMapper で作成した _class 属性を [LAS] ライターフィーチャータイプのフィーチャータイプ名パラメーター (LAS File Name) に設定することにより、分割後の点群の出力先ファイル名を、「地表」などのクラス名とすることができます。

LAS ライターフィーチャータイププロパティ画面 (左) と出力先ファイル名 (右)

補足: 点群の分割は PointCloudFilter でもできます。このトランスフォーマーを使った場合はデータフローが条件別に分岐するので、点群を単に分割するだけでなく、同じワークスペース内で引き続き、分割後の複数の点群に対してそれぞれ異なる変換処理を施す場合に適しています。

PointCloudFilter パラメーター設定例
Expression で条件式、Output Port でその条件に該当する点で構成される点群の出力先ポート名が設定できます。

点群からDEMラスターへの変換

点群 (Point Cloud: ポイントクラウド) は複数の点の集合体を表すデータ構造であり、レーザー測量によって計測された地物の立体的な形状のデータを格納、処理する場合などでよく利用されます。

点群に含まれる個々の点はコンポーネント (components) と呼ばれる特性値を持っています。コンポーネントの種類はアプリケーションや目的に応じて異なりますが、点の座標を示す x, y, z コンポーネントはほとんどの点群データに存在し、レーザー測量の点群データならば多くの場合、x, y, z に加えてレーザーパルスのリターン強度を示す intensity、地物等の種類を示す classification などのコンポーネントが存在しています。また、色 (光の三原色) を示す color_red, color_green, color_blue コンポーネントを持つ点群データもあります。

FMEは ASPRS Lidar Data Exchange Format (LAS), Point Cloud XYZ などのメジャーな点群フォーマットをサポートするとともに、点群を加工、変換するための多数のトランスフォーマーを備えています。

点群とラスターの間の変換も容易であり、これまでラスターを点群に変換する処理を含むワークスペース例を紹介したことがありますが、ここでは逆方向の変換の例として、地表の標高を z コンポーネントに記録している点群データをDEMラスターに変換するワークスペース例を掲げます。

FME Knowledge Base > Creating a Raster from a Point Cloud の後半で説明されている内容とほぼ同じであり、ソースデータセットもこの記事で使用されているものと同じ LAS データを使用しました。

ソースデータセット: LAS 形式, FME Data Inspector による 3D 表示
この点群は地表 (classification = 2) の点のみを含んでおり、点間の間隔は密度が高いところで2m前後です。

FME 2016.1.0.1 build 16494


NumericRasterizer による変換

FMEワークスペース例 1

[LAS] リーダー: LAS形式の点群データを読み込む。
NumericRasterizer: 点群をラスターに変換する。
Inspector: 変換結果を FME Data Inspector に出力する。

NumericRasterizer では、X Cell Spacing, Y Cell Spacing パラメーターによって作成するラスターのセルサイズを指定できます。2m x 2m, 5m x 5m, 10m x 10m の3通りの設定で試した結果は次の図のとおりです。

変換結果: 左からセルサイズ 2m x 2m, 5m x 5m, 10m x 10m

x, y, z コンポーネントを持つ点群データを NumericRasterizer に入力した場合、XY平面に投影したときにセル内に存在する点の z コンポーネントの値をセル値とする数値ラスターが作成されます。

通常、点群データにおける点の配置は不規則でなので、セル内には 0 以上の任意の数の点が存在することになります。NumericRasterizer は、同一のセルに複数の点が存在する場合、データ内に格納されている順番で最後の点の値をそのセルの値として採用し、ひとつも点が存在しないセルには Background Value (背景の値) パラメーターに設定した値を与えます。

ここで、Fill Background with Nodata パラメーターを Yes に設定した場合は背景の値が「データがない」ことを示す Nodata として取り扱われ、FME Data Inspector では透過表示となります。

セルサイズ 2m x 2m の変換結果では、ごま塩のように透過表示 (Nodata) のセルが散在しています。セルサイズを大きくすればセル内に点が存在する確率が高くなり、10m x 10m ではほぼ全域のセルに値が設定されましたが、それでも右下に透過表示の部分が見られます。これは、比較的広い範囲で地表 (classification = 2) の点が存在しない水面であると思われます。


RasterDEMGenerator による変換

FMEワークスペース例 2

[LAS] リーダー: LAS形式の点群データを読み込む。
RasterDEMGenerator: 点群をラスターに変換する。
Inspector: 変換結果を FME Data Inspector に出力する。

RasterDEMGenerator では、Output DEM X Cell Spacing, Output DEM Y Cell Spacing パラメーターによってセルサイズを指定できます。NumericRasterizer と同様、2m x 2m, 5m x 5m, 10m x 10m の3通りの設定で試しました。

変換結果: 左からセルサイズ 2m x 2m, 5m x 5m, 10m x 10m

RasterDEMGeneretor は、指定されたサイズのセルの中心点について、その周辺の点 (セル内にあるとは限らない) の z 値に基づいて値を求めるので、点が存在しないセルにも値が与えられます。

セル値の求め方には AUTO, PLANAR, CONSTANT の3通りがあり、Interpretation Method (補間方法) パラメーターによって選択できます。これらのオプションは、「基盤地図情報DEMによるラインのドレープ (3Dラインへの変換)」で取り上げた SurfaceModeller や SurfaceDraper トランスフォーマーによるドレープと共通です。
----------
AUTO: PLANAR による補間方法で値が求められるセルではその値、PLANAR では Nodata となるセルでは CONSTANT の方法で求められる値。
PLANAR: 点群の点を頂点とする TIN (不規則三角網) を仮定し、XY平面に投影したときにセル中心点を含む3D三角形の3頂点の座標に基づく補間値をそのセルの値とする。ただし、中心点が TIN の外側に位置するセルの値は Nodata とする。
CONSTANT: XY平面に投影したときにセル中心点に最も近い点の z 値をそのセルの値とする。
----------

PLANAR による補間方法は barycentric interpolation と呼ばれ、セル中心点を通り z 軸に平行な直線と TIN を構成する3D三角形が交わる点の z 値が求められます。

上記の変換結果は AUTO によって実行したものです。ラスター外縁部で濃淡の変化が少なく内側のセルの値と不連続な「枠」のように見える部分は、PALANAR による補間方法では値が求められなかったところですが、CONSTANT の適用結果としても少々疑問です。補間方法としては PLANAR を選択し、この「枠」のような部分は Nodata とした方が良いかも知れません。

開口部のある3Dモデルのフットプリント

3Dモデル (ソリッド、サーフェス) を平面に投影することによって形成される2D幾何図形をフットプリント (footprint) と呼び、FME では、SurfaceFootprintReplacer トランスフォーマーによってソリッドやサーフェスをフットプリントに変換することができます (FME 2015.1 以降)。

ただし、現在のバージョン (FME 2016.1) では、上下方向に貫通する開口部がある3Dモデルについて、このトランスフォーマーによる変換後のフットプリントには、その開口部に対応する「穴」 (hole) が設けられないという機能上の制約があります。

将来、開口部に対応する穴を作成するように機能強化されるかも知れませんが、それを待たなくても、いくつかのトランスフォーマーを組み合わせて使用することによって、穴のあるフットプリントを作成することができます。

ここでは、「Esri Shapefile 形式による3Dモデルの保存」で使用した IFC データセットからフィーチャーをひとつ抽出し、それをフットプリント (開口部に対応する穴を含む) に変換するワークスペース例を掲げます。

ソースデータセット (IFC): DC_Riverside_Bldg-LOD_100.ifc

3Dモデルのデータ構造については「Esri Shapefile 形式による3Dモデルの保存」を参照してください。

FME 2016.1.0.1 build 16494

FMEワークスペース例

FME 2016.1 Workbench では、Canvas 上のオブジェクト間の接続ライン (connection) を非表示にすることや、接続ライン上にジャンクション (junction, 分岐・合流点) を設けることができるようになりました。ふたつのジャンクションの間を非表示にした部分をトンネル (tunnel) と呼び、接続ラインの右クリックメニュー > Create Tunnel によって作成することができます。

(テスト用ソリッドフィーチャーの抽出)
[IFC] リーダー: IFC データセットを読み込む。
Tester: テスト用のフィーチャーをひとつ抽出する (GlobalId = 1O62ieKrbD9fMT_ZiaQ9IA)。

(ソリッドのフットプリントへの変換)
2DForcer: ジオメトリを2D化する。
Deaggregator (Mode: Flatten All Levels): 集約ジオメトリを個別のジオメトリに分解する。
GeometryFilter: 面 (ポリゴンまたはドーナツ) のみを抽出する。
Dissolver: 隣接・重複する面を融合して単一の面にする。

IFC データセットから読み込まれたフィーチャーのジオメトリはソリッド (Solid) と Null で構成されている集約ジオメトリ (Aggregate) であり、これに 2DForcer を適用すると、ソリッドの部分が2D化されます (Null の部分には変更が加えられません)。

ソリッドを 2DForcer で2D化することによって得られるものは、「ソリッドの個々の面をそれぞれXY平面に投影することによって形成される複数のフットプリント (2Dジオメトリ) で構成される集約ジオメトリ」です。したがって、2DForcer が出力するフィーチャーは、「複数の2Dジオメトリで構成される集約ジオメトリ」と Null で構成される階層構造を持った集約ジオメトリを持つことになります。

これに含まれる個別の2Dジオメトリは、元のソリッドの上下の面をXY平面に投影した面 (この例では「穴」のある面=ドーナツ) と、垂直の壁面を投影したラインです。ソリッド全体のフットプリントとしてライン (垂直の壁面のフットプリント) は不要なので、Deaggregator によって個別のジオメトリに分解したうえで、GeometryFilter で面のみを抽出しました。階層構造を持った集約ジオメトリは、Deaggregator の Mode パラメーターを "Flatten All Levels" に設定することにより、全ての階層の個別のジオメトリに分解できます。

元のソリッドの上面と下面のフットプリントは重複しているので、Dissolver によってそれらを融合して単一の面にしました。これでソリッドのフットプリント (開口部に対応する穴を持つ) の完成です。

入力フィーチャーが複数ある場合は、Dissolver の Group By パラメーターで各フィーチャーを識別するための属性を指定すれば、元のフィーチャーごとに融合が行われます。そのような属性がない場合は、あらかじめ Counter によって連番を格納する属性を元のフィーチャーに付加しておき、それを ID 属性として使用することができます。

左: 元のソリッド, 右: フットプリント (上記ワークスペース例による変換結果)
どちらも FME Data Inspector によって表示したものです。

実は、Dissolver は、入力フィーチャーが集約ジオメトリを持っている場合は自動的に個別のジオメトリに分解し、面 (ポリゴンまたはドーナツ) のみを対象として処理を行います。したがって、上記ワークスペース例から Deaggregator と GeometryFilter を省いても同じ結果が得られます。

これを実行すると、ソリッドの垂直の壁面を2D化したラインや Null の部分は Dissolver の <Rejected> ポートから出力されるとともに、変換ログにいくつかの警告 (Warning) メッセージが出力されます。それらが処理結果に影響を与えることはありませんが、ジオメトリの変換過程を分かり易く表現するという意味では、Deaggregator, GeometryFilter を省略しないでおいた方が良いかも知れません。

補足: ソリッドの上側と下側がそれぞれ単一の面であり、かつ、全ての壁面が垂直である場合は、2DForcer -> Deaggregator -> GeometryFilter (Area) によって、ひとつのソリッドにつき、同じ形状・位置のフットプリントが2つ (元のソリッドの上下の面を投影した図形) 得られるので、Dissolver によってそれらを融合するのではなく、Sampler や DuplicateFilter などによってそのうちのひとつを選択する方法もあり、その方が効率が良いと予想します。

ワークスペース例で Dissolver を使ったのは、元のソリッドの上側、下側が複数の面で構成されている場合や壁面が垂直でない場合にも対応できるようにするためです。

例: 傾いた柱状のソリッドとそのフットプリント

GeoJSON の変換 - 国土地理院「ベクトルタイル提供実験」地形分類

GeoJSON は JSON (JavaScript Object Notation) の仕様に基づく地理空間データフォーマットです。GML などの XML ベースのフォーマットよりも軽量 (同じ情報量に対するデータサイズが小さい) であるとともに、JavaScript で扱い易いため、ネットワーク経由で地理空間データを配信するために広く利用されています。OGC (Open Geospatial Consortium) の規格ではありませんが、ウェブアプリケーションで地理空間データを扱う場合の事実上の標準フォーマットと言って良いと思います。

FME では、XML と同様に Fragmentation (断片化) と Flattening (平坦化) という2つの手法によって、任意の JSON ドキュメントからメンバーや値を抽出することができます。断片化、平坦化については「XMLの読込 - 断片化と平坦化」も参照してください。

ここでは、国土地理院「ベクトルタイル提供実験」において GeoJSON 形式で公開されている「地形分類」データについて、一般のGISソフトウェア等で利用可能なデータ (ジオメトリと属性) に変換するワークスペース例を掲げます。

データは、国土地理院「ベクトルタイル提供実験」地形分類 (GitHub) に掲載されているサンプルデータの URL から直接読み込むこととします。

なお、「地理院タイルデータの取得 - 西之島付近噴火活動 正射画像の例」と「ベクタージオメトリの色 - 国土地理院「地形分類」」で地理院タイルデータを取得・変換する JpGsiTileFetcher カスタムトランスフォーマーを紹介しましたが、その実装においてベクトルタイルを取得・変換する部分は、基本的には今回のワークスペース例と同じです。

FME 2016.1.0.1 build 16494


FMEワークスペース例

Creator: 処理開始用のフィーチャーを1個作成する。
HTTPCaller: HTTPプロトコル (Get メソッド) によってサンプル GeoJSON ドキュメントを取得する。
JSONFragmenter: GeoJSON ドキュメントを断片化して個別のフィーチャーに分解する。
JSONFlattener: 各フィーチャーを平坦化して "geometry" と "properties" メンバーを抽出する。
GeometryReplacer: "geometry" メンバーに基づいてジオメトリ (ポリゴン) を作成する。
JSONFlattener_2: "properties" メンバーを平坦化して下位のメンバー (属性) "code" を抽出する。

2017-08-25: FME 2017.1.1 では JSONFragmenter に Fragment as Format パラメーターが追加されました。そのパラメーターに "GEOJSON" を選択することにより、GeoJSON ドキュメントに基づいてジオメトリと属性を持つフィーチャーを作成することができます。
FME 2017.1.1 のワークスペース例は、本文末尾に掲載します。

HTTPCaller は、Request URL パラメーターに指定した URL に対して HTTP Method パラメーターで選択したメソッドのリクエストを発行し、レスポンスとして得られたデータを属性 (_response_body) に格納します。この場合は、サンプルデータとして提供されている「地形分類」 GeoJSON ドキュメントが _response_body に格納されます。

ひとつの GeoJSON ドキュメントには "features" という名前の配列が含まれており、その要素としてひとつ以上のフィーチャーオブジェクトが格納されています。各フィーチャーオブジェクトには "goemetry" と "properties" というメンバーオブジェクトが含まれており、それらによってフィーチャーのジオメトリと属性セットが表されます。

国土地理院「ベクトルタイル提供実験」地形分類データ (GeoJSON) の構造
分かり易くするために改行、インデントをつけました。

{
    "type": "FeatureCollection",
    "features": [
        {
            "geometry": {
                "type": "Polygon",
                "coordinates": [ <座標値の配列。詳細略> ]
            },
            "type": "Feature",
            "properties": {
                "code": 10714
            }
        },
        {
            "geometry": {
                "type": "Polygon",
                "coordinates": [ <座標値の配列。詳細略> ]
            },
            "type": "Feature",
            "properties": {
                "code": 11000
            }
        },
        <中略>
    ]
}

断片化 JSONFragmenter

JSONFragmenter の JSON Query パラメーターに次のクエリを設定することにより、このJSON オブジェクトを "features" 配列に格納されている個別のオブジェクト (フィーチャー) に分解することができます。
----------
JSON Query: json["features"][*]
----------
ここで、'json' は JSON ドキュメントのルートを示し、それに続く [ ] で下位の各階層のどのメンバー/要素を処理の対象とするかを指定することができます。[*] (アスタリスク) はその階層における任意のメンバー/要素を示すので、このクエリによって、"features" 配列の全ての要素 (フィーチャー) を JSON オブジェクトの断片として取り出すことができます。本記事作成時点で公開されていたサンプルデータには、58個のフィーチャーが格納されていました。

"features" 配列の要素 (JSON オブジェクトの断片)

{
    "geometry": {
        "type": "Polygon",
        "coordinates": [ <座標値の配列。詳細略> ]
    },
    "type": "Feature",
    "properties": {
        "code": 10714
    }
}

平坦化 JSONFlattener (1)

JSONFlattener によってオブジェクト (フィーチャー) を平坦化して "geometry" オブジェクトと "properties" オブジェクトを属性として抽出するとともに、Attributes to Expose パラメーターでそれらの名前を指定することにより、ワークベンチのインターフェース上に現しました。

Recursively Flatten Objects/Arrays を Yes に設定すると、下位にある全てのオブジェクト/配列を再帰的に平坦化して属性として抽出できますが、ここでは No に設定し、再帰的な平坦化は行いません。

GeoJSON フォーマットで記述されている "geometry" オブジェクトは次の GeometryReplacer によって直接ジオメトリデータに変換することができるので、その下位にある "coordinates" 配列のひとつひとつの要素 (ジオメトリを構成する個々の頂点の座標値) を個別の属性として抽出する必要はないためです。"properties" オブジェクトの平坦化 (個別属性の抽出) は、2番目の JSONFlattener で行います。

ジオメトリの作成 GeometryReplacer と平坦化 JSONFlattener (2)

GeometryReplacer (Geometry Encoding: GeoJSON) によって "geometry" オブジェクトをジオメトリデータ (この例ではポリゴン) に変換し、2番目の JSONFlattener によって "properties" オブジェクトを平坦化して下位のメンバー = 個別属性 (この例では "code" のみ) を抽出しました。GeometryReplacer と 2番目の JSONFlattener は順不同です。


変換結果: FME Data Inspector による表示

座標系について

上記ワークスペース例で作成されるフィーチャーには、座標系が定義されません。

GeoJSON の仕様（2017-08-25 注参照）では、"crs" という名前のメンバーオブジェクトによって空間参照系 (coordinate reference system) を定義することもできるのですが、 サンプルデータには "crs" は含まれていませんでした。"crs" が省略されている場合は、デフォルトの地理空間参照系 = WGS84 測地系 (十進緯度経度) であると解釈すべきとされています。

したがって、上記ワークスペース例に CoordinateSystemSetter を追加して WGS84 緯度経度の座標系 (LL84) を設定すれば、地理空間データとして完全になります。

2017-08-25 注: GeoJSONの最新の仕様では、"crs"メンバーについての規定は削除されており、ジオメトリの座標はWGS84測地系の十進緯度経度で記述することに統一されています。

GEOJSON リーダーによる GeoJSON ドキュメントの読込

FME は GeoJSON フォーマットに対応するリーダー/ライターを備えているので、JSONFragmenter, JSONFlattener などのトランスフォーマーを使わずに GEOJSON リーダーだけでジオメトリデータの作成、属性の抽出を行うこともできます。その場合は、座標系 ("crs" メンバーが省略されている場合は LL84) も自動的に設定されます。

FMEワークスペース例
上: GEOJSON リーダー, 下: FeatureReader トランスフォーマー
FeatureReader では、Format として GEOJSON を指定することにより、GEOJSON リーダーと同じモジュールによるデータ読込処理が行われます。

どちらも Dataset として URL を指定することにより、ウェブサーバーにアクセスして GeoJSON ドキュメントを読み込むこともできます。

この場合、見かけ上はウェブサーバー上のデータを直接読んでいるように見えますが、実際は、一旦 GeoJSON ファイルをダウンロードして作業用の一時フォルダに保存し、それを読み込むという仕組みであるため、ディスクにアクセスするためのオーバーヘッドがかかります。

HTTPCaller によってレスポンスとして得られた GeoJSON ドキュメントを属性に格納し (ファイルには保存せず)、いくつかのトランスフォーマーで変換した方が効率が良い可能性があるので、性能が重視されるプロジェクトでは、実行環境において計測したうえでどの方法を採用するか決定すべきでしょう。

---

2017-08-25: FME 2017.1.1 では、JSONFragmenter (Fragment as Format: GEOJSON) だけで、GeoJSON ドキュメントに基づいてジオメトリと属性を持つフィーチャーを作成することができます。

FMEワークスペース例 (FME 2017.1.1 build 17539)

Tokyo Datum (旧日本測地系) から JGD2000 への座標変換

日本測地系が Tokyo Datum (旧日本測地系) から JGD2000 (Japan Geodetic Datum 2000, 日本測地系2000) に切り替えられてから十数年、さらに JGD2011 に改訂されてからも数年が経ち、さすがに最近では Tokyo Datum によって作成された地理データは少なくなりましたが、それでも時折見かけることがあります。

Tokyo Datum のデータを JGD2000 に座標変換するには、測地系を定義するパラメーターを用いた幾何学的な計算による方法と、国土地理院が公開しているグリッドシフトファイル "TKY2JGD.par" で定義されている「地域ごとの変換パラメーター」を用いた内挿補間による方法の二通りがあります。

FME ワークスペースのデータフローの途中で座標変換を行うには Reprojector トランスフォーマーを使うのが一般的かつ簡便ですが、これはグリッドシフトファイルによる変換はサポートしません。しかし、FME は "TKY2JGD.par" を標準でバンドルしており、CsmapReprojector トランスフォーマーを使うことにより、それを用いた Tokyo Datum から JGD2000 への座標変換も行うことができます。

ここでは、これらのトランスフォーマーによって Tokyo Datum から JGD2000 に座標変換を行うワークスペース例を掲げます。サンプルとしたデータは、ウェブ上で公開されているデータのうち Tokyo Datum で作成された数少ない例のひとつ、国土交通省 50万分の1土地分類基本調査「地形分類図」 (Esri Shapefile 形式) です。

FME 2016.1.0.1 buid 16494

FMEワークスペース例

[ESRISHAPE] リーダー: 地形分類図 (Esri Shape) を読み込む。
Reprojector: 座標変換 (Tokyo Datum -> JGD2000)
CsmapReprojector: 座標変換 (Tokyo Datum -> JGD2000, TKY2JGD.par 使用)
Inspector: 変換結果を FME Data Inspector に出力。

Reprojector パラメーター設定画面

CsmapReprojector パラメーター設定画面

Reprojector、CsmapReprojector ともに、Source Coordinate System パラメーターで変換前の座標系名、Destination Coordinate System パラメーターで変換先の座標系名を指定 (選択) します。ここで指定する座標系名は、FME が定義する座標系の識別子で、"Tokyo" は「Tokyo Datum 緯度経度」、"LL-JGD2K" は「JGD2000 緯度経度」を指します。

# FME がサポートする座標系は、Coordinate System Gallery (Workbench メニュー: Tools > Browse Coordinate System) に一覧表示されます。

ソースデータの内容として座標系が正しく定義されていれば (Esri Shapefile 形式の場合は *.prj ファイルによる)、ソース側の座標系名の指定は省略できます。

CsmapReprojector では、これらに加えて、Transformation パラメーターで変換方法を選択することができます。
"TKY2JGD.par" を用いて Tokyo Datum から JGD2000 への座標変換を行う場合は、パラメーター設定フィールド右側の [...] ボタンをクリックして Geographic Transformation Gallery を開き、"Tokyo-Grid_to_JGD2000" を選択します。

Geographic Transformation Gallery

次の図は、上記ワークスペースの実行結果 (FME Data Inspector への出力) のスクリーンショットです。View 画面上での目視では判別できないほどごくわずかではありますが、Feature Information ウィンドウに表示される頂点の座標値には差があることが確認できます。

Reprojector による変換結果

CsmapReprojector による変換結果 (TKY2JGD.par 使用)

"TKY2JGD.par" で定義されている「地域ごとの変換パラメーター」は、全国の陸部をカバーする範囲について、南北30秒、東西45秒間隔のグリッドポイント、つまり3次メッシュ区画の隅における新旧測地系間の緯度、経度の差分であり、ある地点の座標変換量は、その地点を含む区画四隅のパラメーターの内挿補間 (バイリニア) によって求められます。

この「地域ごとの変換パラメーター」の値は、旧測地系の測量成果における誤差や地殻変動等によって生じた「測地網の歪み」も考慮して求められたものなので、CsmapReprojector (TKY2JGD.par 使用) の方が変換精度が高いとは言えますが、Reprojector による変換結果が「間違い」というわけではありません。実用上、これらのトランスフォーマーによる変換結果の差程度のずれは問題にならないことも多いと思います。

なお、TKY2JGD.par がカバーしていない地点の座標の CsmapReprojector (TKY2JGD.par 使用) による変換結果は、Reprojector による変換結果と同じになります。

ラスター画素値演算の例 - 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)

基盤地図情報DEMに基づく陰影起伏ラスターの作成

「基盤地図情報DEMに基づく傾斜方向ラスターの作成」で触れたように、DEMラスターは RasterHillShader トランスフォーマーによって陰影起伏 (Hillshade) ラスターに変換することができます。このトランスフォーマーが作成する陰影起伏ラスターは白黒の濃淡によって陰影を表現したグレースケールのラスターですが、DEMラスターをカラー画像ラスターに変換したうえで、グレースケールの陰影起伏ラスターから得られる濃淡の比率に基づいて各画素のRGB値を調整することにより、カラーの陰影起伏ラスターを作成することもできます。

ここでは、基盤地図情報DEMデータに基づいて陰影起伏ラスター (グレースケール、カラー) を作成するワークスペース例について、次の2段階に分けて説明します。

1. 基盤地図情報DEMラスターから標高区分カラー画像ラスターへの変換
2. 標高区分カラー画像ラスターから陰影起伏ラスターへの変換

サンプルとしたデータは、基盤地図情報10mメッシュDEM (三宅島の区域を含む2次メッシュ4区画) です。

FME 2016.0.1.2 build 16178


1. 基盤地図情報DEMラスターから標高区分カラー画像ラスターへの変換

FMEワークスペース例 1

[JP_FGD_DEM] リーダー (カスタムフォーマット): 基盤地図情報DEMデータ (XML) をラスターとして読み込む。
RasterMosaicker: 2次メッシュ4区画分のラスターを結合してひとつのラスターを作成する。
Reprojector: 座標系を変換する。
RasterExpressionEvaluator: 標高区分値 (整数) をセル値とするラスターに変換する。
RasterBandMinMaxExtractor: セル値の最小、最大値を求める。
RGBGradationCalculator (FME Hub): 各標高区分値に対応するRGBカラー値を求め、パレットを定義する。
StringConcatenator: パレット定義に Nodata 値用のエントリーを追加する。
RasterPaletteAdder: パレットをラスターに適用する。
RasterPaletteResolover: パレットに基づいて 4バンド (RGBA) ラスターに変換する。
[PNGRASTER] ライター: PNG形式でファイルに出力する。

DEMラスターのセル値は標高値 (実数) ですが、それをある基準によって区分して区分ごとの番号 (整数値) を割り当て、それらの番号をキー値として対応するRGBカラー値を定義したパレットをラスターに与えることによってカラー画像に変換するというアプローチです。

標高値の区分、および、各区分と色の対応づけの具体的な方法は、最終的に得るべき画像の仕様に応じて検討する必要がありますが、この例では、単純に標高値を10m間隔で区分し、対象とする区域 (三宅島周辺) の標高値の範囲で、次の図のように段階的に変化する色を与えることとします (左: 低 -> 右: 高)。

RasterExpressionEvaluator では、次の式によって、DEMラスターのセル値 (標高値) に基づいて10m間隔の区分値 (整数) を求め、それをセル値とするラスターに変換しました。

@if(@isnodata(A[0]), 0, @if(A[0]<10, 1, @uint16(A[0]/10.0)+1))

この式は、標高値が Nodata (海域などでデータがない) の場合は 0, 10m 未満の場合は 1, それ以外の場合は標高値を 10 で除した値の小数部を丸めてから 1 を加えた値を求めています。これによって、Nodata のセルは 0、それ以外のセルは 1 以上の標高区分値 (値が大きいほど標高が高いことを示す連番) のラスターに変換することができます。

RasterMinMaxExtractor は、入力ラスターのバンドごとにセル値の最小値と最大値を求め、リスト属性 _band{}.min, _band{}.max に格納します。Nodata 以外のセル値は 1 から始まる連番の標高区分値としたので、その最大値が、定義すべきパレットのエントリー数 (色の数) ということになります。Nodata 用のエントリーは後で追加します。

RGBGradationCalculator は FME Hub で公開されているカスタムトランスフォーマーで、指定した複数の色 (RGB値) の間を補間して段階的に変化するRGB値を求めるとともに、ラスター用のパレット定義 (文字列) を格納した属性を作成することができます。主なパラメーターの設定内容は次のとおりです。

Seed Colors: #0000ff; #00ffff; #00ff00; #ffff00; #ff0000 (色の範囲の先頭、等分割した中間、末尾のRGB値)
Format of Seed Colors: HTML Color (上記のRGB値の記述形式=HTMLフォーマット)
Number of Destination Color: _band{0}.max (作成する色の数=標高区分値の最大値)
Raster Palette Attribute Name: _palette (パレット定義を格納する属性名)
Raster Palette Interpretation: RGBA32 (パレットの種類=8ビットRGBA4バンド）
Raster Palette Alpha Value: 255 (Alphaバンドの既定値=不透過)
Raster Palette First Key Value: 1 (最初のエントリーのキー値)
Raster Palette Omit Interpretation: Yes (パレット定義先頭行のパレットの種類の記述を省略)

この設定によって、標高区分値 (キー値) 1 以上の部分のパレット定義を格納した属性 _palette が作成できます。さらに、次の StringConcatenator によって、先頭行のパレットの種類の記述 ('RGBA32') と Nodata 用のエントリー (キー値 = 0, RGBA値 = 0,0,0,0) を先頭部分に追加してパレット定義を完成させました。

パレット定義の最初の数行は、次のようになります。パレット定義の書式については、RasterPaletteExtractor トランスフォーマーのヘルプを参照してください。

RGBA32
0 0,0,0,0
1 0,0,255,255
2 0,13,255,255
3 0,26,255,255
4 0,40,255,255
5 0,53,255,255
(以下略)

RasterPaletteAdder によってこのパレット定義をラスターに与えるとともに、陰影起伏ラスターに変換するための準備として、RasterPaletteResolver によってパレットのRGBA値を各バンドの値とするラスター (8ビットRGBA4バンド) に変換しました。ここまでで、次の図のようなカラー画像ラスターが得られます。

結果 1: 標高区分カラー画像ラスター - RasterPaletteResolver からの出力
基盤地図情報10mメッシュDEM (三宅島を含む2次メッシュ4区画) に基づき作成

2. 標高区分カラー画像ラスターから陰影起伏ラスターへの変換

上記のワークスペース例にいくつかのトランスフォーマーを追加することにより、陰影起伏ラスターが作成できます。

FMEワークスペース例 2 (追加部分)

RasterHillshader: DEMラスターを陰影起伏ラスター (グレースケール) に変換する。
RasterBandNodateRemover: 標高区分カラー画像から Nodata 定義を削除する。
RasterExpressionEvaluator: 標高区分カラー画像のRGB値を調整し、陰影起伏ラスター (カラー) に変換する。

RasterExpressionEvaluator では、Mode パラメーターで "Two Rasters" を選択することにより、2つのラスター A, B に基づくセル値の演算を行うこともできます。ここでは、RasterHillshader で作成したグレースケールの陰影起伏ラスターを A、標高区分カラー画像ラスターを B としました。

RasterExpressionEvaluator によって2つのラスターに基づくセル値の演算を行う場合、A, B の範囲、解像度は厳密に同じであるとともに、Nodata 値の定義の有無、定義がある場合はその内容も同じである必要があります。この例では、RasterHillshader で作成した陰影起伏ラスター (A)、カラー画像ラスター (B) ともに同一のDEMラスターから変換した後、範囲、解像度を変更する処理はしていないので、それらが同じであることは保証されていますが、A には Nodata 値の定義がないのに対して、B の方には Nodata 値の定義があります。そのため、RasterBandNodataRemover によって B から Nodata 値の定義を削除しました。

RasterHillshader が作成するグレースケールの陰影起伏ラスターは、デフォルトでは 2 バンド (0-255 の範囲で濃淡を表すバンド、および、Nodata だったセルを 0: 透明、それ以外を 255: 不透明とする Alpha バンド) を持っています。したがって、グレースケール陰影起伏ラスター (A) の最初のバンドのセル値 A[0] を 255 で除した値 (濃淡の比率) を標高区分カラー画像ラスター (B) の R, G, B バンドのセル値にそれぞれ乗じることにより、陰影のついた R, G, B 値が求められ、また、A ラスターの第2バンドのセル値 A[1] を Alpha 値とすることで、Nodata だったセルを透明にすることができます。

以上の考え方に基づくパラメーター設定例は、次のとおりです。

結果 2: 陰影起伏ラスター (グレースケール) - RasterHillshader からの出力
基盤地図情報10mメッシュDEM (三宅島を含む2次メッシュ4区画) に基づき作成

結果 3: 陰影起伏ラスター (カラー) - RasterExpressionEvaluator_2 からの出力
基盤地図情報10mメッシュDEM (三宅島を含む2次メッシュ4区画) に基づき作成