ポイントクラスによる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)