【PythonでGIS】GeoPandasまとめ

はじめに

本記事は、PythonでGISデータを扱うためのライブラリGeoPandasを紹介する記事です
PythonにおけるGISデータの操作全般はこちらの記事で紹介しているので、まずはリンク先をご一読頂ければと思います。

GeoPandasとは

GIS（地理情報システム）データをPythonで扱うためのライブラリです。

※下図はGeoPandasで日本の堤高100m以上のダムを可視化した例

こちらの記事でも紹介しましたが、PythonでのGISデータ操作は多数のライブラリが乱立しています。
GeoPandasはこれらの中でも、データをPandasのDataFrameチックに扱えることが特徴のライブラリです

GeoPandasを選択するメリット

数あるライブラリの中でもGeoPandasは、人気が急上昇しています。
（2021年7月現在のGitHub Star数=2700、2年前から倍増）

個人的に感じたGeoPandasを選択するメリットは、以下となります

・Pandasに慣れ親しんだ人であれば、直感的に操作できる
・非常に多機能なため、入出力、表示、計算など多くの機能を一つのライブラリに統一できる
・座標変換やポリゴン操作などの複雑な処理も、簡単に実装できる
・公式ドキュメントが充実している（英語ですが・・）

実際に使ってみると、こちらの記事の方法を使うよりもシンプルに実装できる処理が多く、非常に使い勝手の良いライブラリだと思いました。

実装手順

実際にShapefileを読み込み、よく使う処理をユースケースとしてひととおり実装してみます。
使用コードは**GitHubにもアップロードしております**

ユースケース一覧

こちらの記事と同じユースケースを、GeoPandasで実装します
また、GeoPandasではジオメトリが自動でShapelyのPoint、LineString、Polygonクラスに格納されるので、
ポイント・ライン・ポリゴンデータへの変換は必要ありません（ユースケースでも紹介しません）

読込

No	内容	追加で使用するライブラリ	実装法リンク
1	Shapefileの読込	-	リンク
2	ジオメトリ情報と属性情報の取得	-	リンク
3	GeoJSONの読込	-	リンク

処理

No	内容	追加で使用するライブラリ	実装法リンク
1	座標変換	-	リンク
2	(ポイントデータへの変換は不要)	-	-
3	ポイント間の距離測定	pyproj	リンク
4	最も近いポイントを探す	scikit-learn	リンク
5	(ラインデータへの変換は不要)	-	-
6	ラインの長さ測定	-	リンク
7	(ポリゴンデータへの変換は不要)	-	-
8	ポリゴンの重心測定	-	リンク
9	ポリゴンの面積測定	-	リンク
10	名称・住所から緯度経度検索 （ジオコーディング）	-	リンク
11	緯度経度から住所検索 （逆ジオコーディング）	-	リンク

保存

No	内容	追加で使用するライブラリ	実装法リンク
1	ポイントのファイル出力	-	リンク
2	ラインのファイル出力	-	リンク
3	ポリゴンのファイル出力	-	リンク

表示

No	内容	使用するライブラリ	実装法リンク
1	ポイントの表示	shapely＋japanmap＋matplotlib	リンク
2	ラインの表示	shapely＋japanmap＋matplotlib	リンク
3	ポリゴンの表示	shapely＋japanmap＋matplotlib	リンク

使用するデータ

こちらの記事と同じデータを使用して、実装を進めます

必要ライブラリのインストール

トラブルシューティング含め、こちらの記事を参照ください

具体的な実装例

Shapefileの読込、処理、保存、表示に分けて、
各ユースケースの実装法を解説します

読込1：Shapefileの読込

Shapefileの読込には、read_fileメソッドを使用します。

Shapefile（ポイント、ライン、ポリゴン）の読込

# 必要ライブラリの読込（読込以外で使用するライブラリも含みます）
import geopandas as gpd
from shapely.geometry import Point
import pyproj
import pandas as pd
import numpy as np
import re

# 入力ファイルのパス
DAM_PATH = './W01-14_GML/W01-14-g_Dam.shp'  # 国交省ダムデータ（ポイント、https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-W01.html）
RIVER_PATH = './W05-09_25_GML/W05-09_25-g_Stream.shp'  # 国交省滋賀県河川データ（ライン、https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-W05.html）
LAKE_PATH = './W09-05_GML/W09-05-g_Lake.shp'  # 国交省湖沼データ（ポリゴン、https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-W09-v2_2.html）
FARM_PATH= './01206釧路市（2020公開）/01206釧路市（2020公開）_5.shp'  # 農水省農場データ（ポリゴン、https://www.maff.go.jp/j/tokei/porigon/hudeporidl.html）

# shpファイル読込
# ポイント（ダムデータ）
gdf_dam = gpd.read_file(DAM_PATH, encoding='cp932')  # Shapefile読込
print(f'number of records = {len(gdf_dam)}')  # レコード数
# ライン（河川データ）
gdf_river = gpd.read_file(RIVER_PATH, encoding='cp932')  # Shapefile読込
# ポリゴン（湖沼データ）
gdf_lake = gpd.read_file(LAKE_PATH, encoding='cp932')  # Shapefile読込
# ポリゴン（農場データ）
gdf_farm = gpd.read_file(FARM_PATH, encoding='cp932')  # Shapefile読込

# 読み込んだダムデータの中身を表示
gdf_dam

表示

number of records = 2749

	W01_001	W01_002	W01_003	W01_004	W01_005	W01_006	W01_007	W01_008	W01_009	W01_010	W01_011	W01_012	W01_013	W01_014	geometry
