Zenn
🗺️

QGISでProject PLATEAUのCityGMLの緯度経度が入れ替わる問題を修正した件

2021/09/20に公開
QGIS
#
gdal
PLATEAU
#
foss4g
tech

Project PLATEAU

Project PLATEAUは国交省が主導する3D都市モデルのオープンデータ化プロジェクトで、CityGMLと呼ばれる形式のデータを2021年9月現在、56都市分公開しています。

https://www.mlit.go.jp/plateau/

緯度と経度が入れ替わる問題

QGIS 3.20.3（2021年9月現在最新版; GDAL 3.0.4使用）にProject PLATEAUのCityGMLをレイヤー追加すると、緯度と経度が取り違えられて、北緯90度よりもさらに北のとんでもないところに地物が追加されてしまいます。
QGIS3.20.3-EPSG6697
データ：plateau-tokyo23ku-citygml-2020/533946_2.zip内 bldg/53394654_bldg_6697_op2.gml; スカイツリー周辺

そのため、EPSG:3857 Webメルカトル図法のような高緯度（85.06以北）を切り捨てる投影法では表示されません。

この問題はデータの緯度と経度を入れ替えてあげることで、表示上で対応する回避策が紹介されていました。

https://qiita.com/Yfuruchin/items/aa3a847db26c7f378d6e

本問題を調査したところ、GDALのジオメトリ読み出し処理に問題があると判明したので、修正を行いpull-requestを申請しました。本修正はGDALに取り込まれ、 3.3.3, および3.4.0で公開されました。

https://github.com/OSGeo/gdal/pull/4505

本修正を含むGDALを使用したQGISでは座標が正しく読み込まれることが確認できます。

QGIS-master-GDAL3.3.2+-EPSG6697
ベースマップ：地理院タイル 標準地図、データ：同上

以下は技術的な詳細になります。

何が問題だったのか？

Project PLATEAUのCityGMLの座標情報はgml:posListで緯度 経度 高さの頂点情報を記録しています。
これはJGD2011の仕様に基づいたシリアライズ方法で、WGS84と同等の仕様です。

Project PLATEAUの3D都市モデル標準作業手順書 (PDF 2.9MB) には次のように定義されています。

B.3.1 gml:pos, gml:posList

gml:Envelope の属性 srsName 座標で指定した空間参照系に基づく座標値を記述する。「日本測地系 2011 における経緯度座標系と東京湾平均海面を基準とする標高の複合座標参照系」の定義に従い、座標値は緯度、経度、標高の順列とし、それぞれを半角スペースで区切る。

この緯度(Y) 経度(X) 高さ(Z)の順列を経度(X) 緯度(Y) 高さ(Z)と読み出しているために、緯度が経度に、経度が緯度に入れ替わってしまいました。

原因1. 複合型CRSが考慮されていない

主な原因はGDALのジオメトリ情報の読み出し処理で、複合型の座標参照系（CRS）が正しく解釈されないことにありました。

GDALはQGISが使用している地理空間データ処理ライブラリで、このライブラリが多岐にわたるデータフォーマットの解釈を担っています。
GDALはデータを読み出す際、緯度経度の順にシリアライズされることが定義された座標系（WGS84やJGD2011）の場合に、読み出した座標のXYをYXに入れ替える処理を持っています。

国土地理院が公開している基盤地図情報のGMLファイルはJGD2011でシリアライズされていますが、QGISでも正しく読むことができます。GDALがJGD2011の緯度経度シリアライズに対応していない、というわけではなさそうです。
それでは何故、基盤地図情報のGMLは読めて、Project PLATEAUのCityGMLは読めないのでしょうか？

そもそも、Project PLATEAUのCityGMLの座標参照系はGMLの冒頭にsrsNameとして定義されています。

<gml:boundedBy>
	<gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/6697" srsDimension="3">
		<gml:lowerCorner>35.70819416724163 139.79979831321313 -0.974</gml:lowerCorner>
		<gml:upperCorner>35.71682557957052 139.8126483718091 100.486</gml:upperCorner>
	</gml:Envelope>
</gml:boundedBy>

