geopandasでみえない交差点を調べる

みえない交差点

データ

オープンデータ

都道府県コード

日本測地系2011（JGD2011）

インストール

pip install pandas
pip install geopandas
pip install folium

プログラム

# 都道府県コード（神奈川県）
pref = 45

# 日本測地系2011（神奈川県）
epsg = 6677

# プログラム

import pandas as pd
import geopandas as gpd

from tqdm.notebook import tqdm

pd.options.plotting.backend = "plotly"
pd.set_option("display.max_columns", None)

df0 = pd.read_csv("https://www.npa.go.jp/publications/statistics/koutsuu/opendata/2021/honhyo_2021.csv", encoding="cp932")

# 緯度経度を度分秒から10進数表記に変換

def dms2deg(df_tmp, col):

    df = df_tmp.copy()
    
    df["d"], df["t"] = df[col].divmod(10000000)
    df["m"], df["s"] = df["t"].divmod(100000)

    result = df["d"] + (df["m"] / 60.0) + (df["s"] / 1000.0 / 3600.0 )

    return result

df0["lat"] = dms2deg(df0, "地点　緯度（北緯）")
df0["lon"] = dms2deg(df0, "地点　経度（東経）")

df0

df0.columns

df0["発生日時　　年"].value_counts()

# 2021年のみ
df1 = df0[df0["発生日時　　年"] == 2021].copy()

df1.shape

# 月ごとの事故数を確認
df1["発生日時　　月"].value_counts()

# 曜日別の事故数を確認
df1["曜日(発生年月日)"].value_counts()

# 昼夜別の事故数を確認
df1["昼夜"].value_counts()

# 道路の形ごとの事故数を確認
df1["道路線形"].value_counts()

# 交差点の大きさごとの事故数を確認
df1["車道幅員"].value_counts()

df2 = df1[df1["車道幅員"].isin([11, 14, 15]) & (df1["都道府県コード"] == pref)]
df2

geo_df = gpd.GeoDataFrame(df2, geometry = gpd.points_from_xy(df2.lon, df2.lat), crs=6668)
geo_df = geo_df.to_crs(epsg)

df3 = geo_df.copy()

# 半径10メートル
df3["buffer"] = df3.buffer(10)

# 10mの範囲
dft = df3["buffer"].apply(lambda x: df3.geometry.within(x))

dft.sum().value_counts().sort_index()

# さらに10mの範囲×2回
dfp = dft.apply(lambda x: dft.loc[dft.loc[x].any()].any(), axis=1)

dfp.sum().value_counts().sort_index()

# 件数
df3["len"] = dfp.sum()

# indexのタプル
df3["idx"] = dfp.apply(lambda x: tuple(dfp.loc[x].index.to_list()))

df3

df3["len"].value_counts()

# 5より多い事故の場所を危険交差点

df4 = df3[df3["len"] > 5].copy()

df4

se = df4.groupby("idx").apply(len).rank(ascending=False, method="first").astype(int).apply(lambda x: f"intersection_{x}")
se.name = "intersection_name"

df5 = pd.merge(df4, se, on="idx")

df5

地図表示

import folium

m = folium.Map(
    location=[df5.lat.mean(), df5.lon.mean()],
    zoom_start=12,
    tiles="https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png",
    attr='&copy; <a href="https://maps.gsi.go.jp/development/ichiran.html">国土地理院</a>',
)

for i, r in df5.iterrows():

    folium.Marker(
        location=[r.lat, r.lon],
        icon=folium.Icon(color="red")
    ).add_to(m)

m