0	駒込	191	堤川	駒込川	7	1,2,6	84.5	270.0	320	7800	4	-	青森県青森市大字駒込	4	POINT (140.88589 40.70977)
1	下湯	192	堤川	堤川	12	1,2,4	70.0	783.5	3723	12600	4	1988	青森県青森市大字荒川字横倉	1	POINT (140.77942 40.69178)
2	浅虫	193	浅虫川	浅虫川	7	1,2	9.0	215.0	10	300	4	2002	青森県青森市大字浅虫字山下	1	POINT (140.87188 40.88522)
3	田の沢溜池	194	小湊川	小湊川	3	3	21.0	205.0	124	1174	5	1945	青森県東津軽郡平内町大字松野木	1	POINT (140.99147 40.89848)
4	清水目	195	野辺地川	野辺地川	7	1	33.5	195.0	75	2630	4	2000	青森県上北郡東北町字清水目	1	POINT (141.07611 40.81095)
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2744	カンジン	2748	仲地川	仲地川	-	3	57.6	1070.0	18	1580	4	2005	沖縄県島尻郡久米島町	4	POINT (126.73585 26.36093)
2745	倉敷脇（再）	2740	比謝川	与那原川	12	1,2,4	15.0	200.0	80	7100	4	1994	沖縄県うるま市石川楚南	2	POINT (127.80719 26.39166)
2746	漢那脇	2733	漢那福地川	漢那福地川	12	1,2,3,4	37.0	500.0	991	8200	2	1993	沖縄県国頭郡宜野座村字漢那地先	2	POINT (127.95243 26.48487)
2747	安波脇	2711	安波川	安波川	12	1,2,4,5	32.0	77.0	74	18600	2	1982	沖縄県国頭郡国頭村字安波小字川瀬原地先	2	POINT (128.26799 26.70919)
2748	安富祖	2726	安富祖川	安富祖川	12	1	37.5	158.1	197	355	5	2018	沖縄県国頭郡恩納村大字安富祖小字上原	4	POINT (127.88598 26.48739)

pyshpでの読込ではジオメトリと属性を別々に読み込む必要がありましたが、
GeoPandasではジオメトリと属性を一緒にDataFrame（正確にはGeoDataFrame）形式で読込でき、実装を簡略化することができます。

また、読み込んだジオメトリは'geometry'列にshapelyのクラス（Point、Line、Polygon）として格納されるため、距離や面積などの算出を簡単に行う事ができます**（詳しくは後述）

読込2：ジオメトリ情報と属性情報を1レコードずつ取得

1レコードずつデータにアクセスしたい場合、PandasのDataFrameと同様にiterrowsメソッドを利用するのが便利です。
iterrowsメソッドは、for文と組み合わせることで、インデックスとレコード（行ごとのデータ）を一緒に取得できます

ジオメトリ情報と属性情報を1レコードずつ取得

# ポイント
for i, row in gdf_dam.iterrows():
    print(f'ポイント位置{row["geometry"]}')  # ジオメトリ情報
    print(f'ダム名:{row["W01_001"]}\n\
          堤高:{row["W01_007"]}m\n\
          総貯水量:{row["W01_010"]}千m3\n\
          所在地:{row["W01_013"]}')  # 属性情報

# ライン
for i, row in gdf_river.iterrows():
    print(f'ライン位置{row["geometry"]}')  # ジオメトリ情報
    print(f'河川名:{row["W05_004"]}')  # 属性情報

# ポリゴン
for i, row in gdf_lake.iterrows():
    print(f'ポリゴン位置{row["geometry"]}')  # ジオメトリ情報
    print(f'湖沼名:{row["W09_001"]}\n\
          最大水深:{row["W09_003"]}千m3\n\
          水面標高:{row["W09_004"]}')  # 属性情報

読込3：GeoJSONの読込

pyshpではShapefile形式の読込しかできませんでしたが、GeoPandasではGeoJSONも読み込むことができます。
Shapefileの読込と同様に、read_fileメソッドを使用します。

GeoJSON読込

dam_over100m_path = './dams_over100m.geojson'
gdf_dam_over100m = gpd.read_file(dam_over100m_path, encoding='cp932')

処理1：座標変換

こちらの記事で紹介しましたが、GeoPandasにおいても座標変換は非常に重要なテクニックです。

座標変換の前提知識として、まずはこちらの節を一読頂ければと思います

一括座標変換（変換前座標を取得できるとき）

osgeoやpyprojライブラリでは、1レコードごとに座標変換を実施（一括変換はリスト内包表記を使用）していましたが、GeoPandasではto_crsメソッドで一括座標変換を実現できます。

読み込んだShapefileに最初から座標系情報が含まれているときは、以下のようにto_crsメソッドに変換後のEPSGコードを指定するだけで、座標変換を実現できます（ポイント、ライン、ポリゴン全て同様です）

農場データを一括座標変換

# 変換後の座標系指定（平面直角座標13系(EPSG2455) → 緯度経度(EPSG4612)）
dst_proj = 4612 # 変換後の座標系を指定

# 座標変換
gdf_farm_transfer = gdf_farm.to_crs(epsg=dst_proj)  # 変換式を作成

for (i1, row1), (i2, row2) in zip(gdf_farm.iterrows(), gdf_farm_transfer.iterrows()):
    print(f'変換前ポリゴン位置 {row1["geometry"]}')  # 位置情報（座標変換前）
    print(f'変換後ポリゴン位置 {row2["geometry"]}')  # 位置情報（座標変換後）

変換前ポリゴン位置 POLYGON ((-12485.05864419129 -77604.35297655607, -12477.82669974685 -77605.23492100052, -12472.71142196907 -77609.99742100052, -12467.47857700221 -77618.30578547425, -12457.36558863574 -77635.92658766719, -12430.50844688025 -77695.9594842001, -12429.12327233178 -77701.77815059376, -12431.23993899844 -77703.54203948265, -12471.10382788733 -77713.06703948266, -12474.12253308018 -77710.23040711162, -12480.17857700221 -77696.19911880759, -12487.37524366887 -77666.14245214092, -12502.52114419129 -77628.16547655607, -12502.52114419129 -77623.75575433385, -12500.7572553024 -77619.34603211163, -12492.03191033554 -77609.27467436314, -12485.05864419129 -77604.35297655607))
変換後ポリゴン位置 POLYGON ((144.0961165921222 43.30135171609329, 144.0962057484647 43.30134389681651, 144.0962689041232 43.30130111011304, 144.0963335891383 43.30122640587214, 144.0964586342828 43.30106795289392, 144.0967910136917 43.30052798773086, 144.0968082177447 43.30047563164845, 144.096782169308 43.30045971838655, 144.0962910556334 43.30037331627315, 144.096253785255 43.30039880129993, 144.0961788250708 43.30052500894768, 144.0960894422639 43.30079545618248, 144.0959019018724 43.30113706905811, 144.0959018016179 43.30117676485103, 144.0959234418833 43.30121648993028, 144.0960307558965 43.30130729589896, 144.0961165921222 43.30135171609329))
:
:

