Helsinki大学のAutomating GIS 2022(https://autogis-site.readthedocs.io/en/latest/index.html)に従ってGeopandasについて学習していたところ、途中で教材になっているフィンランドのGISデータが取得できなかったため、「国土数値情報（https://nlftp.mlit.go.jp/ksj/index.html）」にある「政策区域－行政区域」情報（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html）から「京都府」のShapeFileデータをダウンロード、そして、GoogleDriveにアップロードして、Geopandasで加工・可視化してみました。
与えられた教材を使うのとは異なり、色々と工夫しながら学習を進めていけるので、とても理解が深まります。
まず、「国土数値情報」のHPに行き、左にあるメニューから「2.政策区域」を選択（クリック）し、「行政地域」メニューから「行政区域（ポリゴン）」を選択します。
行政区域データの説明およびダウンロード画面に飛ぶので、画面を下方向にスクロールし、お好みの都道府県を選択して、ShapeFileをダウンロードします。今回、わたしは「京都府」を選択してみました。
ダウンロードしたデータフォルダーを解凍した上で、GoogleDriveのGoogleColabo Notebooksのフォルダにアップロードしてます。
これで準備は完了！
GoogleColaboratoryを立ち上げます。
まず、Geopandasをインポートするとともに、GoogleColaboratoryでGoogleDriveを利用するための以下のおまじないコードを実行します。
!pip install geopandas
import geopandas as gpd
#Google ColabでGoogle Driveを使用するときのおまじない
from google.colab import drive
drive.mount('/content/drive')
早速、「京都府」のShapeFileをGeopandasで読み込みます。
#Define file path
fp = ('/content/drive/MyDrive/Colab Notebooks/N03-20220101_26_GML/N03-22_26_220101.shp')
#Read in the data from the file (startting at row 9)
data = gpd.read_file(fp)
data.head()
以下、ShapeFileの中身が表示されます。
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 0 | 京都府 | None | None | None | None | MULTIPOLYGON (((135.31410 35.45566, 135.31255 ... |
| 1 | 京都府 | None | 京都市 | 北区 | 26101 | POLYGON ((135.72539 35.17080, 135.72552 35.170... |
| 2 | 京都府 | None | 京都市 | 上京区 | 26102 | POLYGON ((135.75062 35.03822, 135.75079 35.038... |
| 3 | 京都府 | None | 京都市 | 左京区 | 26103 | POLYGON ((135.80481 35.31708, 135.80586 35.316... |
| 4 | 京都府 | None | 京都市 | 中京区 | 26104 | POLYGON ((135.73195 35.02251, 135.73195 35.022... |
以下のコードを実行して、座標系を確認してみます。
data.crs
出力です。
<Geographic 2D CRS: EPSG:6668> Name: JGD2011 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: Japan - onshore and offshore. - bounds: (122.38, 17.09, 157.65, 46.05) Datum: Japanese Geodetic Datum 2011 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
座標系はJGD2011（詳しいことはここでは割愛）、この座標系のEPSGコードは6668と表示されます。
ちょっとおためしで、このデータに含まれている座標情報を覗いてみました。以下のコードを実行します。
data.geometry.head()
出力はこんな感じで、各自治体のpolygonを定義する経度・緯度が並んでいます。
0 MULTIPOLYGON (((135.31410 35.45566, 135.31255 ... 1 POLYGON ((135.72539 35.17080, 135.72552 35.170... 2 POLYGON ((135.75062 35.03822, 135.75079 35.038... 3 POLYGON ((135.80481 35.31708, 135.80586 35.316... 4 POLYGON ((135.73195 35.02251, 135.73195 35.022... Name: geometry, dtype: geometry
colmunの名称がわかりにくいので、名前を変更します。以下のコードを実行します。
data = data.rename(
columns={
'N03_001': 'PREF',
'N03_003': 'CITY',
'N03_004': 'TOWN',
'N03_007': 'CODE'
} )
data.head()
出力すると、columnの名称が置き換えられています。
| PREF | N03_002 | CITY | TOWN | CODE | geometry | AREA | |
|---|---|---|---|---|---|---|---|
| 0 | 京都府 | NaN | NaN | NaN | NaN | MULTIPOLYGON (((135.31410 35.45566, 135.31255 ... | 0.456821 |
| 1 | 京都府 | NaN | 京都市 | 北区 | 26101 | POLYGON ((135.72539 35.17080, 135.72552 35.170... | 0.009380 |
| 2 | 京都府 | NaN | 京都市 | 上京区 | 26102 | POLYGON ((135.75062 35.03822, 135.75079 35.038... | 0.000694 |
| 3 | 京都府 | NaN | 京都市 | 左京区 | 26103 | POLYGON ((135.80481 35.31708, 135.80586 35.316... | 0.024415 |
| 4 | 京都府 | NaN | 京都市 | 中京区 | 26104 | POLYGON ((135.73195 35.02251, 135.73195 35.022... | 0.000732 |
PREF列に含まれているデータをuniqueを使って取り出してみます。
data['PREF'].unique()
出力は以下のとおりで、京都府のデータですから都道府県列は「京都府」ということになりました。
array(['京都府'], dtype=object)
次はCITY列で同じことをしてみました。
data['CITY'].unique()
京都府内にある京都市と郡名が出力されました。あれ？舞鶴市とか福知山市がでてきません。それはあとで解決します。
array([None, '京都市', '乙訓郡', '久世郡', '綴喜郡', '相楽郡', '船井郡', '与謝郡'],
dtype=object)このデータを図化してみました。
data.plot()
つぎは、行ごとに面積を計算してみます。
for index, row in data[0:5].iterrows():
polygon_area = row['geometry'].area
print(f'The polygon in row {index} has a surface area of {polygon_area:.4f}.')
一番最初の行が京都府全体になりますが、0.4568（ネットでは京都府の面積は4612平方キロメートルとあります）という訳のわからない数字が出てきました。JGD2011では面積は面積は正確に計算できませんし、メートル単位ではないということが理解できました。
The polygon in row 0 has a surface area of 0.4568. The polygon in row 1 has a surface area of 0.0094. The polygon in row 2 has a surface area of 0.0007. The polygon in row 3 has a surface area of 0.0244. The polygon in row 4 has a surface area of 0.0007.
※実用精度の面積の計算方法は別途メモしようと思います。
面積の単位や精度は一度脇において、GeoDataFlameの処理に戻ります。新たに設定した「Area」列に、各行の面積を計算して入力します。
data['AREA'] = data.area
data.head()
各行の最終列に面積らしき数値が入りました。
Area列の統計諸量を表示してみます。
#面積の平均値、最大値、最小値などを表示
data['AREA'].describe()
count 1.000000e+03 mean 9.136425e-04 std 1.494483e-02 min 9.440544e-11 25% 4.549586e-09 50% 9.340953e-09 75% 2.564530e-08 max 4.568212e-01 Name: AREA, dtype: float64
各数値をみても、やはり数字の桁数がメートル系とは全く異なっています。
次は、京都市だけをピックアップして、「kyoto_city」に入力します。
#京都市のみを抽出
kyoto_city = data[data.CITY == '京都市']
kyoto_city.tail()
以下、出力です。tail()を使っているので、表示データは末尾からのピックアップになっています。
| PREF | N03_002 | CITY | TOWN | CODE | geometry | AREA | REGION | |
|---|---|---|---|---|---|---|---|---|
| 7 | 京都府 | NaN | 京都市 | 南区 | 26107 | POLYGON ((135.72575 34.98546, 135.72577 34.985... | 0.001560 | 京都市 |
| 8 | 京都府 | NaN | 京都市 | 右京区 | 26108 | POLYGON ((135.72989 35.26630, 135.72995 35.266... | 0.028891 | 京都市 |
| 9 | 京都府 | NaN | 京都市 | 伏見区 | 26109 | POLYGON ((135.78157 34.97374, 135.78224 34.973... | 0.006084 | 京都市 |
| 10 | 京都府 | NaN | 京都市 | 山科区 | 26110 | POLYGON ((135.81113 35.01549, 135.81151 35.015... | 0.002833 | 京都市 |
| 11 | 京都府 | NaN | 京都市 | 西京区 | 26111 | POLYGON ((135.65553 35.02394, 135.65551 35.023... | 0.005848 | 京都市 |
kyoto_cityを使って、京都市だけを図化してみました。
kyoto_city.plot()
おっ、この大きさだと、各行政地域の境界が識別できますね。
さて、このへんで行方不明のままの福知山市や舞鶴市を探そうと思います。データベースの最初の方は京都市なので、その後ろくらいを覗いてみました。
data[10:15]
京都市のうしろの行に福知山市出現。２列目はNoneで12行目の３列目に福知山市、４列目以降に舞鶴市が出現しました。
| PREF | N03_002 | CITY | TOWN | CODE | geometry | AREA | |
|---|---|---|---|---|---|---|---|
| 10 | 京都府 | NaN | 京都市 | 山科区 | 26110 | POLYGON ((135.81113 35.01549, 135.81151 35.015... | 2.833339e-03 |
| 11 | 京都府 | NaN | 京都市 | 西京区 | 26111 | POLYGON ((135.65553 35.02394, 135.65551 35.023... | 5.847774e-03 |
| 12 | 京都府 | NaN | NaN | 福知山市 | 26201 | POLYGON ((135.12664 35.46949, 135.12692 35.469... | 5.477151e-02 |
| 13 | 京都府 | NaN | NaN | 舞鶴市 | 26202 | POLYGON ((135.31410 35.45566, 135.31255 35.455... | 1.628898e-07 |
| 14 | 京都府 | NaN | NaN | 舞鶴市 | 26202 | POLYGON ((135.33622 35.47246, 135.33616 35.472... | 3.938604e-09 |
データベースの構造は１列目が都道府県名、２列目が政令指定都市名or郡名、３列目が一般の市名or町名のようです。
データベースの内容が少しずつ分かってきました。
市町村名を大ぐくりにするために、市域名と郡域名で整理してみます。
’REGION'という列を追加して、もし’CITY’がNaN（isna()を使用します）ならば、’TOWN'の文字列を、’CITY'に有意な文字列があれば’CITY'の文字列を、’REGION'列に入力します。
import pandas as pd
data["REGION"] = "" #"REGION"をstring型として仮設定
for index,row in data.iterrows():
if pd.isna(row["CITY"]):
data.at[index, "REGION"] = row["TOWN"]
else:
data.at[index, "REGION"] = row["CITY"]
data[10:15]
出力は以下のようになります。表示は市域ばかりですので、’REGION'は市の名前ばかりです。京都市や舞鶴市は複数行に分かれています。
| PREF | N03_002 | CITY | TOWN | CODE | geometry | AREA | REGION | |
|---|---|---|---|---|---|---|---|---|
| 10 | 京都府 | NaN | 京都市 | 山科区 | 26110 | POLYGON ((135.81113 35.01549, 135.81151 35.015... | 2.833339e-03 | 京都市 |
| 11 | 京都府 | NaN | 京都市 | 西京区 | 26111 | POLYGON ((135.65553 35.02394, 135.65551 35.023... | 5.847774e-03 | 京都市 |
| 12 | 京都府 | NaN | NaN | 福知山市 | 26201 | POLYGON ((135.12664 35.46949, 135.12692 35.469... | 5.477151e-02 | 福知山市 |
| 13 | 京都府 | NaN | NaN | 舞鶴市 | 26202 | POLYGON ((135.31410 35.45566, 135.31255 35.455... | 1.628898e-07 | 舞鶴市 |
| 14 | 京都府 | NaN | NaN | 舞鶴市 | 26202 | POLYGON ((135.33622 35.47246, 135.33616 35.472... | 3.938604e-09 | 舞鶴市 |
groupby()を活用して、’REGION'の内容ごとに整理してみます。
grouped_data = data.groupby("REGION")
grouped_data.head()
以下、少し難解な出力になるので省略します。おそらくグループごとに先頭５行程度が出力されているのだと理解します。
グループ化されたデータを図化してみました。これで市域・郡域が整理された図が出てくると思ったら、大違い。市域・郡域別の図が出力されました。逆に面白い出力だったので、掲載しておきます。
grouped_data.plot(column="REGION", legend=True, cmap='tab20')
気になるのは凡例が文字化けしていることです。
よくよく考えると、grouped()したデータでは無く、grouped()しない生データを使えばうまく図化してくれるのではないかという結論に至り、以下のコードを実行しました。なんと結果は思い通りになりました。
data.plot(column="REGION", legend=True, cmap='tab20')
残念なのは凡例が文字化けしたママであるということと、凡例が図と重なっていて、図が見えないという点です。
matplotlibおよびgeopandasで日本語が文字化けするという現象はよくあるようで、ネット上にたくさんのヒントが落ちていました。一番効率がよさそうなのは、ズバリ、以下のコードだとわたしは思います。
#matplotlibで日本語を表示するためのおまじない
!pip install japanize_matplotlib
import japanize_matplotlib
これを冒頭において先ほどの図化コードを実行します。
#matplotlibで日本語を表示するためのおまじない
!pip install japanize_matplotlib
import japanize_matplotlib
data.plot(column="REGION", legend=True, cmap='tab20')
ごいごいすー！
このおまじないで日本語が表示できました。
次に凡例(legend)と図の重なりを解消するために、凡例をグラフ外に出す方法を検索しました。しかし、matplotlib側のコードは出ているのですが、df.plot(legend=True)で図化した場合の凡例位置の命令方法は探し出すことができませんでした。
なんとか方法はないものかと悪戦苦闘していたところ、グラフのサイズを変えると、図と凡例の重なりは回避できることを発見。根本的な解決にはなりませんが、まあ、実用上は許容範囲でしょうか？
data.plot(column="REGION", cmap='tab20', figsize=(10,10), legend=True)
ちなみに、例えば与謝郡だけを取り出して図化するならば、以下のようなコードだと思います。
yosa_cou = grouped_data.get_group("与謝郡")
yosa_cou.plot()
コメント
コメントを投稿