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

最終更新日

Comments: 0

ブックオフでまれにGISの本を見かけることがあります。古本なので情報は古いものの、分析手法はまだまだ現役のものもあります。今回紹介するのは古今書院「GISを利用した社会・経済の空間分析」です。この本にある近接性の測定を現在のデータで挑戦してみました。

この本は日本大学文理学部の先生方が書かれているので、おそらく教科書として使用されていたと思います。10年以上前のものになりますが、分析課題がわかりやすく、図はカラー印刷になっていて見やすかったです。なお、GISの使用ソフトはArcGISです。(残念ながらQGISは登場しません。)

 

 

チャレンジする項目

第10章「交通条件と交通規制を考慮した近接性の測定」という項目があり、3番目に「施設からの近接性の測定」という項目があります。ここでは世田谷区のある地点から1キロ、2キロ圏内にどれだけの人口が居住しているかGISを使用して測定するものです。今回はこの本の方法を参考に、最新データを使用して再現してみようと思います。

留意点として、本ではArcGISを使っているだけあって(研究のためかスペック高め?)、統計データや電子地図データも有料のものが使用されています。これらは極力無料で手に入る物で置き換えます。また、この章に限らず他の場所でも出てくるのですが、「データはプログラム(VB)で処理した」といった箇所があるのですが、そのコードなどが公開されていないため、その処理については推測となります。

 

データ入手(国勢調査の境界データ)

GISを使ったことがあれば必ず通る?小地域データのシェープファイルをe-Statからダウンロードします。対象は東京都は大田区、世田谷区、目黒区、神奈川県は川崎市高津区、川崎市中原区です。

pythonでそれぞれのシェープファイルを結合します。QGISでもできますが、こちらの方が楽だと思います。

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')

QGISで確認してみます。 

 

対象物件の設定

本に掲載されている図が小さいのでわかりにくいのですが、世田谷区野毛3丁目あたりです。QGIS上で設定しました。タイルマップを使うと便利です。(使わないと設定できないかも)

本ではこの物件を「用地A」としているので、以降「用地」と呼ぶことにします。

 

 

分析1「直線距離圏内の人口数の測定」

用地から直線1キロ圏内と2キロ圏内の人口を測定します。本では町丁目の中心点をクリップする方法となっていますが、その後の人口計算方法がそのままなのか、面積按分しているのか記載されていません。(面積按分した方が精度があがるという注意書きはありますが)

ここでは面積按分をする方法で行います。バッファの作成はQGISでもできるのですが、その後の面積按分が面倒なので、ここもPythonを使用して行います。

※分析は2キロまでですが、あとの処理のため、3キロのバッファも作成します。

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')

できあがったシェープファイルをQGISで確認してみます。バッファにより分割された町丁目は人口が面積比で按分されていることがわかります。

あとはDBFファイルをExcelで編集して集計してみます。結果は・・・。

<本に掲載されているデータ>

<今回算出したデータ>

1キロ圏2キロ圏
面積3.1㎢12.6㎢
居住地区数1559
人口32,171140,009
人口密度10,378/㎢11,112/㎢
1キロ圏2キロ圏
面積3.1㎢12.5㎢
居住地区数2583
人口35,675166,138
人口密度11,508/㎢13,291/㎢

使用するデータの出所が違うのか随分と居住地区数が違っています。人口も少し多いような気がします。刊行当時の人口と比較するとこのくらいかも知れません。

 

データ入手(道路のデータ)

次の分析に移る前に道路のデータが必要なので用意します。どれを使おうか迷ったのですが、今回は国土地理院のベクトルタイル提供実験からダウンロードできるデータを使います。データの仕様はこちら

ダウンロードするときはPythonでまとめてできるようにしています。

①必要なタイルを調べる

こちらにタイルの番号一覧があるので、取得したい地域の左上の番号と右下の番号を控えます。道路データはズームレベルが16なので、各メッシュに表示されているコードが「16/××/××」と表示されるレベルまで地図をズームして番号を調べます。

②Pythonのコードでダウンロード&シェープファイルの出力

下のコードを使用してダウンロードします。(他に流用される方は座標系と範囲の大きさに注意)

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')

試しにQGISで表示してみます。すこし広く取得してしまいましたが問題なさそうです。

 

分析2「道路距離圏内の人口数の測定」

用地から道路距離で1キロ、2キロの圏域を作成します。本ではArcGISを使用しているので簡単にできてしまうようですが、QGISを使用するのでプラグインを使用して作成します。