メートル単位の平面直角座標（EPSG2455）から、緯度経度座標（EPSG4612）に変換されていることが分かります

一括座標変換（変換前座標を取得できない）

Shapefileに座標系情報が含まれていないときは、to_crsメソッドでの座標変換を実施する前に、crsプロパティに変換前の座標系を指定する必要があります。

ダム、河川、湖沼データを一括座標変換

# 変換前後の座標系指定（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
src_proj = 4612  # 変換前の座標系を指定
dst_proj = 3099  # 変換後の座標系を指定

# ポイント（ダムデータ、TransformPointの引数は緯度,経度の順番で指定）
gdf_dam_utm = gdf_dam.copy()
gdf_dam_utm.crs = f'epsg:{src_proj}'  # 変換前座標を指定
gdf_dam_utm = gdf_dam_utm.to_crs(epsg=dst_proj)  # 変換後座標に変換

# ライン（河川データ）
gdf_river_utm = gdf_river.copy()
gdf_river_utm.crs = f'epsg:{src_proj}'  # 変換前座標を指定
gdf_river_utm = gdf_river_utm.to_crs(epsg=dst_proj)  # 変換後座標に変換

# ポリゴン（湖沼データ）
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{src_proj}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=dst_proj)  # 変換後座標に変換

処理3：ポイントデータ間の距離測定

こちらの記事と同様に、平面座標系のときと緯度経度座標系のときで大きく方法が異なるので、ご注意ください

平面座標系のとき

'geometry'列にdistanceメソッドを適用することで、距離が測定できます
（1レコードずつ求める必要があります）

ポイント間の距離を測定（平面座標系）

# 2点が格納されたGeoDataFrame作成
gdf_new = gpd.GeoDataFrame(crs = 'epsg:3099')  # UTM座標53N系(EPSG3099)のGeoDataFrameを作成
gdf_new['geometry'] = None
gdf_new.at[0, 'geometry'] = Point(0, 0)
gdf_new.at[1, 'geometry'] = Point(1, 1)
print(gdf_new)
# 距離測定
dist = gdf_new.at[0, 'geometry'].distance(gdf_new.at[1, 'geometry'])  # 0行目と1行目の距離を測定
print(f'距離={dist}')

                  geometry
0  POINT (0.00000 0.00000)
1  POINT (1.00000 1.00000)
距離=1.4142135623730951

緯度経度座標系のとき

緯度経度座標系においては、こちらの記事と同様に、以下の2通りの方法で距離算出できます

1. UTM座標に変換してPoint.distanceメソッドで距離算出
2. pyproj.Geodメソッドを使用して距離算出

1の方法は、GeoPandasにおいては距離の起点、終点それぞれの座標変換が必要となり、一括座標変換のメリットを活かせないため、
今回は2の方法（pyproj.Geod）での実装法を紹介します

こちらの記事と同様に、こちらのデータを使用して、各ダムから都道府県庁までの距離を求めてみます

データの読込

# 県庁所在地の座標を読込
df_prefecture = pd.read_csv('./prefecture_location.csv', encoding='cp932')
# 緯度経度の分秒を小数に変換
df_prefecture['longitude'] = df_prefecture['都道府県庁 経度'].apply(lambda x: float(x.split('°')[0])
                                                                   + float(x.split('°')[1].split("'")[0]) / 60
                                                                   + float(x.split("'")[1].split('"')[0]) / 3600)
df_prefecture['latitude'] = df_prefecture['都道府県庁 緯度'].apply(lambda x: float(x.split('°')[0])
                                                                  + float(x.split('°')[1].split("'")[0]) / 60
                                                                  + float(x.split("'")[1].split('"')[0]) / 3600)
# 県庁所在地の緯度経度をDictionary化
dict_pref_office = {row['都道府県']: (row['longitude'], row['latitude']) for i, row in df_prefecture.iterrows()}

各ダムから都道府県庁までの距離をpyproj.Geodで算出

# 距離測定用のGRS80楕円体
grs80 = pyproj.Geod(ellps='GRS80')

# 都道府県と都道府県庁の位置を紐づけ
gdf_dam['prefecture'] = gdf_dam['W01_013'].apply(
    lambda x: re.match('..*?県|..*?府|東京都|北海道', x).group())
gdf_dam['pref_office_point'] = gdf_dam['prefecture'].apply(
    lambda x: Point(*dict_pref_office[x]))

# ダムデータを1点ずつ走査
for i, row in gdf_dam.iterrows():
    # 距離を計算（pyprojライブラリ使用）
    azimuth, bkw_azimuth, dist = grs80.inv(row['geometry'].x, row['geometry'].y,
                                           row['pref_office_point'].x, row['pref_office_point'].y)
    print(f'{row["W01_001"]}ダム {row["prefecture"]}庁まで{dist/1000}km')

駒込ダム 青森県庁まで17.717080913897718km
:
中城ダム 沖縄県庁まで14.200787571626435km
石垣ダム 沖縄県庁まで406.95145701089456km

上記の都道府県庁の紐づけにはapplyメソッドを使用しています。
PandasのDataFrameと同様に、一括処理にはapplyメソッドが便利です

距離を一括算出

「距離の最大値」のような統計値を計算したいとき、applyメソッドで距離を一括算出すると便利です。
下の例では、pyplot.Geodで緯度経度座標における都道府県庁からの距離を一括算出し、最も遠いダムを求めています。

県庁までの距離を一括算出し、最も遠いダムを表示

# 距離測定用のGRS80楕円体
grs80 = pyproj.Geod(ellps='GRS80')

# 都道府県と都道府県庁の位置を紐づけ
gdf_dam['prefecture'] = gdf_dam['W01_013'].apply(
    lambda x: re.match('..*?県|..*?府|東京都|北海道', x).group())
