GISで空間分析をしてみる－近接性の測定

import glob
import geopandas as gpd
import pandas as pd

# 読み込みと結合
region = pd.concat([gpd.read_file(file, encoding='cp932') for file in glob.glob('input\*.shp')])

# CRSの変換
region = region.to_crs('EPSG:6677')

# 保存
region.to_file('01region.shp', encoding='cp932')

import geopandas as gpd

# 小地域の読み込み
region = gpd.read_file('01region.shp', encoding='cp932')

# 面積計算（小地域にもあるが、按分計算用に改めて計算）
region['area_all'] = region.area

# 用地の読み込み
site = gpd.read_file('02site.shp', encoding='cp932')

# バッファ作成とオーバーレイ
for r in range(1000, 3001, 1000):
    # バッファ作成
    buffer = gpd.GeoDataFrame({'distance': [r]}, geometry=site.buffer(r).geometry)
    buffer.to_file('03buffer-{}.shp'.format(r))

    # オーバーレイ
    zone = gpd.overlay(region, buffer, how='intersection')
    
    # 人口の按分計算
    zone['area_cut'] = zone.area # バッファにより刻まれたポリゴンの面積
    zone['pop'] = zone['JINKO'] * (zone['area_cut'] / zone['area_all']) # 按分計算
    zone.to_file('04zone-{}.shp'.format(r), encoding='cp932')

import geopandas as gpd
import pandas as pd

# ズームレベル(z)
z = 16

# タイル座標（x,y）
u = (58181, 25816) # 左上
d = (58197, 25831) # 右下

# 読み込み
pieces = []
for x in range(u[0], d[0]+1):
    for y in range(u[1], d[1]+1):
        try:
            gdf = gpd.read_file('https://cyberjapandata.gsi.go.jp/xyz/experimental_rdcl/{}/{}/{}.geojson'.format(z, x, y))
            pieces.append(gdf)
        except:
            pass

# 結合
road = pd.concat(pieces) 

# CRSの変換
road = road.to_crs('EPSG:6677')

# 保存
road.to_file('05road.shp', encoding='cp932')

import geopandas as gpd

# 小地域の読み込み
region = gpd.read_file('01region.shp', encoding='cp932')

# 面積計算（小地域にもあるが、按分計算用に改めて計算）
region['area_all'] = region.area

# 道路圏域の読み込み
r_buffer = gpd.read_file('06road-buffer.shp', encoding='cp932')

# 道路圏域とオーバーレイ
for r in [1000, 2000]:
    # 距離別の道路圏域を抽出
    buffer = r_buffer[r_buffer['cost_level']==r]

    # オーバーレイ
    zone = gpd.overlay(region, buffer, how='intersection')
    
    # 人口の按分計算
    zone['area_cut'] = zone.area # バッファにより刻まれたポリゴンの面積
    zone['pop'] = zone['JINKO'] * (zone['area_cut'] / zone['area_all']) # 按分計算
    zone.to_file('07road-zone-{}.shp'.format(r), encoding='cp932')

import requests
import geopandas as gpd
import json
import time

# 3キロ圏域の小地域データを読み込み
region = gpd.read_file('04zone-3000.shp', encoding='cp932')

# 中心点（重心）を作成
center = gpd.GeoDataFrame(region.copy(), geometry=region.centroid)

# ナビタイムでは緯度経度を使用するのでCRSを変換
center = center.to_crs('EPSG:4612')

# ナビタイムAPIの設定
url = 'https://navitime-route-car.p.rapidapi.com/route_car'
headers = {'x-rapidapi-key': "xxxxxxxxxx",
           'x-rapidapi-host': "navitime-route-car.p.rapidapi.com"}

pieces = []
for i, lon , lat in zip(center.index, center['geometry'].x, center['geometry'].y):
    # 始点・終点の設定
    st = '{},{}'.format(lat, lon)
    gr = '35.60700155540336,139.6399715455508'

    # ルート検索方法
    querystring = {"start":st,"goal":gr,"datum":"wgs84","coord_unit":"degree","condition":"free_distance"} # 距離優先・有料道路つかわない
    #querystring = {"start":st,"goal":gr,"datum":"wgs84","coord_unit":"degree","condition":"free_time"} # 時間優先・有料道路つかわない

    # データ取得
    try:
        res = requests.request("GET", url, headers=headers, params=querystring)
        j = json.loads(res.text)
        distance = j['items'][0]['summary']['move']['distance'] 

        # 確認用でjsonファイルを保存しておく
        with open('save/navi{}.json'.format(i), 'w') as f:
            json.dump(j, f, indent=4)

    except:
        print('エラー', i, lat, lon)
        distance = -1

    pieces.append(distance)
    time.sleep(3) # 1分あたり50コールまでなので間隔をあける

# 距離データを追加
center['navi'] = pieces

# 検索できなかったデータを除いたポイントのシェープファイルを保存
center[center['distance']>0].to_file('08navi-point.shp', encoding='cp932')

import geopandas as gpd

# 小地域の読み込み
region = gpd.read_file('01region.shp', encoding='cp932')

# 面積計算（小地域にもあるが、按分計算用に改めて計算）
region['area_all'] = region.area

# 道路圏域の読み込み
n_buffer = gpd.read_file('09navi-buffer.shp', encoding='cp932')

# 道路圏域とオーバーレイ
for r in [1000, 2000]:
    # 距離別の道路圏域を抽出
    buffer = n_buffer[n_buffer['navi']==r]

    # オーバーレイ
    zone = gpd.overlay(region, buffer, how='intersection')
    
    # 人口の按分計算
    zone['area_cut'] = zone.area # バッファにより刻まれたポリゴンの面積
    zone['pop'] = zone['JINKO'] * (zone['area_cut'] / zone['area_all']) # 按分計算
    zone.to_file('10navi-zone-{}.shp'.format(r), encoding='cp932')