携帯基地局のカバー面積を集計

!pip install geopandas
!pip install rtree
!pip install pygeos

!wget https://www.esrij.com/cgi-bin/wp/wp-content/uploads/2012/10/japan_ver83.zip

url = "https://docs.google.com/spreadsheets/d/e/2PACX-1vTHEggpuuf5ratyC4d3PFpsYaw46tVopeFglSIUUb47TcEpCOQmjv4OapFDTfBrLqbwDByiJ3Nc4DZ9/pub?gid=0&single=true&output=csv"

import pandas as pd
import geopandas as gpd

# 地理座標系

japan = gpd.read_file("japan_ver83.zip!japan_ver83")
japan

nara = japan[japan["KEN"] == "奈良県"]
nara

nara.crs

nara.plot()

# 基地局データ

df = pd.read_csv(url, index_col=0)

df

pt_df = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df.lng, df.lat), crs="EPSG:6668")

base = nara.plot(color="white", edgecolor="black")
pt_df.plot(ax=base, marker="o", color="red", markersize=5)

# 平面直角座標系

# https://www.gsi.go.jp/sokuchikijun/jpc.html
# https://homata.gitbook.io/geodjango/hajimeteno/coordinate


# 地域によって変更
nara.to_crs(epsg=6674, inplace=True)
pt_df.to_crs(epsg=6674, inplace=True)

pt_df.crs

nara.crs

circles = pt_df.geometry.buffer(710)
circles

circles.plot(cmap="tab20b")

mp = circles.unary_union
mp

nara_map = nara["geometry"].difference(mp)

nara.area

nara_map.area

nara_map.plot()

df1 = pd.concat([nara["SIKUCHOSON"], nara.area.round(3), nara_map.area.round(3)], axis=1)

df1.to_csv("nara.tsv", header=False, sep="\t")