gdf_dam['pref_office_point'] = gdf_dam['prefecture'].apply(
    lambda x: Point(*dict_pref_office[x]))
# 距離を計算（pyproj.Geod使用）
gdf_dam['dist'] = gdf_dam.apply(
    lambda rec: grs80.inv(rec['geometry'].x, rec['geometry'].y,
                          rec['pref_office_point'].x, rec['pref_office_point'].y)[2], axis=1)

# 都道府県庁から最も遠いダムを表示
farthest_index = gdf_dam['dist'].idxmax(axis=1)
print(f'{gdf_dam.at[farthest_index, "prefecture"]}庁から最も遠いダム={gdf_dam.at[farthest_index, "W01_001"]}ダム\
      距離={gdf_dam.at[farthest_index, "dist"]/1000}km')

東京都庁から最も遠いダム=乳房ダム  距離=1029.7692879063898km

処理4：最も近いポイントを探す

最も近い点を探す問題は「最近傍探索」と呼ばれており、様々なライブラリが存在します。
今回はこちらの記事と同様に、Scikit-LearnのNearestNeighborsクラスを利用します。

最も近いダムを探索

from sklearn.neighbors import NearestNeighbors
# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_dam_utm = gdf_dam.copy()
gdf_dam_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_dam_utm = gdf_dam_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 全点の位置関係を学習
dam_points_array = np.array([[point.x, point.y] for point in gdf_dam_utm['geometry']])  # ndarray化
nn = NearestNeighbors(algorithm='ball_tree')
nn.fit(dam_points_array)

# ダムデータを1点ずつ走査
for i, row in gdf_dam_utm.iterrows():
    # 最近傍点を探索
    point_utm = np.array([[row['geometry'].x, row['geometry'].y]])  # ndarrayに変換
    dists, result = nn.kneighbors(point_utm, n_neighbors=3)  # 近傍上位3点を探索
    # 見付かった最近傍点(1番目は自身なので、2番目に近い点が最近傍点)を表示
    row_nearest = gdf_dam_utm.loc[result[0][1]]  # 最近傍点のデータ取得
    dist_nearest = dists[0][1]/1000  # 最近傍点までの距離(m単位なのでkm単位に変換ん)
    print(f'{row["W01_001"]}ダムから最も近いダム: {row_nearest["W01_001"]}ダム  距離={dist_nearest}km')

駒込ダムから最も近いダム: 下湯ダム  距離=9.241586379700022km
:

最も近いポイントを一括算出

NearestNeighborsクラスのkneighborsメソッドに全点の座標を渡せば、
各ポイントから最も近いポイントを一括算出できます。

下のサンプルコードでは、最も近いダムまでの距離が最大のダムを求めています

最も近いダムまでの距離が最大のダムを探索

from sklearn.neighbors import NearestNeighbors
# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_dam_utm = gdf_dam.copy()
gdf_dam_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_dam_utm = gdf_dam_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 全点の位置関係を学習
dam_points_array = np.array([[point.x, point.y] for point in gdf_dam_utm['geometry']])  # ndarray化
nn = NearestNeighbors(algorithm='ball_tree')
nn.fit(dam_points_array)

# 最近傍点を一括探索
dists, results = nn.kneighbors(dam_points_array)
# 見付かった最近傍点(1番目は自身なので、2番目に近い点が最近傍点)と距離を保持
dist_nearest = [dist/1000 for dist in dists[:,1]]  # 最近傍点までの距離(m単位なのでkm単位に変換)

# 最も近いダムからの距離が最大のダムを表示
farthest_index = np.argmax(dist_nearest)
print(f'最近傍ダムからの距離が最大のダム={gdf_dam_utm.at[farthest_index, "W01_001"]}ダム\
    最近傍ダム={gdf_dam_utm.at[results[farthest_index][1], "W01_001"]}ダム\
    距離={dist_nearest[farthest_index]}km')

最近傍ダムからの距離が最大のダム=新堤ダム    最近傍ダム=青野大師ダム    距離=203.02523976527996km

処理6：ラインの長さ測定

ラインの長さは、GeoPandasにおいてはlengthプロパティから簡単に取得することができます。
※緯度経度座標のときは、事前にUTM座標or平面直角座標に変換する必要があります

河川の長さ測定

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_river_utm = gdf_river.copy()
gdf_river_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_river_utm = gdf_river_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 河川データを1点ずつ走査
for i, row in gdf_river_utm.iterrows():
    # 長さを算出して表示
    length = row.length
    print(f'{row["W05_004"]}  長さ={length}m')

金勝川  長さ=102.29698101602632m
:

長さを一括算出

上位ランキング等の統計値を求めたい際は、一括座標変換 → lengthプロパティで一括長さ算出が便利です。
lenghtプロパティは、個々のレコードだけでなく全レコードの長さを一括取得することもできます

下の例では、滋賀県の河川長さランキングを求めています。
※河川が分割して登録されているので、河川名でグルーピングして長さの合計値算出しています
※琵琶湖と名称不明はランキング対象として不適なので、除外しています

河川の長さ上位5位を算出

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_river_utm = gdf_river.copy()
gdf_river_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_river_utm = gdf_river_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 長さを一括算出(km単位にするため1000で割る)
gdf_river_utm['river_length'] = gdf_river_utm.length / 1000

# 河川名グルーピングして長さ降順でソート
df_length = gdf_river_utm[['W05_004', 'river_length']].groupby('W05_004').sum(
    ).sort_values('river_length', ascending=False)
print(df_length.head(5))
# 名称不明と琵琶湖を除外
df_length = df_length[~df_length.index.isin(['名称不明', '琵琶湖'])]
print(df_length.head(5))

         river_length
W05_004              
野洲川         67.495013
安曇川         55.069210
愛知川         53.348960
日野川         49.800100
高時川         47.912388

処理8：ポリゴンの重心測定

ポリゴンデータの座標をどこか1つの代表点にまとめたいとき、重心算出が便利です。

GeoPandasにおいてはcentroid.coordsプロパティから重心を取得できます。

緯度経度座標の場合、精度を重視するのであれば事前に平面座標に変換する事が望ましいですが、
ライン長さとは異なり、緯度経度座標のままでもそこそこの精度で重心を求めることができます。

