以前に書いた、国土数値に公開されている GIS データを Mathematica 7 で利用する記事では、Windows 上で動作する変換ソフトを介したやり方を検討した。もちろん、Windows 限定のやり方だ。
→参照:Mathematica で日本の GIS データを (1)
僕はもともと ひどく Mac な人で、これまで会社では、不本意を抱えつつ Windows を使っていた。が、前回の記事のあと、長年の念願かなって、僕の会社のメインマシンを Mac に変更してもらうことができた(ひゃっほぅ!)。そこで今回は、Windows の変換ソフトを介さずに、国土数値の JPGIS データ(XML ファイル)を Mathematica で直接利用する方法を検討してみた。
JPGIS の XML 構成については、データと一緒に公開されているドキュメントもあるし、国土地理院のサイトには JPGIS の詳細な仕様書も公開されている。が、これらのドキュメント、僕がお気楽に試すにはちょっと敷居が高い。というか、こんなの読みこなすほどの根気ないデス。
というわけで、とりあえず Mathematica に読み込んでから、手探りでいろいろ試してみることにした。Mathematica の XML 機能で XML ファイルを読み込むと、XMLElement というシンボルで各種属性値をまとめたデータ構造が構築されるので、これを調べて行くわけだ。
今回は、僕の故郷の滋賀県のデータを使ってみる。世間では県の面積のほとんどが琵琶湖だと思われているらしい不憫(?)な県だ。
まず、国土数値の行政区域データから、平成 20 年の滋賀県のデータをダウンロードして展開。今回は、このデータを展開したフォルダに Mathematica ノートブックを保存してから作業を始めることにする。ノートブックを保存したら、まずは読み込み。
In[1]:=
SetDirectory[NotebookDirectory[]];
kuikixml = Import["N03-090320_25.xml"];NotebookDirectory は ver.6 で追加された関数で、ノートブックファイルを保存しているフォルダのファイルパスが分かる。SetDirectory と組み合わせれば、ノートブックと同じフォルダのファイルはファイル名だけで扱えるようになる。(ちなみに、ver.6 より前のバージョンで同じことをするには、SetDirectory["FileName" /. NotebookInformation[EvaluationNotebook[]] // First // ToFileName] など)
以前も少し書いたが、Mathematica の Import や Export で内部的に呼び出される外部プログラムには日本語を含むファイルパスに対応していない物がままあるので、Import / Export 関連の作業をするときには、フォルダ名やファイル名は無難な半角文字のみにしておこう。
前回、JPGIS を変換して読み込んだときには、データに行政区域の名称が含まれていた。ということは、XML からの読み込みでも、どこかに行政区域名があるはずだ。そこでまず、読み込んだ XML データから「彦根市」という文字列を含む階層を探してみる。
In[3]:=
Cases[kuikixml, _[___, "彦根市", ___], Infinity]
Out[3]=
{{彦根市}}おー、都合良く1個しかヒットしない。
この Cases での検索で、_[___, ..., ___] の階層を増やしていけば、「彦根市」を含むデータ階層をひとつずつ広げてみていくことができる。上で見つけた {"彦根市"} の部分を固定してさらに3階層広げると次のようになった。
In[4]:=
Cases[kuikixml, _[___, _[___, _[___, {"彦根市"}, ___], ___], ___],
Infinity]
Out[4]=
{XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, EC01}, {id -> ec00100013}, {XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, ARE}, {idref -> a00100013}, {}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, PRN}, {}, {滋賀県}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, CN2}, {}, {彦根市}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, AAC}, {codeSpace -> AdminAreaCd.xml}, {25202}]}]}id で示されているのが彦根市の ID だとすると、idref -> a00100013 は他のデータを参照するキーに違いない。ARE は Area かな? じゃあ次はこの a00100013 で同じことをやってみる。
In[5]:=
Cases[kuikixml, _[___, _[___, "a00100013", ___], ___], Infinity]
Out[5]=
{{id -> a00100013}, {idref -> a00100013}, {idref -> a00100013}}ここでは、Cases の結果として3つの場所が検出されている。idref が参照キーだとすると、本体データは {id -> a00100013} のものに違いないので、こちらを追って...、
In[6]:=
Cases[kuikixml, _[___, {"id" -> "a00100013"}, ___], Infinity]
Out[6]=
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Surface}, {id -> a00100013},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_OrientablePrimitive.orientation}, {}, {+}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_OrientablePrimitive.primitive}, {idref -> a00100013}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Surface.patch}, {},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Polygon}, {},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_SurfacePatch.interpolation}, {}, {planar}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Polygon.boundary}, {idref -> sfBd00100013},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_SurfaceBoundary}, {id -> sfBd00100013},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Complex.element}, {idref -> rn00100013}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_SurfaceBoundary.exterior}, {},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_Ring}, {id -> rn00100013},
{XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_OrientablePrimitive.orientation}, {}, {+}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_OrientablePrimitive.primitive}, {idref -> rn00100013}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_CompositeCurve.generator}, {idref -> _c00100048}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_CompositeCurve.generator}, {idref -> _c00100059}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_CompositeCurve.generator}, {idref -> c00100063}, {}],
XMLElement[{http://www.gsi.go.jp/GIS/jpgis/standardSchemas, GM_CompositeCurve.generator}, {idref -> c00100069}, {}],
.....<以下略>最後に並んでいる GM_CompositeCurve.generator というのは、名前からすると曲線データを示しているっぽい。ということは、ここの idref で示された ID の曲線データをまとめる形で行政区域の面輪郭が定義されているのに違いない。
...とかいうのを繰り返して、芋づる式にデータを辿って、地図を描くのに必要なデータ構造を調べた結果は次の通り。
- EC01 (区域名など)
- ARE
- idref ― 面輪郭(GM_Surface)の ID
- CN2
- 区域名文字列
- ARE
- GM_Surface (面輪郭データ)
- id ― 面輪郭の ID
- GM_SurfaceBoundary.exterior または GM_SurfaceBoundary.interior (繰り返しあり)
- GM_CompositeCurve.generator (繰り返しあり)
- idref ― 曲線 GM_Curve の ID。GM_Curve データを逆順に使うときは ID の前にアンダースコア(_)
- GM_CompositeCurve.generator (繰り返しあり)
- GM_Curve (曲線データ)
- id ― 曲線の ID
- GM_PointRef.point (DirectPosition.coordinate と合わせて繰り返しあり)
- idref ― ポイント GM_Point の ID
- DirectPosition.coordinate (GM_PointRef.pointと合わせて繰り返しあり)
- "緯度 経度" 形式の文字列
- GM_Point (ポイントデータ)
- id ― ポイントの ID
- DirectPosition.coordinate
- "緯度 経度" 形式の文字列
ここまで分かればあとは一本道。構造未知の XMLElement の深いところからデータを抽出するので、深さ無限を Infinity オプションで指定した Cases を多段にしてデータを抽出していく。上のまとめを逆順に辿る形。
In[7]:=
(* "緯度 経度" 文字列を {経度,緯度} リストに変換する関数 *)
convCoord[s_] := ToExpression[StringReplace[s,
y : NumberString ~~ Whitespace ~~ x : NumberString :>
"{" <> x <> "," <> y <> "}"]]
(* ポイントデータ *)
pts = Cases[kuikixml, XMLElement[{_, "GM_Point"},
{"id" -> id_}, {_[__, {_[__, {_["DirectPosition.coordinate",
_, {yx_}], _}]}]}] :> Rule[id, convCoord[yx]], Infinity];
(* 曲線データ *)
cvs = Cases[kuikixml,
XMLElement[{_, "GM_Curve"}, {"id" -> id_}, data__] :>
Rule[id, Cases[data, _[kind : "GM_PointRef.point" |
"DirectPosition.coordinate", {pid___}, {yx___}] :>
If[kind === "GM_PointRef.point", pid[[2]] /. pts,
convCoord[yx]], Infinity]], Infinity];
(* 面輪郭データ *)
sfs = Cases[kuikixml,
XMLElement[{_, "GM_Surface"}, {"id" -> id_}, data__] :>
Rule[id, Cases[data, _[{_, "GM_SurfaceBoundary.exterior" |
"GM_SurfaceBoundary.interior"}, _, {cvdata__}] :>
Polygon[Join @@ Cases[cvdata, _[{_,
"GM_CompositeCurve.generator"},
{"idref" -> cvid_}, _] :> If[StringMatchQ[cvid, "_" ~~ __],
Reverse[StringDrop[cvid, 1] /. cvs], cvid /. cvs],
Infinity]], Infinity]], Infinity];
(* 区域名データ *)
twns = Cases[kuikixml, XMLElement[{_,
"EC01"}, _, {_[{_, "ARE"}, {"idref" -> sfid_}, {}], __, _[{_,
"CN2"}, _, {twn_}], _}] :> {twn, sfid}, Infinity];さぁ、準備はできた。いよいよ描画。
In[12]:=
kuiki = Graphics[{EdgeForm[Thin], FaceForm[LightGreen],
Tooltip[#[[2]] /. sfs, #[[1]]] & /@ twns}]
あらららら。琵琶湖が...ない。
どうやら国土数値で公開されている行政区域のデータには、湖沼のデータは含まれていないらしい。湖沼は、日本全国の分をひとつにまとめた大容量(圧縮ファイルで 8.4MB)のものが公開されていた。複数の都府県にまたがる湖沼が少なからずあって、都道府県別にはできないとか、そういう事情なのかしら。今回みたいに「滋賀県の白地図を描きたい」みたいなときにはちょっと不便。
でも、まぁ、ここまで来たので、琵琶湖の抽出もやってみようじゃないの。
データにはたぶんそれぞれの湖沼の名称が含まれているだろうから、ダウンロードしたファイルを同じ階層に展開してから読み込んで、「琵琶湖」をキーにさっきと同じ要領で探してみる。XML ファイルが大きい(展開後のサイズが 183MB)ので、かなり時間がかかる。
In[13]:=
kosyouxml = Import["W09-05.xml"];
Cases[kosyouxml, _[___, _[___, _[___, _[___, "琵琶湖", ___], ___],
___], ___], Infinity]
Out[14]=
{XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, GC01}, {id -> gc01_419}, {XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, ARE}, {idref -> s419000001}, {}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, LPN}, {}, {琵琶湖}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, AAC}, {}, {25201}], XMLElement[{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app, AAC}, {}, {25202}],
.....<以下略>行政区域と同じく、idref を持つ ARE があるので、これを同じ要領で辿っていけば OK。今回は、滋賀県にある琵琶湖以外の小さな湖沼は無視して、琵琶湖のデータだけを決め打ちで拾うことにする。
In[15]:=
(* ポイントデータ *)
kosyoupts = Cases[kosyouxml, XMLElement[{_, "GM_Point"}, {"id" ->
id_}, {_[__, {_[__, {_["DirectPosition.coordinate", _, {yx_}],
_}]}]}] :> Rule[id, convCoord[yx]], Infinity];
(* 曲線データ *)
kosyoucvs = Cases[kosyouxml, _[{_, "GM_Curve"}, {"id" -> id_},
data__] :> Rule[id, Cases[data, _[kind : "GM_PointRef.point" |
"DirectPosition.coordinate", {pid___}, {yx___}] :>
If[kind === "GM_PointRef.point", pid[[2]] /. kosyoupts,
convCoord[yx]], Infinity]], Infinity];
(* 琵琶湖の面輪郭データ *)
biwakosfc = Cases[kosyouxml,
XMLElement[_, {"id" -> "s419000001"}, data__] :>
Cases[data, _[{_, "GM_SurfaceBoundary.exterior" |
"GM_SurfaceBoundary.interior"}, _, {cvdata_}] :>
Polygon[Join @@ Cases[cvdata, _[{_,
"GM_CompositeCurve.generator"},
{"idref" -> cvid_}, _] :> StringDrop[cvid, 1], Infinity] /.
kosyoucvs], Infinity], Infinity][[1]];準備完了。
In[18]:=
biwako = Graphics[{EdgeForm[Thin], FaceForm[LightBlue],
Insert[biwakosfc, FaceForm[LightGreen], 2]}]
あとは、行政区域の白地図にこれを重ねれば完成。行政区域を描くときに Tooltip を入れておいたので、マウスでポイントした区域の名称が表示される。メデタシメデタシ(?)。
In[19]:=
Show[{kuiki, biwako}]
ちなみに、滋賀県の面積のうち琵琶湖の占める割合は、約6分の1です。


