Python(psycopg2) + PostgreSQLで特定の地物と重なる地物を抽出する(POLYGON)
/
Python(psycopg2) + PostgreSQLで特定の地物と重なる地物を抽出します(POLYGON)
# ファイル構成
project_dir
├── /data
│ ├── /R2_boundary/04_miyagi/04_miyagi.shp
│ ├── /osm/tohoku/gis_osm_roads_free_1.shp
│ └── (省略)
└── merge_bound.ipynb <- 実行用ノートブック
# ライブラリのインポート
import os
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.ops import unary_union
import psycopg2
from shapely import wkt
from sqlalchemy import create_engine
from geoalchemy2 import Geometry, WKTElement
# OSMのDB作成
# osmデータのshpファイル読込
osm_gdf = gpd.read_file("./data/osm/tohoku/gis_osm_roads_free_1.shp")
# WKTに変換
osm_gdf['geometry'] = osm_gdf['geometry'].apply(lambda x: WKTElement(x.wkt, srid=4326))
# DBのURL
DATABASE_URL='postgresql://postgre:postgre@workspace-postgres-1:5432/postgres'
# テーブル作成のDB起動
engine = create_engine(DATABASE_URL)
# テーブル作成 if_exists='replace' or 'append'
osm_gdf.to_sql('osm',con=engine,if_exists='replace',index=None,
dtype={'geometry': Geometry(geometry_type='LINESTRING', srid= 4326)})
# 境界データのDB作成(最大面積のPOLYGON)
MULTIPOLYGONからPOLYGONへ分解し、最大面積のPOLYGONを抽出します。
# 境界データのshpファイル読込
pref_gdf = gpd.read_file("./data/R2_boundary/04_miyagi/04_miyagi.shp")
# 境界データをmerge
bound = gpd.GeoSeries(unary_union(pref_gdf.geometry))
bound.crs = "epsg:4326"
# MULTIPOLYGONからPOLYGONへ分解
bound_df = bound.explode().reset_index()
bound_df = bound_df.rename(columns={0: "geom"})
# 面積の算出
bound_df["area"] = bound_df["geom"].area
# 面積でソート
bound_df = bound_df.sort_values('area', ascending=False).reset_index(drop=True)
# 宮城県でデータフレーム作成
pgsl_miyagi = pd.DataFrame([["miyagi", bound_df["geom"][0]]], columns=["pref", "geometry"])
# WKTに変換
pgsl_miyagi['geometry'] = pgsl_miyagi['geometry'].apply(lambda x: WKTElement(x.wkt, srid=4326))
# DBのURL
DATABASE_URL='postgresql://postgre:postgre@workspace-postgres-1:5432/postgres'
# テーブル作成のDB起動
engine = create_engine(DATABASE_URL)
# テーブル作成 if_exists='replace' or 'append'
pgsl_miyagi.to_sql('pref_miyagi',con=engine,if_exists='replace',index=None,
dtype={'geometry': Geometry(geometry_type='POLYGON', srid= 4326)})
# POLYGON内の道路を抽出
PostGISを使う: データ管理とクエリ (opens new window)
import os
import psycopg2
import pandas as pd
# DATABASE_URL='postgresql://postgre:postgre@localhost:5432/postgres'
DATABASE_URL='postgresql://postgre:postgre@workspace-postgres-1:5432/postgres'
# postgresの接続
conn = psycopg2.connect(DATABASE_URL)
# conn.autocommit = True # 操作の重複を防ぐ(databaseの操作)呪文
# cur = conn.cursor()
# テーブル名一覧を取得するSQL
query = """
SELECT
t1.osm_id, t1.code, t1.fclass, t1.name,
ST_AsText(t1.geometry) AS geometry,
ST_Intersects(t1.geometry, t2.geometry)
FROM
public.osm AS t1,
public.pref_miyagi AS t2
WHERE
t2.pref='miyagi' AND
ST_Intersects(t1.geometry, t2.geometry);"""
df = pd.read_sql(query, con=conn)
# dbとカーソルを閉じる
conn.close()
display(df.head())
osm_id | code | fclass | name | geometry | st_intersects |
---|---|---|---|---|---|
15450378 | 5121 | unclassified | None | LINESTRING(140.6271628 37.9942567,140.627269 3... | True |
115450383 | 5121 | unclassified | None | LINESTRING(140.6130333 37.991994,140.6131692 3... | True |
115450398 | 5121 | unclassified | None | LINESTRING(140.621601 37.990825,140.6218304 37... | True |
115450693 | 5122 | residential | None | LINESTRING(140.6343659 37.994537,140.6347219 3... | True |
115450727 | 5121 | unclassified | None | LINESTRING(140.6315763 37.9931608,140.6327262 ... | True |
# shpファイルで出力
※時間がかかる
# geopandasに変換
df_shp = df.copy()
geometry = df_shp["geometry"].apply(wkt.loads)
df_shp.drop('geometry', axis=1, inplace=True)
geo_df = gpd.GeoDataFrame(df_shp, geometry=geometry)
#Shape出力
!mkdir -p ./data/shp_dir
geo_df.to_file(driver='ESRI Shapefile', filename="./data/shp_dir/miyagi_osm.shp")
# まとめ
Python(psycopg2) + PostgreSQLで特定の地物と重なる地物を抽出しました(POLYGON)
# 参考サイト
The Shapely User Manual (opens new window)
PostGISを使う: データ管理とクエリ (opens new window)
geopandasを使ってShapeファイルを作成しよう! ~Airbnbのデータを可視化してみよう~ (opens new window)
pandas-drop (opens new window)
TypeError: Input geometry column must contain valid geometry objects (opens new window)
PostGISを使う: データ管理とクエリ (opens new window)