緯度経度座標のまま重心算出

緯度経度座標のまま重心を算出してみます

緯度経度座標のまま湖沼の重心算出

# 湖沼データを1点ずつ走査
for i, row in gdf_lake.iterrows():
    # 重心を算出
    center = list(row['geometry'].centroid.coords)[0]
    print(f'{row["W09_001"]}  重心={center}')

さっぽろ湖  重心=(141.14867614403863, 43.00018728073415)
:
霞ヶ浦  重心=(140.37323343888406, 36.042847739825625)
:
琵琶湖  重心=(136.08878954009663, 35.28527410588072)
:

こちらの記事と比べると、短い行数で重心算出できることが分かります。

UTM座標に変換して重心算出

精度向上のため、UTM座標に変換して重心を算出してみます

UTM座標に変換して湖沼の重心算出

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 湖沼データを1点ずつ走査
for i, row in gdf_lake_utm.iterrows():
    # 重心を算出
    center_utm = row['geometry'].centroid
    # 緯度経度座標に戻す（pyprojを使用。TransformPointは緯度→経度の順で返すので、元の座標系に合わせ経度を先に反転させる）
    transformer = pyproj.Transformer.from_crs('EPSG:3099', 'EPSG:4612')
    center = transformer.transform(center_utm.x, center_utm.y)[1::-1]
    print(f'{row["W09_001"]}  重心={center}')

さっぽろ湖  重心=(141.14867860778017, 43.00018596200766)
:
霞ヶ浦  重心=(140.37327735295187, 36.042836688005195)
:
琵琶湖  重心=(136.08860483036526, 35.28517147865989)
:

両算出法の差

UTM座標に変換しない場合と変換した場合の差は、こちらを参照ください

重心位置の一括算出

位置関係をランク付けしたい際などには、重心位置を一括算出すると効率的です。
centroid.coordsプロパティは、個々のレコードだけでなく全レコードの重心を一括取得することもできます

緯度経度のまま、UTM座標に変換した2パターンで、湖沼の重心位置を一括算出してみます

緯度経度座標のまま一括算出

緯度経度座標のまま重心一括算出し、最も北にある湖沼を抽出

# 重心位置を一括計算
centers = gdf_lake['geometry'].centroid
# 重心が最も北にある湖を表示
northest_index = np.argmax([center.y for center in centers])
print(f'重心が最も北にある湖={gdf_lake.at[northest_index, "W09_001"]}\
      北緯{centers.at[northest_index].y}度')

重心が最も北にある湖=久種湖      北緯45.431622616632914度

UTM座標に変換して一括算出

UTM座標に変換して高精度に重心を算出する場合、算出した重心を最終的に緯度経度座標に戻す必要があります。
このとき、set_geometryメソッドを使ってジオメトリをポリゴン→重心位置のポイントに変更し、to_crsメソッドで一括座標変換すると便利です

UTM座標に変換して重心一括算出し、最も北にある湖沼を抽出

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換
# 重心位置を一括計算
gdf_lake_utm['centroid_column'] = gdf_lake_utm['geometry'].centroid
# 重心位置をジオメトリに設定
gdf_lake_center = gdf_lake_utm.set_geometry('centroid_column')
# 緯度経度座標に戻す
gdf_lake_center = gdf_lake_center.to_crs(epsg=4612)
# 重心が最も北にある湖を表示
northest_index = gdf_lake_center['centroid_column'].y.idxmax()
print(f'重心が最も北にある湖={gdf_lake_center.at[northest_index, "W09_001"]}\
      北緯{gdf_lake_center.at[northest_index, "centroid_column"].y}度')

重心が最も北にある湖=久種湖      北緯45.4316225877112度

処理9：ポリゴンの面積測定

ポリゴンの面積は、GeoPandasにおいてはareaプロパティから簡単に取得することができます。
緯度経度座標系のときは事前に座標変換が必要なので、ご注意ください

湖沼の面積測定

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換

# 湖沼データを1点ずつ走査
for i, row in gdf_lake_utm.iterrows():
    # 面積を算出（m2 → km2に単位変換）
    area = row['geometry'].area/1000000
    print(f'{row["W09_001"]}  面積={area}km2')

さっぽろ湖  面積=1.7843064500985144km2
:
霞ヶ浦  面積=168.46126253416477km2
:
琵琶湖  面積=669.4501123974004km2
:

こちらの記事と比べると、短い行数で面積が取得できることが分かります。

ポリゴン面積を一括算出

上位ランキング等の統計値を求めたい際は、一括座標変換 → 一括面積算出が便利です。

下の例では、面積最大の湖を求めています。

面積最大の湖沼を算出

# 座標変換（緯度経度(EPSG4612) → UTM座標53N系(EPSG3099)）
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換
# 面積を一括計算
areas = gdf_lake_utm['geometry'].area/1000000
# 面積最大の湖を表示
biggest_index = np.argmax(areas)
print(f'面積最大の湖={gdf_lake_utm.at[biggest_index, "W09_001"]}  面積={areas[biggest_index]}km2')

面積最大の湖=琵琶湖  面積=669.4501123974004km2

処理10：名称・住所から緯度経度検索（ジオコーディング）

名称や住所から緯度経度を検索することを、「ジオコーディング」と呼びます。
GeoPandasにおいてジオコーデイングを行いたい際は、geopandas.tools.geocodeメソッドを利用します。

geocodeメソッドは、内部的にはジオコーディングを提供する外部サービスのAPIを叩いており、'provider'引数でサービスを選択することができます。
各APIサービスに関する詳細は、こちらをご覧ください（リクエスト数の上限にご注意ください）

下のサンプルコードでは、堤高150m以上のダムの緯度経度をNominatimのジオコーディングで検索し、Shapefileの座標と比較しています

堤高150m以上のダムの緯度経度をジオコーディングで検索

import geopandas.tools as gts
# 堤高150m以上のダムに絞る
gdf_dam_over150m = gdf_dam[gdf_dam['W01_007'] > 150]
# ジオコーディング実行
gdf_dam_geo = gts.geocode(gdf_dam_over150m['W01_001'].apply(lambda x: f'{x}ダム'),
                          provider='nominatim', user_agent='test')
