#始めに
大阪府を対象に、過去（２０００年〜２０１０年）にあった土地利用や人口数などの街の変化と、駅利用者数の間の関係を記述する試みをしました。利用するデータのほとんどはEstatと国土数値情報で公開されているものです。ツールは可視化用のQGISとデータ処理用のR言語を使います。
#地理的回帰加重分析について
地理空間加重回帰分析法（Geographically weighted regression,GWR）は複数の変数間の関係を調べるために使われている。線形回帰が一つの回帰係数で対象地を全般的に説明するのに対して、GWRは空間内の異なる場所にに異なる関係を持つことを認め、空間的非定常性（Exploring Spatial）が考慮されるので、回帰係数やt-valueの空間的分布を地図化することができます。
#各メッシュに対してこのようなデータを準備する
###メッシュごとに各土地利用の面積を計算する
元データ：
国土地理院数値地図5000（土地利用）：https://www.gsi.go.jp/kankyochiri/lum-5k.html
５００mメッシュ：https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=Q&coordsys=2&format=shape
重ねてみるとこのような感じになります。
Clipすると、土地利用のポリゴンがメッシュで分割されます。また、メッシュコードも全てのポリゴンに付与されます。
メッシュコード、面積、土地利用分類コードを持ったデータを元にメッシュごとに各土地利用分類の割合（イメージは以下のような）を計算します。
|1.商業・業務用地|2.密集低層住宅地|...|７.その他|メッシュコード|
|-|---|---|---|---|---|表示が省略されている列であれば、アラインメント行は好きに書いてよい。|
|....
|
# データを読み込み
landuse2001=read.csv('landuse2001_main.txt!',sep = ',',header = T)
# 重複しているメッシュコードを削除する
code=landuse2001$mesh_500m_keycode[!duplicated(landuse2001$mesh_500m_keycode)]
t.area.f2001=data.frame()
m=1
i=1
for (m in 1:length(code)) {
k=code[m]
i=i
for (i in 1:7) {
index=which(landuse2001$mesh_500m_keycode==k&landuse2001$lu_name_1==i)
t.area.f2001[as.character(k),i]=sum(landuse2001[index,]$Shape_Area)
i=i+1
}
m=m+1
}
t.area.f2001$code=rownames(t.area.f2001)
names(t.area.f2001)=c('2001_1','2001_2','2001_3','2001_4','2001_5','2001_6','2001_7','code')
メッシュごとの各土地利用の面積が得られます。
> head(t.area.f2001)
2001_1 2001_2 2001_3 2001_4 2001_5 2001_6 2001_7 code
513530284 0 0 423.9092 0 0 0.0000 1308.2131 513530284
513530293 0 0 0.0000 0 0 0.0000 1444.3341 513530293
513530294 0 0 0.0000 0 0 0.0000 525.2058 513530294
513530382 0 0 1387.9418 0 0 0.0000 3399.8467 513530382
513530384 0 0 1259.7997 0 0 0.0000 4938.8111 513530384
513530391 0 0 4057.4594 0 0 316.5574 9943.3534 513530391
QGISで計算した各メッシュの面積を使って絶対値を割合に変換し、土地利用データの準備は完了します。
###メッシュごとに土地利用以外のデータをマージする
Estatで５００mメッシュに合わせた人口数、事業所数、従業者数のデータもマージしてみました。（https://www.e-stat.go.jp/gis/statmap-search?type=1）
マージする流れは全く同じなので、事業所数データの例だけ出します。
alldata = t.area.f2001
mesh.jigyou20015135=read.csv("jigyou2001H5135_500m.txt",sep = ",",header = T)
colnames(mesh.jigyou20015135)=c("code","事業所数2001","従業者数2001")
mesh.jigyou20015135=mesh.jigyou20015135[-1,]
i=1
for (i in 1:dim(alldata)[1]) {
t=alldata$code[i]
alldata$事業所数2001[i]=mesh.jigyou2001[which(mesh.jigyou2001$code==t),][2]
alldata$従業者数2001[i]=mesh.jigyou2001[which(mesh.jigyou2001$code==t),][3]
i=i+1
}
alldata$事業所数2001=as.numeric(alldata$事業所数2001)
###最寄り駅と間の距離を求める
国土数値情報の鉄道駅情報を利用します。（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html）
QGISでメッシュデータと重ねて表示し、水色のラインは駅となります。プロセシングツールの「最寄りハブの距離」を使うことによって、入力レイヤの地物から最も近いハブレイヤの地物を特定し、距離を計算することもできます。ポリゴンとポリラインに対して距離は地物の中心から計算されます。
さらに、駅別乗降客数データ（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-S12-v2_7.html）を追加します。
各メッシュの最寄り駅を特定することによって、駅の利用者数と対応させます。
###出来上がったデータ
２００１年データの処理手順と同じく２００８年のデータもマージしたので、最後に出来上がったデータの行列このようになります。一つのメッシュに対してメッシュ内の情報と最寄り駅の情報が格納されています。
|メッシュコード|人口総数2000|男2000|女2000|世帯総数2000|人口総数2010|男2010|女2010|世帯総数2010|事業所数2001|従業者数2001|事業所数2006|従業者数2006|scode_2000|neardis_2000|scode_2010|neardis_2010|駅名_2000|運営会社_2000|乗車_2000|降車_2000|駅名_2010|運営会社_2010|乗車_2010|降車_2010|X2001_1|X2001_2|X2001_3|X2001_4|X2001_5|X2001_6|X2001_7|X2008_1|X2008_2|X2008_3|X2008_4|X2008_5|X2008_6|X2008_7|
-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|
|....
|513530384|NA|NA|NA|NA|2336|1249|988|782|NA|NA|NA|NA|21840|4111.64|8000060140|4111.64|孝子|南海|194|224|孝子|20|141|94|21840|4111.64|8000060140|4111.64|孝子|南海|194|224|孝子|20|141|94|0|0|0.0078|0|0|0|0.0192|0|0|0.0087|0|0|0|0.0192|
|513530391|1054|393|3|1343|480|1450|1534|1021|NA|NA|1|2|21840|3585.59|8000060140|3585.59|孝子|南海|194|224|孝子|20|141|94|0|0|0.0153|0|0|0.0012|0.0374|0|0|0.0153|0|0|0.0012|0.0374|
|....
#作成したデータを利用して地理的回帰加重分析（GWR）を行う
解析用のモデルGWmodelをインストールします。
install.packages('GWmodel')
library(GWmodel)
GWRは空間データをベースに行うので、データをメッシュShapeにマージし、Shapefileとして再度読み込む必要があります。GISToolsモデルを使うと便利です。
install.packages('GISTools')
library(GISTools)
gwrshp=readShapePoly('GWR.shp')
駅利用者数は従属変数で、その他は独立変数として、GWR分析を行います。
gwrshp2000mesh_res=gwr.basic(利用者2000~X2001_1+X2001_2+X2001_3+X2001_4+X2001_5+X2001_6+X2001_7+人口総数2000+事業所数2001+従業者数2001+neardis_2000,data=gwrshape,kernel='gaussian',bw=1000)
###分析結果
GWRの特徴である、距離によってt-valueの空間的分布ということは、各メッシュに対して一つの回帰式が作成され、メッシュごとのt-valueも出力されます。ｔ値の絶対値が大きいほど、説明変数が目的変数に与える影響が強いことを意味しています。今回t-valueで評価する際に、ｔの絶対値が1.64より大きな場合は、モデル式が90%の精度を有するという基準で変数と説明変数の相互関係を判断します。
例えば、「商業・業務用地」の割合と最寄り駅利用者数の完成性を知りたい場合、「商業・業務用地」独立変数のt-valueを使ってQGISで色つけてみました。赤色で示されたところにおいて、最寄り駅の利用者数との正の相関関係、青色は負の相関関係、黄色の場合は有意義な相関関係がないことが示されています。
正の相関関係（赤色）がある区域を見れば、2000年には大阪市の福島区、中央区、西区、浪速区を主とした区域と、東淀川区、吹田市南部を主とした区域が北と南の２エリアに分けており、真ん中の負の相関関係（青色）がある北区と明確な相関関係のない淀川区等が挟まれていた。2010年になると、東淀川区から浪速区までの中心エリアが目立つようになり、正の相関関係を持っていた吹田市周辺エリアが、負に転換しており、商業・業務用地の増加が駅利用者数の減少をもたらすことが示されています。駅に対する地域の機能の大きな転換が見られ、市中心部への集中に伴う外縁部の衰退が考えられます。また、2000年と比べると、2010年の方が大阪市より南の沿岸部において、幾つかのエリア新たに現れたことが分かります。南沿岸部、特に乗り換え駅の近くにおいて、商業・業務用地と駅利用者の関係が強まっていったことが考えられます。
###終わりに
今回の記事では、メッシュデータに基づいて地理的加重回帰の手法を介し、駅の利用者数と駅周辺土地利用などの状況の関係の強さと関連する範囲を記述してみました。メッシュを利用しているので、今回のデータに限らず、メッシュに格納できるデータを同じ手法での分析を行うことも期待できます。