![蝉蹿パッケージで読み込み
国?数値情報の公?地価(H29)のデータ(point data)を?います。
library(sf)
library(dplyr)
#フォルダ名に日本語が入っていると使えないことが多い
d <- st_read(dsn = "L01-17_13_GML", layer = "L01-17_13",stringsAsFactors = FALSE)
dat <- select(d, LPH = L01_006, CITY = L01_018, STATUS = L01_021) %>%
mutate_at("LPH",as.integer) %>%
filter(!CITY %in% c("神津島", "新島", "東京大島", "東京三宅", "八丈", "小笠原")) %>%
filter(STATUS == "住宅")
dat_mat <- t(matrix(unlist(dat$geometry), nrow=2))
dat$Y <- dat_mat[,2]
dat$X <- dat_mat[,1]
dat$ID <- 1:1360
dat2 <- as.data.frame(dat)
sfのメリットの1つはdata frameのようにdplyrで操作できることです。](https://image.slidesharecdn.com/random-170902075222/85/-6-320.jpg)
![plot(dat[1])
今?は蝉蹿パッケージに挑戦し、割と挫折した話です。](https://image.slidesharecdn.com/random-170902075222/85/-7-320.jpg)
![普遍型クリギング
LPH.gu <- gstat(id="ID", formula=LPH/10000~X+Y, locations=~X+Y, data=dat2, model=LPH.model1)
#LPH.pu <- predict(LPH.gu, mesh)
#spplot(LPH.pu[1])
sfc_POLYGONは使えないとエラーが出る](https://image.slidesharecdn.com/random-170902075222/85/-16-320.jpg)
![蝉蹿パッケージを諦める
mesh2 <- st_read(dsn = "L03-a-14_5339-jgd_GML", layer = "L03-a-14_5339") %>%
st_make_grid() %>%
st_coordinates()
mesh2 <- coordinates(mesh2)%>%
as.data.frame()
coordinates(mesh2) <- c("X","Y")
mesh2 <- as(mesh2, "SpatialPixelsDataFrame")
LPH.gu <- gstat(id="ID", formula=LPH/10000~X+Y, locations=~X+Y, data=dat2, model=LPH.model1)
LPH.pu <- predict(LPH.gu, mesh2)
## [using universal kriging]](https://image.slidesharecdn.com/random-170902075222/85/-17-320.jpg)
![推計結果
spplot(LPH.pu[1])](https://image.slidesharecdn.com/random-170902075222/85/-18-320.jpg)
![誤差
spplot(LPH.pu[2])](https://image.slidesharecdn.com/random-170902075222/85/-19-320.jpg)
![隣接?列の作成
coords <- matrix(0, nrow=length(kanto$LPH), ncol=2)
coords[,1] <- kanto$X
coords[,2] <- kanto$Y
lph.tri.nb <- tri2nb(coords)](https://image.slidesharecdn.com/random-170902075222/85/-22-320.jpg)
![夜間??密度の回帰係数
q_gwr <- classIntervals(round(kanto$gwr.popd,2), n=4, style="quantile")
q_gwr_Col <- findColours(q_gwr,pal1)
plot(kanto[20],col=q_gwr_Col,main="")
legend("bottomleft",fill=attr(q_gwr_Col,"palette"),legend=names(attr(q_gwr_Col,"table")), cex=1.1, bty="n")](https://image.slidesharecdn.com/random-170902075222/85/-30-320.jpg)
![第3次産業就業??の回帰係数
q_gwr <- classIntervals(round(kanto$gwr.emp3d,2), n=4, style="quantile")
q_gwr_Col <- findColours(q_gwr,pal1)
plot(kanto[21],col=q_gwr_Col,main="")
legend("bottomleft",fill=attr(q_gwr_Col,"palette"),legend=names(attr(q_gwr_Col,"table")), cex=1.1, bty="n")](https://image.slidesharecdn.com/random-170902075222/85/-31-320.jpg)
![標準誤差
q_gwr.se <- classIntervals(round(kanto$gwr.pred.se,2), n=4, style="quantile")
q_gwr.se_Col <- findColours(q_gwr.se,pal1)
plot(kanto[22],col=q_gwr.se_Col,main="")
legend("bottomleft",fill=attr(q_gwr.se_Col,"palette"),legend=names(attr(q_gwr.se_Col,"table")), cex=1, bty="n")](https://image.slidesharecdn.com/random-170902075222/85/-32-320.jpg)
![localR2
q_gwr.localR2 <- classIntervals(round(kanto$gwr.localR2,2), n=4, style="quantile")
q_gwr.localR2_Col <- findColours(q_gwr.localR2,pal1)
plot(kanto[23],col=q_gwr.localR2_Col,main="")
legend("bottomleft",fill=attr(q_gwr.localR2_Col,"palette"),legend=names(attr(q_gwr.localR2_Col,"table")), cex=1, bty="n")](https://image.slidesharecdn.com/random-170902075222/85/-33-320.jpg)
空间统计を使って地価分布図を描いてみる
- 3. 空間統計って? Tolber(1970)の地理学の第?法則(?rst law of geography) 全ては他の全てに関連しているが、近いものほど密接に関連している ?データ?点データ?ラスタデータ とりあえず最初は点データを使った空间内挿 具体的にはクリギング その後地理的加重回帰分析 ※今回は厳密な分析よりも?法の紹介をメインとしています。 どの?法が良いかなどは各?で調整してください。
- 4. 蝉蹿パッケージ
- 5. 今まで主流だったパッケージを使うと警告が出る library(maptools) map <- readShapeSpatial("L01-17_13_GML/L01-17_13.shp") ## Warning: use rgdal::readOGR or sf::st_read ## Warning: use rgdal::readOGR or sf::st_read
- 6. 蝉蹿パッケージで読み込み 国?数値情報の公?地価(H29)のデータ(point data)を?います。 library(sf) library(dplyr) #フォルダ名に日本語が入っていると使えないことが多い d <- st_read(dsn = "L01-17_13_GML", layer = "L01-17_13",stringsAsFactors = FALSE) dat <- select(d, LPH = L01_006, CITY = L01_018, STATUS = L01_021) %>% mutate_at("LPH",as.integer) %>% filter(!CITY %in% c("神津島", "新島", "東京大島", "東京三宅", "八丈", "小笠原")) %>% filter(STATUS == "住宅") dat_mat <- t(matrix(unlist(dat$geometry), nrow=2)) dat$Y <- dat_mat[,2] dat$X <- dat_mat[,1] dat$ID <- 1:1360 dat2 <- as.data.frame(dat) sfのメリットの1つはdata frameのようにdplyrで操作できることです。
- 8. メッシュの読み込み 国?数値情報にある3次メッシュ(raster data) mesh <- st_read(dsn = "L03-a-14_5339-jgd_GML", layer = "L03-a-14_5339") %>% st_make_grid()
- 9. ヴァリオグラムの计算
- 11. 緯度経度によるトレンド LPH.var1 <- variogram(LPH/10000~X+Y, locations=~X+Y, data=dat2) plot(LPH.var1)
- 12. ヴァリオグラム雲 plot(variogram(LPH/10000~X+Y, locations=~X+Y, data=dat2, cloud=TRUE))
- 13. 異?性モデル plot(variogram(LPH/10000~X+Y, locations=~X+Y, data=dat2, alpha=0:4*90))
- 14. 指数モデル LPH.model1 <- vgm(psill=400, model="Exp", range=0.175, nugget=35) plot(LPH.var1, LPH.model1)
- 15. 空间内挿
- 16. 普遍型クリギング LPH.gu <- gstat(id="ID", formula=LPH/10000~X+Y, locations=~X+Y, data=dat2, model=LPH.model1) #LPH.pu <- predict(LPH.gu, mesh) #spplot(LPH.pu[1]) sfc_POLYGONは使えないとエラーが出る
- 17. 蝉蹿パッケージを諦める mesh2 <- st_read(dsn = "L03-a-14_5339-jgd_GML", layer = "L03-a-14_5339") %>% st_make_grid() %>% st_coordinates() mesh2 <- coordinates(mesh2)%>% as.data.frame() coordinates(mesh2) <- c("X","Y") mesh2 <- as(mesh2, "SpatialPixelsDataFrame") LPH.gu <- gstat(id="ID", formula=LPH/10000~X+Y, locations=~X+Y, data=dat2, model=LPH.model1) LPH.pu <- predict(LPH.gu, mesh2) ## [using universal kriging]
- 20. 地理的加重回帰分析(骋奥搁)
- 21. できなかっただけで終わるのも寂しいので地理的加重回帰分析します。 国?数値情報より公?地価(H27)及び都道府県地価調査(H27)の平均 (円/m2)と国勢調査(H27)から夜間??密度(?/km2)と第3次産業就業 ??密度(?/km2)を?います。 library(spdep) library(spgwr) library(classInt) # データの読み込み poly <- read_sf(dsn="chika", layer="kanto_area") %>% mutate_at(vars(c("N03_007")),as.integer) %>% filter(FIRST_N0_2 != "所属未定地") chika <- read.csv("chika/chika.csv",header = T) #データの結合 kanto <- left_join(poly, chika, by = c("N03_007" = "CODE")) %>% mutate(area = area / 1000000) %>% mutate(EMP3D = (EMP3) / area) %>% mutate(POPD = (POP) / area)
- 22. 隣接?列の作成 coords <- matrix(0, nrow=length(kanto$LPH), ncol=2) coords[,1] <- kanto$X coords[,2] <- kanto$Y lph.tri.nb <- tri2nb(coords)
- 23. Moran?s I 地価(LPH)の空間的??相関 moran.test(kanto$LPH, nb2listw (lph.tri.nb,style="W")) ## ## Moran I test under randomisation ## ## data: kanto$LPH ## weights: nb2listw(lph.tri.nb, style = "W") ## ## Moran I statistic standard deviate = 27.025, p-value < 2.2e-16 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.867311806 -0.003105590 0.001037349
- 24. 夜間??密度(POPD)の空間的??相関 moran.test(kanto$POPD, nb2listw (lph.tri.nb,style="W")) ## ## Moran I test under randomisation ## ## data: kanto$POPD ## weights: nb2listw(lph.tri.nb, style = "W") ## ## Moran I statistic standard deviate = 25.484, p-value < 2.2e-16 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.817623729 -0.003105590 0.001037183
- 25. 第3次産業就業??密度(EMP3D)の空間的??相関 moran.test(kanto$EMP3D, nb2listw (lph.tri.nb,style="W")) ## ## Moran I test under randomisation ## ## data: kanto$EMP3D ## weights: nb2listw(lph.tri.nb, style = "W") ## ## Moran I statistic standard deviate = 25.323, p-value < 2.2e-16 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.811922251 -0.003105590 0.001035854
- 26. LPHをPOPDとEMP3Dで説明するモデルを考えます。 #バンド幅の計算 lph.bw <- gwr.sel(LPH~POPD+EMP3D, data=kanto, coords=coords) ## Bandwidth: 1.179038 CV score: 313170872880 ## Bandwidth: 1.905819 CV score: 319400361928 ## Bandwidth: 0.7298629 CV score: 297809422747 ## Bandwidth: 0.4522572 CV score: 270222406281 ## Bandwidth: 0.2806875 CV score: 240274206853 ## Bandwidth: 0.1746515 CV score: 204186527593 ## Bandwidth: 0.1091177 CV score: 170026780494 ## Bandwidth: 0.06861559 CV score: 176626753765 ## Bandwidth: 0.1014936 CV score: 167831743790 ## Bandwidth: 0.09480816 CV score: 166921437479 ## Bandwidth: 0.09173094 CV score: 166710380251 ## Bandwidth: 0.08290166 CV score: 166620328185 ## Bandwidth: 0.08627638 CV score: 166497884612 ## Bandwidth: 0.08671787 CV score: 166503113471 ## Bandwidth: 0.08602755 CV score: 166496528230 ## Bandwidth: 0.08585765 CV score: 166496319666 ## Bandwidth: 0.08589834 CV score: 166496314613 ## Bandwidth: 0.08593903 CV score: 166496344502 ## Bandwidth: 0.08589834 CV score: 166496314613
- 27. #地理的加重回帰モデルの推定 lph.gwr <- gwr(LPH~POPD+EMP3D, data=kanto, coords=coords, bandwidth=lph.bw, hatmatrix=TRUE) lph.gwr ## Call: ## gwr(formula = LPH ~ POPD + EMP3D, data = kanto, coords = coords, ## bandwidth = lph.bw, hatmatrix = TRUE) ## Kernel function: gwr.Gauss ## Fixed bandwidth: 0.08589834 ## Summary of GWR coefficient estimates at data points: ## Min. 1st Qu. Median 3rd Qu. Max. ## X.Intercept. -1757.119 13468.046 23515.427 63176.813 224307.500 ## POPD -938.010 -42.229 -20.460 10.582 249.412 ## EMP3D -686.948 18.996 109.246 195.089 3445.297 ## Global ## X.Intercept. 23177.1185 ## POPD 20.8749 ## EMP3D 0.2208 ## Number of data points: 323 ## Effective number of parameters (residual: 2traceS - traceS'S): 144.1372 ## Effective degrees of freedom (residual: 2traceS - traceS'S): 178.8628 ## Sigma (residual: 2traceS - traceS'S): 22088.53 ## Effective number of parameters (model: traceS): 110.9879 ## Effective degrees of freedom (model: traceS): 212.0121 ## Sigma (model: traceS): 20288.34 ## Sigma (ML): 16437.11 ## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 7532.024 ## AIC (GWR p. 96, eq. 4.22): 7298.536 ## Residual sum of squares: 87267699421 ## Quasi-global R2: 0.9614499
- 28. 结果を地図で表?
- 29. 推定結果を地図上で表?できるようにする。 kanto <- kanto %>% mutate(gwr.popd = lph.gwr$SDF$POPD) %>% mutate(gwr.emp3d = lph.gwr$SDF$EMP3D) %>% mutate(gwr.pred.se = lph.gwr$SDF$pred.se) %>% mutate(gwr.localR2 = lph.gwr$SDF$localR2) pal1 <- gray.colors(n=4,start=1,end=0.3)
- 30. 夜間??密度の回帰係数 q_gwr <- classIntervals(round(kanto$gwr.popd,2), n=4, style="quantile") q_gwr_Col <- findColours(q_gwr,pal1) plot(kanto[20],col=q_gwr_Col,main="") legend("bottomleft",fill=attr(q_gwr_Col,"palette"),legend=names(attr(q_gwr_Col,"table")), cex=1.1, bty="n")
- 31. 第3次産業就業??の回帰係数 q_gwr <- classIntervals(round(kanto$gwr.emp3d,2), n=4, style="quantile") q_gwr_Col <- findColours(q_gwr,pal1) plot(kanto[21],col=q_gwr_Col,main="") legend("bottomleft",fill=attr(q_gwr_Col,"palette"),legend=names(attr(q_gwr_Col,"table")), cex=1.1, bty="n")
- 32. 標準誤差 q_gwr.se <- classIntervals(round(kanto$gwr.pred.se,2), n=4, style="quantile") q_gwr.se_Col <- findColours(q_gwr.se,pal1) plot(kanto[22],col=q_gwr.se_Col,main="") legend("bottomleft",fill=attr(q_gwr.se_Col,"palette"),legend=names(attr(q_gwr.se_Col,"table")), cex=1, bty="n")
- 33. localR2 q_gwr.localR2 <- classIntervals(round(kanto$gwr.localR2,2), n=4, style="quantile") q_gwr.localR2_Col <- findColours(q_gwr.localR2,pal1) plot(kanto[23],col=q_gwr.localR2_Col,main="") legend("bottomleft",fill=attr(q_gwr.localR2_Col,"palette"),legend=names(attr(q_gwr.localR2_Col,"table")), cex=1, bty="n")
- 34. 参考?献など Arbia, Giuseppe. 2016. 『Rで学ぶ空間計量経済学??』. 勁草書房. 古?知之. 2011. 『Rによる空間データの統計分析』. 朝倉書店. 岩波データサイエンス刊?委員会. 2016. 『岩波データサイエンス vol.4 : 特集地理空間情報処理』. 岩波書店. 瀬?創 ? 堤盛?. 2014. 『空間統計学 : ?然科学から???社会科学 まで』. 朝倉書店. ?村晋. 2010. 『地理空間データ分析』. 共?出版. Tobler, W. R. 1970. “A Computer Movie Simulating Urban Growth in the Detroit Region.” Economic Geography 46: 234–40. 堤盛? ? 村上?輔 ? 嶋?章. 2014. 「我が国の三?都市圏を対象とし た住宅地価分布図の作成.」 『GIS-理論と応?』 22(2): 69–79.
- 36. ご清聴ありがとうございました