# ジオコーディング結果と元のポイントを比較
gdf_dam_geo.insert(0, 'geometry_shp', gdf_dam_over150m['geometry'])  # 元のポイントを結合
print(gdf_dam_geo)

結果

                    geometry_shp                    geometry  \
587   POINT (139.07974 36.88256)  POINT (139.07863 36.88135)   
698   POINT (139.24936 35.54239)  POINT (139.24958 35.54189)   
709   POINT (139.05186 35.95399)  POINT (139.05290 35.95361)   
820   POINT (139.24977 37.15398)  POINT (139.24989 37.15436)   
833   POINT (137.66209 36.56647)  POINT (137.66364 36.56669)   
948   POINT (136.64386 36.26528)  POINT (136.64494 36.26498)   
1045  POINT (137.71840 36.13266)  POINT (137.72042 36.13263)   
1057  POINT (137.69003 36.47381)  POINT (139.24958 35.54189)   
1129  POINT (136.50203 35.66364)  POINT (136.50140 35.66599)   
1200  POINT (137.79410 35.09914)    GEOMETRYCOLLECTION EMPTY   
1211  POINT (137.79403 35.09916)    GEOMETRYCOLLECTION EMPTY   
1958  POINT (132.29948 34.63432)  POINT (132.30202 34.63473)   

                                          address  
587                    奈良俣ダム, みなかみ町, 利根郡, 群馬県, 日本  
698                     宮ケ瀬ダム, 愛川町, 愛甲郡, 神奈川県, 日本  
709                            浦山ダム, 秩父市, 埼玉県, 日本  
820                     奥只見ダム, 魚沼市, 南会津郡, 新潟県, 日本  
833                      黒部ダム, 立山町, 中新川郡, 富山県, 日本  
948                           手取川ダム, 白山市, 石川県, 日本  
1045                          奈川渡ダム, 松本市, 長野県, 日本  
1057                    宮ケ瀬ダム, 愛川町, 愛甲郡, 神奈川県, 日本  
1129                     徳山ダム, 揖斐川町, 揖斐郡, 岐阜県, 日本  
1200                                         None  
1211                                         None  
1958  中国地方整備局温井ダム管理所, 国道186号, 安芸太田町, 山県郡, 広島県, 日本

処理11：緯度経度から住所を取得（逆ジオコーディング）

ジオコーディングとは逆に、緯度経度から住所を取得する方法を「逆ジオコーディング」と呼びます。
GeoPandasにおいて逆ジオコーデイングを行いたい際は、geopandas.tools.geocodeメソッドを利用します。

堤高150m以上のダムの住所を逆ジオコーディングで検索

import geopandas.tools as gts
# 堤高150m以上のダムに絞る
gdf_dam_over150m = gdf_dam[gdf_dam['W01_007'] > 150]
# 逆ジオコーディング実行
gdf_dam_rgeo = gts.reverse_geocode(gdf_dam_over150m['geometry'],
                                   provider='nominatim', user_agent='test')
gdf_dam_rgeo.insert(0, 'W01_001', gdf_dam_over150m['W01_001'])  # ダム名を結合
print(gdf_dam_rgeo)

結果

     W01_001                    geometry  \
587      奈良俣  POINT (139.08041 36.88052)   
698      宮ヶ瀬  POINT (139.25002 35.54164)   
709       浦山  POINT (139.05178 35.95367)   
820      奥只見  POINT (139.24979 37.15367)   
833       黒部  POINT (137.66364 36.56669)   
948      手取川  POINT (136.64380 36.26518)   
1045     奈川渡  POINT (137.71827 36.13263)   
1057      高瀬  POINT (137.69014 36.47397)   
1129      徳山  POINT (136.49965 35.66185)   
1200  佐久間（再）  POINT (137.79420 35.09957)   
1211  佐久間（元）  POINT (137.79413 35.09958)   
1958      温井  POINT (132.29960 34.63417)   

                                                address  
587                       みなかみ町, 利根郡, 群馬県, 379-1721, 日本  
698   ダム下 山麓駅, 玉すだれの滝トンネル, 半原, 緑区, 愛川町, 愛甲郡, 神奈川県, 2...  
709                                        秩父市, 埼玉県, 日本  
820                                  魚沼市, 南会津郡, 新潟県, 日本  
833                     黒部ダム, かんぱ谷橋, 立山町, 中新川郡, 富山県, 日本  
948                              白山市, 石川県, 920-2336, 日本  
1045                               国道158号, 松本市, 長野県, 日本  
1057                                       大町市, 長野県, 日本  
1129                         国道417号, 揖斐川町, 揖斐郡, 岐阜県, 日本  
1200                  飯田富山佐久間線, 天竜区, 浜松市, 北設楽郡, 静岡県, 日本  
1211                  飯田富山佐久間線, 天竜区, 浜松市, 北設楽郡, 静岡県, 日本  
1958                                安芸太田町, 山県郡, 広島県, 日本

保存1：ポイントのファイル出力

Shapefile形式での保存と、GeoJSON形式での保存の2種類を紹介します。

こちらの記事では保存形式ごとにライブラリを使い分けていましたが、GeoPandasでは1つのライブラリでどちらの形式も出力できます。

ファイル出力には、to_fileメソッドを使用します。

サンプルとして、堤高100m以上のダムに絞ったデータを、フィールドを絞って出力します

Shapefile形式での出力

堤高100m以上のダムをShapefile形式で保存

# 出力用のデータ（堤高100m以上のダム）作成
gdf_dam_over100m = gdf_dam[gdf_dam['W01_007'] > 100]

# フィールドをダム名、堤高、総貯水量に絞る
gdf_dam_over100m = gdf_dam_over100m[['W01_001', 'W01_007', 'W01_010', 'geometry']]
# フィールド名を変更
gdf_dam_over100m = gdf_dam_over100m.rename(columns={'W01_001': 'ダム名',
                                                    'W01_007': '堤高',
                                                    'W01_010': '総貯水量'})