ここで定義されているEPSG:6697というのは「JGD2011 + JGD2011 (vertical) height」という座標参照系で、EPSG:6668（JGD2011）とEPSG:6695（JGD2011 (vertical) height）からなる座標系であると定義されています。

GDALでEPSG:6697のWell-known textを出力すると次のようになります。

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> ref.ImportFromEPSG(6697)
>>> print(ref.ExportToPrettyWkt())
COMPD_CS["JGD2011 + JGD2011 (vertical) height",
    GEOGCS["JGD2011",
        DATUM["Japanese_Geodetic_Datum_2011",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","1128"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AXIS["Latitude",NORTH],
        AXIS["Longitude",EAST],
        AUTHORITY["EPSG","6668"]],
    VERT_CS["JGD2011 (vertical) height",
        VERT_DATUM["Japanese Geodetic Datum 2011 (vertical)",2005,
            AUTHORITY["EPSG","1131"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Gravity-related height",UP],
        AUTHORITY["EPSG","6695"]],
    AUTHORITY["EPSG","6697"]]

EPSG:6668はJGD2011としてお馴染みですが、これをGEOGCS（PROJ内部ではhorizontal CRS - 水平方向CRS）として継承（or 再利用）し、新たにVERT_CS（PROJ内部ではvertical CRS - 垂直方向CRS）として高さ方向の軸（EPSG:6695）を追加するという、既存のCRSを拡張する作りになっているようです。

この形式のCRSをGDALが使用する座標参照系ライブラリPROJの内部ではCompoundCRSと呼んでいるようです。（先述のProject PLATEAUの座標系の定義には「複合座標参照系」とあるので、本記事では「複合」で統一しました。）

さて、GDALはWGS84やJGD2011のような緯度経度形式のシリアライズフォーマットを次のような処理で判定し、ret == trueであれば、読み出したジオメトリのXYを交換する仕組みになっていました。

auto cs = proj_crs_get_coordinate_system(d->getPROJContext(),
                                            d->m_pj_crs);
if ( cs )
{
    const char* pszDirection = nullptr;
    if( proj_cs_get_axis_info(
        d->getPROJContext(), cs, 0, nullptr, nullptr, &pszDirection,
        nullptr, nullptr, nullptr, nullptr) )
    {
        if( EQUAL(pszDirection, "north") )
        {
            ret = true;
        }
    }

    proj_destroy(cs);
}

proj_crs_get_coordinate_systemはPROJの関数で、この関数にEPSG:6697のCRSオブジェクトが渡されると、csがnullptrになります。

PROJの本関数を調べてみると次のように実装されています。

PJ *proj_crs_get_coordinate_system(PJ_CONTEXT *ctx, const PJ *crs) {
    SANITIZE_CTX(ctx);
    if (!crs) {
        proj_context_errno_set(ctx, PROJ_ERR_OTHER_API_MISUSE);
        proj_log_error(ctx, __FUNCTION__, "missing required input");
        return nullptr;
    }
    auto l_crs = dynamic_cast<const SingleCRS *>(crs->iso_obj.get());
    if (!l_crs) {
        proj_log_error(ctx, __FUNCTION__, "Object is not a SingleCRS");
        return nullptr;
    }
    return pj_obj_create(ctx, l_crs->coordinateSystem());
}

注目すべきはCRSオブジェクトをSingleCRSオブジェクトにダイナミックキャストするdynamic_cast<const SingleCRS *>(...);の部分で、EPSG:6697のオブジェクトが渡された場合、この部分でnullptrが返されます。

そもそもEPSG:6697のような複合型のCRSはSingleCRSではなく、CompoundCRSという別のクラスとして作成されているため、複合型のCRSをproj_crs_get_coordinate_systemに渡すのは不適当ということになります。

複合型のCRSであることはproj_get_type関数の返り値がPJ_TYPE_COMPOUND_CRSになることで判定ができます。
そこで、入力が複合型のCRSの場合は水平方向CRS（今回の場合はEPSG:6668、つまり純粋なJGD2011）を取り出してから軸判定を、そうでないCRSの場合は従来どおりの判定を、という修正を行いました。

原因2. デフォルトCRSが反映されない

複合型のCRSへの対応を行うだけでは緯度経度の入れ替えが行われませんでした。

基盤地図情報のGMLの場合、CRSは次のように定義されています。

<?xml version="1.0" encoding="utf-8"?>
<Dataset xsi:schemaLocation="http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema FGD_GMLSchema.xsd"
	xmlns:gml="http://www.opengis.net/gml/3.2"
	xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
	xmlns:xlink="http://www.w3.org/1999/xlink"
	xmlns="http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema"
	gml:id="Dataset1">
<BldA gml:id="K17_5031012840_2">
<area>
<gml:Surface gml:id="K17_5031012840_2-g" srsName="fguuid:jgd2011.bl">
<gml:patches>
<gml:PolygonPatch>
<gml:exterior>
<gml:Ring>
<gml:curveMember>
<gml:Curve gml:id="K17_5031012840_2-3">
<gml:segments>
<gml:LineStringSegment>
<gml:posList>...</gml:posList>
</gml:LineStringSegment>
</gml:segments>
</gml:Curve>
</gml:curveMember>
</gml:Ring>
</gml:exterior>
</gml:PolygonPatch>
</gml:patches>
</gml:Surface>
</area>
</BldA>
...
</Dataset>

注目すべきはsrsName="fguuid:jgd2011.bl"の位置で、この属性値はファイル内の全てのgml:Surface要素に定義されています。
これをツリー構造にまとめると次のようになります。

  • Dataset
    • BldA
      • area
        • gml:Surface [srsName="fguuid:jgd2011.bl"]
    • BldA
      • area
        • gml:Surface [srsName="fguuid:jgd2011.bl"]
    • ...

（本筋からは外れますが、このfguuid:jgd2011.blはGDALにハードコーディングされており、あまり美しくない作りになっています…。これ以外に方法はないと思われるので仕方がないですが…。）

一方で、Project PLATEAUのCityGMLのsrsNameは次のように定義されており、基盤地図情報のGMLとは異なります。

<core:CityModel xmlns:grp="http://www.opengis.net/citygml/cityobjectgroup/2.0" xmlns:core="http://www.opengis.net/citygml/2.0" xmlns:pbase="http://www.opengis.net/citygml/profiles/base/2.0" xmlns:smil20lang="http://www.w3.org/2001/SMIL20/Language" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:smil20="http://www.w3.org/2001/SMIL20/" xmlns:bldg="http://www.opengis.net/citygml/building/2.0" xmlns:xAL="urn:oasis:names:tc:ciq:xsdschema:xAL:2.0" xmlns:uro="http://www.kantei.go.jp/jp/singi/tiiki/toshisaisei/itoshisaisei/iur/uro/1.4" xmlns:luse="http://www.opengis.net/citygml/landuse/2.0" xmlns:app="http://www.opengis.net/citygml/appearance/2.0" xmlns:gen="http://www.opengis.net/citygml/generics/2.0" xmlns:dem="http://www.opengis.net/citygml/relief/2.0" xmlns:tex="http://www.opengis.net/citygml/texturedsurface/2.0" xmlns:tun="http://www.opengis.net/citygml/tunnel/2.0" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:sch="http://www.ascc.net/xml/schematron" xmlns:veg="http://www.opengis.net/citygml/vegetation/2.0" xmlns:frn="http://www.opengis.net/citygml/cityfurniture/2.0" xmlns:gml="http://www.opengis.net/gml" xmlns:tran="http://www.opengis.net/citygml/transportation/2.0" xmlns:wtr="http://www.opengis.net/citygml/waterbody/2.0" xmlns:brid="http://www.opengis.net/citygml/bridge/2.0" xsi:schemaLocation="http://www.kantei.go.jp/jp/singi/tiiki/toshisaisei/itoshisaisei/iur/uro/1.4 http://www.kantei.go.jp/jp/singi/tiiki/toshisaisei/itoshisaisei/iur/schemas/uro/1.4/urbanObject.xsd http://www.opengis.net/citygml/2.0 http://schemas.opengis.net/citygml/2.0/cityGMLBase.xsd http://www.opengis.net/citygml/landuse/2.0 http://schemas.opengis.net/citygml/landuse/2.0/landUse.xsd http://www.opengis.net/citygml/building/2.0 http://schemas.opengis.net/citygml/building/2.0/building.xsd http://www.opengis.net/citygml/transportation/2.0 http://schemas.opengis.net/citygml/transportation/2.0/transportation.xsd http://www.opengis.net/citygml/generics/2.0 http://schemas.opengis.net/citygml/generics/2.0/generics.xsd http://www.opengis.net/citygml/cityobjectgroup/2.0 http://schemas.opengis.net/citygml/cityobjectgroup/2.0/cityObjectGroup.xsd http://www.opengis.net/gml http://schemas.opengis.net/gml/3.1.1/base/gml.xsd http://www.opengis.net/citygml/appearance/2.0 http://schemas.opengis.net/citygml/appearance/2.0/appearance.xsd">
	<gml:boundedBy>
		<gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/6697" srsDimension="3">
			<gml:lowerCorner>35.70819416724163 139.79979831321313 -0.974</gml:lowerCorner>
			<gml:upperCorner>35.71682557957052 139.8126483718091 100.486</gml:upperCorner>
		</gml:Envelope>
	</gml:boundedBy>
	<core:cityObjectMember>
		...
	</core:cityObjectMember>
	...
</core:CityModel>

同様にツリー構造にすると次のようになり、各要素にsrsNameがかかっていない、グローバルな定義であることが分かります。

  • core:CityModel
    • gml:boundedBy
      • gml:Envelope [srsName="http://www.opengis.net/def/crs/EPSG/0/6697" srsDimension="3"]
    • core:cityObjectMember
    • core:cityObjectMember
    • ...

（このwww.opengis.netのURIは単体でも動作し、XMLの内容を解釈して座標系の関係を理解することができて楽しいです。）

GDALの内部でgml:boundedBy/gml:Envelopeが検知されると、IGMLReader::SetGlobalSRSNameが呼ばれ、グローバルなSRSNameとして登録される実装になっていましたが、設定されたグローバルなSRSNameを参照する直前にIGMLReader::CanUseGlobalSRSNameが呼ばれ、この関数が常にfalseを返してしまうために、グローバルなSRSNameが使われない、という状態になっていました。

そこで、SetGlobalSRSNameでグローバルなSRSNameを設定したタイミングで、CanUseGlobalSRSNameがtrueを返すよう設定し、ジオメトリを読み出すタイミングでグローバルな座標系定義を参照できるようにしました。

まとめ

以上の修正がGDALのmasterブランチおよびGDAL 3.3ブランチにマージされています。

QGISの地理空間データと座標系の扱いの関係はQGIS -> GDAL -> PROJという多段構成になっており、問題を追跡するのに横断的な知識が必要なようです。
QGISやGDALのメンテナの方々はそれぞれのプロジェクトに横断的にコントリビュートされているようですね。GDALへのpull-requestは非常に素早く対応いただきました。感謝します。

今回はGDALライブラリをPythonインタプリタ上で呼び出し、基盤地図情報のGMLとProject PLATEAUのCityGMLをそれぞれ開いてジオメトリ読み込み状態の比較を行うことでGDALの問題として切り分け、原因箇所の特定を行いました。
GDALもPROJも素直なC++で書かれているため、ライブラリのコードを追いかけるのは楽な部類だと感じました。
lldb等で関数の呼び出し階層を取得すると、Reader等の状態把握が非常に楽に行えます。

pull-requestのレポートの大半をDeepL翻訳に頼り切っています。そのため少々不自然な英語があるかもしれません。

まさかGDALにコントリビュートできる日がくるとは思っていませんでしたが、コードは素直で手がつけやすく、メンテナの方々の応答も非常に素早いのでコントリビュートしやすいプロジェクトという印象を抱きました。
皆さんも何か問題に当たったときは気軽にコントリビュートしてみてはいかがでしょうか？

Discussion