最近、社会情勢もあって白地図を目にすることが多くなっていますが、この機会にGISデータを扱う、シェープファイルやGeoJSONにふれてみようと思いました。
PythonであればシェープファイルやGeoJSONに関する便利ライブラリもいっぱいあるだろう、ひょっとするとPandasでもシェープファイルを使えるんじゃなかろかと思って軽い気持ちでググったら、GeoPandasというのがあるようです。
ということで、GeoPandasを少しさわってみました。
（この後いろいろやってみるにあたり）最初におさえておくと便利かなといういくつかのAPIをつまみ食いしてみます。
利用したGeoPandasのバージョンは0.7.0です。 また、ここではインストール方法についてはふれませんが、私の環境の都合でmatplotlibとrtreeを個別にインストールしました。
GeoPandasについて
◆公式 geopandas.org
私のこの記事の自己否定ではないですが（笑）、Pandasに慣れた方であれば、公式サイトのExample Galleryと USER GUIDEにスニペットが豊富なので、それをたどっていくのが近道だと思いました。
◆参考にさせていただいた日本語のサイト・ブログ
サンプルコード
データの前準備
この手の話でやっかいなのは、身近な例に対応したシェープファイルなどの「データ」です。
今回は、次のリンクで日本の様々なシェープファイルがダウンロードできますので、こちらを使います。
www.gsi.go.jp 2020年4月時点では、第2.2版ベクタ(2016年公開)というのがあるので、これの「全レイヤ」のZIPファイルをダウンロードしたものを使います。
※本件に限りませんが、この類のものは特に、データの権利や利用規約についてご注意ください。
サンプルコード
上記の全球地球データを展開して、そのディレクトリで以下のPythonのプログラムを起動してみてください。
import geopandas as gpd import pandas as pd from matplotlib import pyplot as plt #https://www.gsi.go.jp/kankyochiri/gm_jpn.html # 地球地図日本 # ファイル名と実際のデータの対応はGlobal Map Specifications # https://github.com/globalmaps/specifications/blob/master/gmspec-2.2.pdf # でいうと、p.3の表（だと思われる） # これによると、例えば次のようなもの。 """ geopandas.read_file(シェープファイル名) # GeoJSONもいける """ bnd = gpd.read_file('polbnda_jpn.shp') # Political Boundary Area ? riverl = gpd.read_file('riverl_jpn.shp') # Water Course ? rstatp = gpd.read_file('rstatp_jpn.shp') # Railroad Station ? bnd riverl rstatp # geometryというフィールドをみてみると、ああそんな感じなんだというところが察せられる """ なんと、to_json()で、GeoJSON形式で出力できる。つまり、シェープファイルを慣れたGeoJSON形式に出力することもできる。 もちろん、Pythonにはいくつかこの用途の有名ライブラリがあると思われるが、read_fileからのto_jsonはその中でも手軽な部類と思われる。 """ # rstatp.to_json() """ plotに対応 """ ax = bnd.plot(color='#fffacd',edgecolor='gray') plt.show() # ↓ # ひとまず、このレベルのコーディングで、bndデータを元にした、市区町村レベルの境界線の白地図が表示される。 """ geometryフィールド ← ポイント！！ geometryフィールドをdirして、なんとなく察しておく """ bnd['geometry'].head(1).geom_type #Polygon riverl['geometry'].head(1).geom_type #LineString rstatp['geometry'].head(1).geom_type #Point dir(bnd.geometry) # bnd.geometryでもアクセス可能。以降、geometryはこのアクセス方法で記載（単に好みとこの記事での目的からすると、この記法の方がgeopandasの約束事っぽい感じがするためわかりやすいのではと考えたため) """ geometryフィールドのtypeは？(ここでは説明しないが、Pandasの仕組み) ↓ GeoSeries """ # #In [241]: type(bnd.geometry) #Out[241]: geopandas.geoseries.GeoSeries """ geopandas.geoseries.GeoSeriesの便利メソッド ※いろいろあるが、データ全体のオーバービューとして私が面白いと思っているのは、以下のもの。 """ bnd.geometry.bounds #各シェープレコードごとの四隅の座標（相当) bnd.geometry.total_bounds # bndデータ全体の四隅の座標（相当)→ この例で言うとこのデータで定めるところの日本の最(東|西|南|北)端の座標 # ↓こんな感じ #In [251]: bnd.geometry.total_bounds #Out[251]: array([122.93349315, 20.42274033, 153.98686576, 45.55733109]) """ 基本的にはPandas(のDataFrame)なので、ブールインデックスによるselectが可能 """ rstatp.head() rstatp[rstatp['nam'] == 'HACHINOHE'] rstatp[rstatp['nam'] == 'HACHINOHE'].geometry # 八戸駅のデータ hachi = rstatp[rstatp['nam'] == 'HACHINOHE'].iloc[0].geometry hachi.xy # XY座標 hachi.x # X座標 hachi.wkt # WKT形式の座標出力 """ distanceメソッド（GeSeriesデータに対して、あるレコードの距離一覧 """ rstatp.distance(hachi) #全国の各駅からの八戸駅への距離一覧（単位は度数ベース？ 八戸駅と仙台駅の間が約14.12という数値だった。) """ contains intersects """ #八戸駅を含むエリア bnd[bnd.contains(hachi)] # 八戸市が戻る #利根川 tonegawa = riverl[riverl['nam'] == 'TONE G.'] print(tonegawa) #LineStringではあるが、複数レコードの模様 bnd.intersects(tonegawa.iloc[0].geometry) #複数レコードなのでひとまず、「利根川の1件目」のレコードを含む地域 """ unary_unionによるフィーチャーの結合を使い、「利根川全体」データを使う。 → 「利根川全体」に対するintersectsにより、利根川が流れる都市の一覧を取得 """ bnd[bnd.intersects(tonegawa[0:].geometry.unary_union)] """ dissolveによるフィーチャーの結合 GeoDataFrameのある列でgroup化して、geometryフィールドをgroupごとの塊のの図形に結合する。 cf.GeoSeriesをunary_unionで結合 """ # bndデータは、namフィールドに、都道府県名のローマ字が入っている。都道府県名ごとにdissolveすると、各都道府県の pref = bnd.dissolve(by='nam') pref.plot(); plt.show() """ contains再び(dissolve確認用) """ sapporo = bnd.head(1) bnd.contains(sapporo) bnd[bnd.contains(sapporo)] bnd[bnd.contains(sapporo.iloc[0].geometry)] pref[pref.contains(sapporo.iloc[0].geometry)] """ shapelyおよびgeopandas.GeoSeriesによる自作のGeoDataFrame #Polygon以外のLineStringなども同様 """ from shapely.geometry import Polygon polys1 = gpd.GeoSeries([Polygon([(0,0), (2,0), (2,2), (0,2)]),Polygon([(2,2), (4,2), (4,4), (2,4)])]) polys2 = gpd.GeoSeries([Polygon([(1,1), (3,1), (3,3), (1,3)]),Polygon([(3,3), (5,3), (5,5), (3,5)])]) df1 = gpd.GeoDataFrame({'geometry':polys1, 'foo':[1,2]}) df2 = gpd.GeoDataFrame({'geometry':polys2, 'bar':[1,2]}) df1 df2 """ sjoin (Spacial JOIN) 2つのGeoSeries（を含むDataFrame)について、「op=intersects/contains/within」の関係にある場合に、これらを結合したDataFrameを得ることができます。 例えば、川のシェープデータ一覧と都道府県のシェープデータの一覧がある場合、これらを結合して、ある都道府県を流れる川という関係を求めます。 """ # 川と都道府県の関係 _named_riverl = riverl[ ~ riverl['nam'].str.contains('UNK')] #UNKNOWN named_riverl = _named_riverl.dissolve(by='nam') # 支流ごとに別レコードなのでそのようなものをMULTILINESTRING gpd.sjoin(named_riverl,pref,how='inner',op='intersects')[['index_right']] # 東京を流れる川の一覧 gpd.sjoin(named_riverl,pref,how='inner',op='intersects').query('index_right == "Tokyo To"')[['index_right']]
上記を保存して、pythonをシェルから起動するか、ipythonやjupyterにコピペ起動すると、データの中身なども確認できます。
↓ 途中、2回ほどデータをmatplotlibで可視化したウィンドウが出ますが、そのうちの一つ
解説
以下、上記のサンプルコードの部分部分の解説です。
read_fileでシェープファイルやGeoJSONファイルの読み込み
ファイル名を指定するとあっさり読み込みできます。
bnd = gpd.read_file('polbnda_jpn.shp')
規格化されていることの強みですね〜。
オーバービューしてみます。
In [492]: bnd Out[492]: f_code coc nam laa pop ypc adm_code salb soc geometry 0 FA001 JPN Hokkai Do Sapporo Shi 1930496 2014 01100 UNK JPN POLYGON ((141.44980 43.16333, 141.44769 43.157... 1 FA001 JPN Hokkai Do Hakodate Shi 274485 2014 01202 UNK JPN POLYGON ((140.86501 42.01013, 140.86800 42.008... 2 FA001 JPN Hokkai Do Otaru Shi 127224 2014 01203 UNK JPN POLYGON ((141.24820 43.15973, 141.24680 43.158... 3 FA001 JPN Hokkai Do Asahikawa Shi 349057 2014 01204 UNK JPN POLYGON ((142.43280 43.94814, 142.43790 43.944... 4 FA001 JPN Hokkai Do Muroran Shi 91276 2014 01205 UNK JPN POLYGON ((140.99080 42.43800, 140.99680 42.436... ... ... ... ... ... ... ... ... ... ... ... 2909 FA001 JPN Tokyo To UNK -99999999 0 UNK UNK JPN POLYGON ((139.80595 35.58558, 139.79480 35.592... 2910 FA001 JPN Tokyo To UNK -99999999 0 UNK UNK JPN POLYGON ((140.05172 31.43784, 140.05053 31.437... 2911 FA001 JPN Tokyo To UNK -99999999 0 UNK UNK JPN POLYGON ((140.29860 30.47233, 140.28740 30.482... 2912 FA001 JPN Tokyo To UNK -99999999 0 UNK UNK JPN POLYGON ((140.34214 29.79279, 140.34153 29.792... 2913 FA001 JPN Aichi Ken UNK -99999999 0 UNK UNK JPN POLYGON ((136.81190 34.99609, 136.80595 35.000... [2914 rows x 10 columns] In [495]: riverl Out[495]: f_code hyc lit nam soc geometry 0 BH140 8 1 UNK JPN LINESTRING (127.73975 26.35407, 127.74320 26.3... 1 BH140 8 1 UNK JPN LINESTRING (127.74320 26.35579, 127.73947 26.3... 2 BH140 8 1 UNK JPN LINESTRING (127.79680 26.35945, 127.80429 26.3... 3 BH140 8 1 UNK JPN LINESTRING (127.74320 26.35579, 127.75616 26.3... 4 BH140 8 1 UNK JPN LINESTRING (127.76031 26.36250, 127.75810 26.3... ... ... ... ... ... ... ... 4790 BH140 8 2 ONUMA JPN LINESTRING (141.77698 45.38849, 141.77331 45.3... 4791 BH140 8 1 KOETOI G. JPN LINESTRING (141.77080 45.40520, 141.76999 45.4... 4792 BH140 8 1 UNK JPN LINESTRING (141.85770 45.43480, 141.86020 45.4... 4793 BH140 8 1 UNK JPN LINESTRING (148.61110 45.49347, 148.61604 45.4... 4794 BH140 8 2 UNK JPN LINESTRING (136.91471 36.12237, 136.91835 36.1... [4795 rows x 6 columns] In [496]: rstatp Out[496]: f_code nam soc geometry 0 AQ125 HACHINOHE JPN POINT (141.43359 40.51220) 1 AQ125 NINOHE JPN POINT (141.28420 40.25346) 2 AQ125 IWATENUMAKUNAI JPN POINT (141.22236 39.96242) 3 AQ125 MORIOKA JPN POINT (141.14149 39.69653) 4 AQ125 OMAGARI JPN POINT (140.48230 39.46413) .. ... ... ... ... 103 AQ125 SHICHINOHETOWADA JPN POINT (141.15408 40.71989) 104 AQ125 SHINAOMORI JPN POINT (140.69348 40.82801) 105 AQ125 OKUTSUGARUIMABETSU JPN POINT (140.51524 41.14521) 106 AQ125 KIKONAI JPN POINT (140.43393 41.67762) 107 AQ125 SHINHAKODATEHOKUTO JPN POINT (140.64658 41.90539) [108 rows x 4 columns]
↑ 前者が行政区の境界線、川、鉄道駅のデータです。
シェープファイルもGeoJSONも良くわかっていませんが、geometryフィールドに格納されているPOLYGON、LINESTRING、POINTといったそれらしきデータ形式なんだろうなというところが察せられます。
また、「geometry」は予約語、...なのかはわかりませんが、GeoSeriesというこれまたPandasのSeriesの拡張らしきデータ型のようです。
plot してみる
Pandasの血を引くだけあってか、matplotlibでの可視化に対応しています。DataFrame(この場合は、GeoDataFrame).plot()でラフにデータを可視化しながら遊ぶことで、理解が深まるような気がします。
ax = bnd.plot(color='#fffacd',edgecolor='gray') plt.show() # ↓ # ひとまず、このレベルのコーディングで、bndデータを元にした、市区町村レベルの境界線の白地図が表示される。
to_json() : GeoJSONを出力
サンプルプログラムではコメントアウトしていますが、GeoDataFrameのメソッドでto_json()というものがありますが、こいつは、いわゆるGeoJSON形式のデータを出力できます。
私がこの関数を初めて試した時に出力されたGeoJSONファイルをそのまま、これまた私が良くわかっていないGeoJSONを取り扱える別のソフトに読み込ませてみたところ、特にデータ内容・形式の知識がなくても当たり前のように取り扱いできました。
シェープファイル形式のデータでしか情報公開されていないデータも多くありますので、GeoPandasのread_fileでシェープファイルを読み込み、to_jsonでGeoJSONを出力して、あとはGeoJSONの扱いが得意なWeb系にデータを引き渡すツールとしての使い方もあると思います。
GeoSeriesの便利メソッド
geometryフィールドに格納された、シェープデータ(GeoSeries型)に関する各種演算が可能です。
面積を求めたり、シェープデータどおしの包含関係を調べたり、複数のレコードを結合して大きな形状を得たり、座標どおしの関係性に基づいて結合したり、上記のデータ例で言うと、都道府県名に「HOKKAIDO」を持つ行政区（市区町村、郡もかな）を結合して、「北海道」全体のデータを生成したりといった、こんなこともできるのかすげーなてなこともできます。
サンプルコードの中に説明もつけているので、あらためて見てみてください。
なお、これらの便利メソッド名だけもう一度だけ箇条書きしておきます。
- bounds
- total_bounds
- wkt (Well-Known Text) Well-known text - Wikipedia
- contains、intersects .... : 「geometry」フィールドを元にしたシェープデータ包含関係の演算（ ex. 多摩川が東京都を横切るか? ）
- dissolve (「geometry」ではなく、他のフィールドを元にした、groupby相当〜 結合されたシェープデータを生成）
- sjoin : 精度の課題と使いどころに気をつける必要はありそうですが、これは特にすごいですね。↓
最後の sjoin(Spacial JOIN) について、サンプルプログラム中の解説コメントを再掲しておきます。
2つのGeoSeries（を含むDataFrame)について、「op=intersects/contains/within」の関係にある場合に、これらを結合したDataFrameを得ることができます。 例えば、川のシェープデータ一覧と都道府県のシェープデータの一覧がある場合、これらを結合して、ある都道府県を流れる川という関係を求めます。