# Shapefileを出力
outpath = './dams_over100m/dams_over100m.shp'
gdf_dam_over100m.to_file(outpath, encoding='cp932')

GeoJSON形式での出力

GeoJSON形式での出力には、'driver'引数に'GeoJSON'を指定します

堤高100m以上のダムをGeoJSON形式で保存

# 出力用のデータ（堤高100m以上のダム）作成
gdf_dam_over100m = gdf_dam[gdf_dam['W01_007'] > 100]

# フィールドをダム名、堤高、総貯水量に絞る
gdf_dam_over100m = gdf_dam_over100m[['W01_001', 'W01_007', 'W01_010', 'geometry']]
# フィールド名を変更
gdf_dam_over100m = gdf_dam_over100m.rename(columns={'W01_001': 'ダム名',
                                                    'W01_007': '堤高',
                                                    'W01_010': '総貯水量'})

# GeoJSONを出力
outpath = './dams_over100m.geojson'
gdf_dam_over100m.to_file(outpath, driver='GeoJSON', encoding='cp932')

出力ファイルは以下のGitHubリポジトリにもアップロードしました
（GeoJSONはGitHubで表示できるのが素晴らしいですね！）
https://github.com/c60evaporator/shapefile_tutorials/blob/master/dams_over100m.geojson

保存2：ラインのファイル出力

Shapefile形式での保存と、GeoJSON形式での保存の2種類を紹介します。
ポイントと同様、ファイル出力にはto_fileメソッドを使用します。

サンプルとして、長さ上位5位までの河川データを、フィールドを絞って出力します

Shapefile形式での出力

長さ上位5位の河川をShapefile形式で保存

# 出力用のデータ（上位5位の河川に絞る）作成
gdf_river_top5 = gdf_river[gdf_river['W05_004'].isin(['野洲川', '安曇川', '愛知川', '日野川', '高時川'])]
# フィールドを河川名、河川コードに絞る
gdf_river_top5 = gdf_river_top5[['W05_004', 'W05_002', 'geometry']]
# フィールド名を変更
gdf_river_top5 = gdf_river_top5.rename(columns={'W05_004': '河川名',
                                                    'W05_002': '河川コード'})

# Shapefileを出力
outpath = './river_top5/river_top5.shp'
gdf_river_top5.to_file(outpath, encoding='cp932')

GeoJSON形式での出力

GeoJSON形式での出力には、'driver'引数に'GeoJSON'を指定します

長さ上位5位の河川をGeoJSON形式で保存

# 出力用のデータ（上位5位の河川に絞る）作成
gdf_river_top5 = gdf_river[gdf_river['W05_004'].isin(['野洲川', '安曇川', '愛知川', '日野川', '高時川'])]
# フィールドを河川名、河川コードに絞る
gdf_river_top5 = gdf_river_top5[['W05_004', 'W05_002', 'geometry']]
# フィールド名を変更
gdf_river_top5 = gdf_river_top5.rename(columns={'W05_004': '河川名',
                                                    'W05_002': '河川コード'})

# GeoJSONを出力
outpath = './river_top5.geojson'
gdf_river_top5.to_file(outpath, driver='GeoJSON', encoding='cp932')

出力したファイルは以下のGitHubリポジトリにアップロードしております
https://github.com/c60evaporator/shapefile_tutorials/blob/master/river_top5.geojson

保存3：ポリゴンのファイル出力

Shapefile形式での保存と、GeoJSON形式での保存の2種類を紹介します。
ポイントと同様、ファイル出力にはto_fileメソッドを使用します。

サンプルとして、面積100km2以上の湖沼データを、フィールドを絞って出力します

Shapefile形式での出力

面積100km2以上の湖沼データをShapefile形式で保存

# 出力用のデータ（面積100km2以上の湖沼）作成
# UTM座標変換
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換
# 面積100km2以上のデータ抽出
gdf_lake_utm['lake_area'] = gdf_lake_utm.area / 1000000
gdf_lake_over100km2 = gdf_lake_utm[gdf_lake_utm['lake_area'] > 100]
# 座標を緯度経度に戻す
gdf_lake_over100km2 = gdf_lake_over100km2.to_crs(epsg=4612)  # 変換後座標に変換
# フィールドを湖沼名、最大水深、面積に絞る
gdf_lake_over100km2 = gdf_lake_over100km2[['W09_001', 'W09_003', 'lake_area', 'geometry']]
# フィールド名を変更
gdf_lake_over100km2 = gdf_lake_over100km2.rename(columns={'W09_001': '湖沼名',
                                                'W09_003': '最大水深',
                                                'lake_area': '面積'})

# Shapefileを出力
outpath = './lake_over100km2/lake_over100km2.shp'
gdf_lake_over100km2.to_file(outpath, encoding='cp932')

こちらの記事と比べると、コードが短く済むことが分かりますね！

GeoJSON形式での出力

GeoJSON形式での出力には、'driver'引数に'GeoJSON'を指定します

面積100km2以上の湖沼データをGeoJSON形式で保存

# 出力用のデータ（面積100km2以上の湖沼）作成
# UTM座標変換
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換
# 面積100km2以上のデータ抽出
gdf_lake_utm['lake_area'] = gdf_lake_utm.area / 1000000
gdf_lake_over100km2 = gdf_lake_utm[gdf_lake_utm['lake_area'] > 100]
# 座標を緯度経度に戻す
gdf_lake_over100km2 = gdf_lake_over100km2.to_crs(epsg=4612)  # 変換後座標に変換
# フィールドを湖沼名、最大水深、面積に絞る
gdf_lake_over100km2 = gdf_lake_over100km2[['W09_001', 'W09_003', 'lake_area', 'geometry']]
# フィールド名を変更
gdf_lake_over100km2 = gdf_lake_over100km2.rename(columns={'W09_001': '湖沼名',
                                                'W09_003': '最大水深',
                                                'lake_area': '面積'})
# GeoJSONを出力
outpath = './lake_over100km2.geojson'
gdf_lake_over100km2.to_file(outpath, driver='GeoJSON', encoding='cp932')

出力したファイルは以下のGitHubリポジトリにアップロードしております
https://github.com/c60evaporator/shapefile_tutorials/blob/master/lake_over100km2.geojson

