Pythonで特定の地物と重なる地物を抽出する
/
Pythonで特定の地物と重なる地物を抽出します.
※ただし、時間がかかる
# ファイル構成
project_dir
├── /data
│ ├── /R2_boundary/04_miyagi/04_miyagi.shp
│ ├── /osm/tohoku/gis_osm_roads_free_1.shp
│ └── (省略)
└── merge_bound.ipynb <- 実行用ノートブック
# shpファイルを読み込み
shpファイルを読み込みます.
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.ops import unary_union
# osmデータのshpファイル読込
osm_gdf = gpd.read_file("./data/osm/tohoku/gis_osm_roads_free_1.shp")
display(osm_gdf.head())
# 境界データのshpファイル読込
pref_gdf = gpd.read_file("./data/R2_boundary/04_miyagi/04_miyagi.shp")
# display(pref_gdf.head(2))
# fig = plt.figure()
# ax = fig.add_subplot()
# 境界データをmerge
bound = gpd.GeoSeries(unary_union(pref_gdf.geometry))
bound.crs = "epsg:4326"
# # 描画
# bound.plot(ax=ax, color='none', edgecolor='b', linewidth=1)
# plt.show()
# 地物の抽出(within)
import time
start = time.time()
# 一時データ
osm_gdf_0_10 = osm_gdf[0:10]
# 県の境界データに含まれるか
# LINESTRINGの重心で判定
geometry = osm_gdf_0_10["geometry"]
_bool = (geometry.centroid).within(bound.iloc[0])
tmp_gdf = osm_gdf_0_10[_bool]
elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")
# elapsed_time:0.4750847816467285[sec]
# 地物の抽出(contains)
import time
start = time.time()
# 一時データ
osm_gdf_0_10 = osm_gdf[0:10]
# 県の境界データに含まれるか
# LINESTRINGの重心で判定
geometry = osm_gdf_0_10["geometry"]
_bool = geometry.apply(lambda x: bound.contains(x.centroid))
tmp_gdf = osm_gdf_0_10[_bool]
elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")
# elapsed_time:0.5308980941772461[sec]
# まとめ
Pythonで特定の地物と重なる地物を抽出しました.
# 参考サイト
Geopandas - ValueError: Cannot transform naive geometries. Please set a crs on the object first (opens new window)
【PythonでGIS】GeoPandasまとめ (opens new window)
Geopandas contains working for one point but not many (opens new window)
Get min/max lat and long values from geodataframe (opens new window)