こちら(QGISについての多くの解説があります)の方法を参照し、QNEAT3というプラグインを使用して到達圏域を作成します。

圏域をポリゴンデータで欲しいので、「Iso-Area as Polygons(from Point) 」で道路距離圏のレイヤを作成し、これを保存します。

分析1と同様に、今回作った道路距離圏と小地域を使って面積按分した人口データを算出します。

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')

分析1と同様にシェープファイルを見てみます。道路距離圏でオーバーレイされたシェープファイルができあがっています。

 

同じくDBFファイルをExcelで編集して集計してみます。

<本に掲載されているデータ>

<今回算出したデータ>

1キロ圏2キロ圏
面積1.9㎢
(61.3%)
9.5㎢
(75.4%)
居住地区数9
(60.0%)
46
(78.0%)
人口21,628
(67.2%)
102,428
(73.2%)
1キロ圏2キロ圏
面積1.6㎢
(51.6%)
7.5㎢
(60.0%)
居住地区数18
(72.0%)
57
(68.7%)
人口18,726
(52.5%)
93,746
(56.4%)

※()内は直線距離のデータに対する割合

道路距離圏が本よりも狭く作成されました。プラグイン(QNEAT3)の設定はデフォルトのままで行っています。もしかしたら速度などを変更した方が良かったのかも?

ですが、著しく傾向が異なるようには見えませんので良いかなぁと。ついでにシェープファイルの形はなんとなくですが似た形になっていました。

 

分析3「ナビゲーション距離の測定」

実際に本では「複雑な交通条件を考慮したナビゲーション距離の測定」という項目になっています。

分析2では道路網の条件がなかったため、短い区間でも高速道路を利用していたり、渋滞の存在がなかったり、一方通行、信号なども反映されてなかったりします。そこで、そういった条件を考慮した道路ネットワークで検証しましょうというものです。

ここでは「シミュレーションにあたては、現実の通行を再現するために、人と機械が相互に対話して経路のナビゲーションを行った」と記載されたうえで、条件設定の例は記載されているのですが、詳しい手法や反映されたGISのデータがどういったものか不明です。

ここで、乗り換え案内などのサービスを駆使すれば複雑な交通条件を考慮した経路が算定できるのでは?と思いチャレンジしてみました。

 

①経路データの入手

今回はナビタイムのAPI(Rakuten RapidAPI:一定数のアクセスまでは無料)を使って小地域の中心点から用地までの経路を検索します。

小地域のポイントデータに検索した距離をデータとして持たせ、それをシェープファイルとして出力します。

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')

作成したシェープファイルを開いて確認してみます。距離別に色を変えてみました。

②等値線を作成

ポイントデータを元にした等値線(等高線)をContourというプラグインで作成します。

※分析1で作成した3キロ圏はここで使います。広めにしておけば等値線が欠けない(閉じた図形になる)ためです。

本でも2キロ圏内は多摩川を越えていないので、似たような圏域となりました。

作成できるデータはラインデータになるので、これをポリゴンデータに変換します。できあがったポリゴンのシェープファイルを保存します。

③結果

前の2つと同様にできたシャープファイルでオーバーレイして面積按分し、人口を算出します。

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')

一応シェープファイルも見てみます。

<本に掲載されているデータ>

<今回算出したデータ>

最後の結果をみてみましょう。

1キロ圏2キロ圏
面積0.7㎢
(22.6%)
2.8㎢
(22.2%)
居住地区数4
(26.7%)
9
(15.3%)
人口7,302
(22.7%)
23,298
(16.6%)
1キロ圏2キロ圏
面積0.9㎢
(%)
4.1㎢
(%)
居住地区数11
(44.0%)
32
(38.6%)
人口10,142
(29.0%)
51,489
(31.0%)

※()内は直線距離のデータに対する割合

1キロでは似たような傾向でしたが、2キロだと乖離がありました。今回作成したものは本より西側の区域が広くなっていました。まあ、検証しようがないで。ここまで来たら結果を出せただけでも十分ということにしましょう。

 

 

まとめ

今回は少し古めの本になりますが、古今書院「GISを利用した社会・経済の空間分析」より近接性の測定にチャレンジしてみました。自分なりに着色しながらという点もありましたが、ひとまず結果を出力するところまではできたと思います。いずれ別の項目もチャレンジしてみようかと思います。

 

シェアする

コメントを残す

メールアドレスが公開されることはありません。

コメントする