表示1：ポイントの表示

GeoPandasには、Matplotlibを使ってジオメトリを表示する、plotメソッドが存在します。
実装が非常にシンプルで便利ですが、インタラクティブに地図を表示したい場合はfoliumライブラリを使用するのがお薦めです

下の例では、plotメソッドを使って堤高100m以上のダムを表示しています

堤高100m以上のダムを表示（ポイントのみ）

# 表示用のデータ（堤高100m以上のダム）作成
gdf_dam_over100m = gdf_dam[gdf_dam['W01_007'] > 100]
# ポイントをプロット
gdf_dam_over100m.plot(column = 'W01_007',  # 色分け対象の列
                      cmap = 'OrRd'  # 色分けのカラーマップ
                      )

日本地図上でポイント表示

ポイントだけだと位置関係が分かりづらいので、日本地図上に表示ができると便利です。
日本地図は、**japanmapライブラリ**を使って読み込んだポリゴンをGeoDataFrame形式に変換すると、plotメソッドで簡単に表示できてオススメです。

下の例では、堤高100m以上のダムを堤高によって色分けして日本地図上に表示しています
（色分け対象列には'column'引数を、カラーマップには'cmap'引数を、カラーバー表示には'legend'引数を使用します）

堤高100m以上のダムを表示（日本地図上に表示）

from japanmap import get_data, pref_points
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
# 表示用のfigure作成
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
# 日本地図をプロット
gdf_pref.plot(ax = ax,
              color = 'gray'  # 塗りつぶし色を指定
              )

# 表示用のデータ（堤高100m以上のダム）作成
gdf_dam_over100m = gdf_dam[gdf_dam['W01_007'] > 100]

# ポイントをプロット
gdf_dam_over100m.plot(ax = ax,  # 描画先のax
                      column = 'W01_007',  # 色分け対象の列
                      cmap = 'OrRd',  # 色分けのカラーマップ
                      legend = True,  # 色分けのカラーバー表示
                      legend_kwds = {'label': 'dam height',  # カラーバーにラベル設定
                                     'shrink': 0.6},  # カラーバーが長すぎるので短く
                      s = 6  # 点マーカーのサイズ
                      )

表示2：ラインの表示

ポイントデータのときとほぼ同様の手順です。

下のサンプルでは、滋賀県の長さ上位5河川を日本地図上で表示しています
（日本地図は滋賀県のみでフィルタリングしています）

滋賀県の長さ上位5河川を表示（日本地図上に表示）

from japanmap import pref_names, get_data, pref_points
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
# 表示用のfigure作成
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
gdf_pref['prefecture'] = pref_names[1:]  # 県名を格納
# 滋賀県に絞る
gdf_pref = gdf_pref[gdf_pref['prefecture'] == '滋賀県']
# 日本地図をプロット
gdf_pref.plot(ax = ax,
              color = 'gray'  # 塗りつぶし色を指定
              )

# 表示用のデータ（上位5位の河川に絞る）作成
gdf_river_top5 = gdf_river[gdf_river['W05_004'].isin(['野洲川', '安曇川', '愛知川', '日野川', '高時川'])]
# フィールドを河川名、河川コードに絞る
gdf_river_top5 = gdf_river_top5[['W05_004', 'W05_002', 'geometry']]
# フィールド名を変更
gdf_river_top5 = gdf_river_top5.rename(columns={'W05_004': '河川名',
                                                    'W05_002': '河川コード'})

# ラインをプロット
gdf_river_top5.plot(ax = ax,  # 描画先のax
                    color = 'blue'
                    )

表示3：ポリゴンの表示

ポイントデータのときとほぼ同様の手順です。

下のサンプルでは、面積50km2以上の湖を最大水深で色分けし、日本地図上で表示しています
（日本地図は県コード30以下（関西以東）でフィルタリングしています）

面積50km2以上の湖沼を表示（日本地図上に表示）

from japanmap import pref_names, get_data, pref_code
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
# 表示用のfigure作成
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
gdf_pref['prefecture'] = pref_names[1:]  # 県名を格納
gdf_pref['pref_code'] = gdf_pref['prefecture'].apply(lambda x: pref_code(x))  # 県コードを格納
# 県コード30以下に絞る
gdf_pref = gdf_pref[gdf_pref['pref_code'] <= 30]
# 日本地図をプロット
gdf_pref.plot(ax = ax,
              color = 'gray'  # 塗りつぶし色を指定
              )

# UTM座標変換
gdf_lake_utm = gdf_lake.copy()
gdf_lake_utm.crs = f'epsg:{4612}'  # 変換前座標を指定
gdf_lake_utm = gdf_lake_utm.to_crs(epsg=3099)  # 変換後座標に変換
# 面積100km2以上のデータ抽出
gdf_lake_utm['lake_area'] = gdf_lake_utm.area / 1000000
gdf_lake_over100km2 = gdf_lake_utm[gdf_lake_utm['lake_area'] > 100]
# 座標を緯度経度に戻す
gdf_lake_over100km2 = gdf_lake_over100km2.to_crs(epsg=4612)  # 変換後座標に変換
# フィールドを湖沼名、最大水深、面積に絞る
gdf_lake_over100km2 = gdf_lake_over100km2[['W09_001', 'W09_003', 'lake_area', 'geometry']]
# フィールド名を変更
gdf_lake_over100km2 = gdf_lake_over100km2.rename(columns={'W09_001': '湖沼名',
                                                'W09_003': '最大水深',
                                                'lake_area': '面積'})
# 色分け用の最大水深フィールドを文字型→数値型に修正
gdf_lake_over100km2 = gdf_lake_over100km2.astype({'最大水深': float})

# ラインをプロット
gdf_lake_over100km2.plot(ax = ax,  # 描画先のax
                         column = '最大水深',  # 色分け対象の列
                         cmap = 'Blues',  # 色分けのカラーマップ
                         legend = True,  # 色分けのカラーバー表示
                         legend_kwds = {'label': 'depth',  # カラーバーにラベル設定
                                     'shrink': 0.6},  # カラーバーが長すぎるので短く
                         )