半径2kmのジオフェンスで地物抽出します

• ライブラリのインポート
• データの読み込み
• 駅周辺を抽出(半径1km)
  • データの読み込み
  • 2kmのジオフェンス(円)を作成
• ジオフェンス内のデータを抽出
  • keplerで描画
• まとめ

ライブラリのインポート

import base64
import pandas as pd
import numpy as np

from shapely import wkt
from shapely.geometry import Point

try:
  import geopandas as gpd
  import leafmap.kepler as leafmap
  from keplergl import KeplerGl
except ImportError:
  !pip install geopandas
  import geopandas as gpd
  !pip install keplergl
  !pip install git+https://github.com/giswqs/leafmap.git
  import leafmap.kepler as leafmap
  from keplergl import KeplerGl

データの読み込み

国勢調査 (opens new window)からデータをダウンロードして，GeoPandasで読み込む．

import glob

paths = glob.glob("./data/*/*.shp")

gdf_pop = pd.DataFrame()
for path in paths:
  # 駅データの読み込み
  gdf_pop_wk = gpd.read_file(path)
  gdf_pop = pd.concat([gdf_pop, gdf_pop_wk],axis=0)
  gdf_pop.head(2)

# 人口密度
gdf_pop["PDE"] = gdf_pop["JINKO"]/gdf_pop["AREA"] 

駅周辺を抽出(半径1km)

下記のサイトから、鉄道データをダウンロードします。

そして、下記のコードで読み込みます。 駅拠点 (opens new window)

データの読み込み

# 新幹線駅の読み込み
df_station_geo = pd.read_csv("./data/station.csv")
# geometryに変換
df_station_geo["geometry"] = df_station_geo["geometry"].apply(wkt.loads)
df_station_geo.head(2)

2kmのジオフェンス(円)を作成

def circle(center, meter):
  local_azimuthal_projection = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(
      lat, lon
  )
  wgs84_to_aeqd = partial(
      pyproj.transform,
      pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
      pyproj.Proj(local_azimuthal_projection),
  )
  aeqd_to_wgs84 = partial(
      pyproj.transform,
      pyproj.Proj(local_azimuthal_projection),
      pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
  )

  point_transformed = transform(wgs84_to_aeqd, center)
  buffer = point_transformed.buffer(meter)
  # Get the polygon with lat lon coordinates
  circle_poly = transform(aeqd_to_wgs84, buffer)
  return circle_poly

df_station_geo["geometry"] = df_station_geo["geometry"].apply(lambda x: circle(x, 2000))
df_station_geo.head(2)

ジオフェンス内のデータを抽出

df_plot = pd.DataFrame()
for i in range(len(df_station_geo)):
  # 境界データ
  bound = df_station_geo["geometry"][i]

  # ジオフェンスに含まれるか
  geometry = gdf_pop["geometry"].centroid
  _bool = geometry.within(bound)
  tmp_gdf = gdf_pop[_bool]
  # plot用のデータフレーム
  df_plot = pd.concat([df_plot, tmp_gdf],axis=0)

# 重複削除
df_plot = df_plot.drop_duplicates()
df_plot.head()

keplerで描画

# データを描画
df_station_geo_wk = df_station_geo.copy()
df_plot_wk = df_plot.copy()
df_station_geo_wk["geometry"] = df_station_geo_wk["geometry"].astype(str)
# df_plot_wk["geometry"] = df_plot_wk["geometry"].astype(str)
map1 = KeplerGl(height=600)
map1.add_data(data=df_plot_wk, name='df_plot')
map1.add_data(data=df_station_geo_wk, name='df_station_geo')
map1

まとめ

半径2kmのジオフェンスで地物抽出しました。