QGIS, Rを用いた公示地価データの空間統計分析2

ここでは、QGIS, Rを用いて空間統計学の基礎的な演習を説明しています。

R, QGISのインストール

以下のサイトからインストーラーをダウンロードし、インストールする。

使用するデータ

データの前処理

空間統計解析の前に、上記のデータを処理する必要がある。

（市町村ポリゴン）京都市及び周辺都市の抽出 (QGIS)

QGISを起動する。本来は不要かもしれないが、念のために以下の処理を済ませておく。

新規プロジェクトを選択し、何もレイヤを追加していない段階で、座標参照系 (coordinate reference system: CRS)を設定する。

「japan_ver83.shp」は日本測地系 2011(JGD2011)の座標系で作成されている。これを平面直角座標系VI系に変換する。 注意：座標系を変換するには、別名で保存する必要あり。

その後、右クリックで「レイヤの削除」を選んで一旦画面から消去する。その後、「japan_ver83_VI.shp」を読み込めば、平面直角座標系VI系で表示されるようになる。

注意：上記画像の右下のように、座標が平面直角座標系（m単位）で表示されているか確認する。度の単位だとこの時点で何か誤りが発生していると思われる。

その後、選択した複数の市を別のshapeファイルとして保存する。

以下の処理のために、「japan_ver83_Kyoto.shp」は残し、平面直角座標系第VI系で表示されているはずである。「japan_ver83_VI.shp」は「レイヤの削除」で削除しておく。

CRSの設定後、ウィンドウが真っ白になる場合がある。その場合には焦らず、「（メニューバー）ビュー」→「全域表示」下記のように指定すれば所定のデータが表示される。

（公示地価データ）選択地域で切り抜き (QGIS)

公示地価データ「L01-21_26.shp」は京都府全域の公示地価データが含まれている。上記の対象地域だけのデータを抽出するために下記の処理を行う。

この後、CRSに平面直角座標系第IV系を指定した別のshapeファイルとして保存する。

（公示地価データ）緯度、経度座標の出力 (QGIS)

上記の平面直角座標系第IV系の公示地価データ（L01-21_26_VI）に対して処理を進めていく。属性テーブルを開き、下記の処理をフィールド名「X」と「Y」（共にm単位）について行う。

（公示地価データ）特定の地点までの距離の計算 (QGIS)

後の重回帰分析において特定地点との直線距離（本来なら移動に要する時間距離が望ましい）を用いる為、ここで計算する。以下、四条駅の平面直角座標系VI系での座標値（単位はm）を記す。これを利用して、上記で作成した「L01-21_26_VI.shp」の属性テーブル内に新たなフィールドを設け、距離を格納する。

（公示地価データ）不要なデータの削除、列名の変更 (Excel)

最後に、Excel上で今回の演習に関係のない列を削除して、列の項目名を変更する。下記に示す列だけを残し、1行目に記された列名を変更した後に、CSVファイル「L01-21_26_Kyoto_edited.csv」として保存する。

重回帰モデルのパラメータ推定

install.packages("mgcv")
install.packages("RColorBrewer")
install.packages("sp") 
library(mgcv)
library(RColorBrewer)
library(sp) 

# CSVファイルの読み込み
dat <- read.csv("L01-21_26_Kyoto_edited.csv")

# 奇数：モデル作成用, 偶数：検証用
dat_tmp <- split(dat, rep(c("odd", "even")))
dat1 <- dat_tmp$odd # モデル作成用
dat2 <- dat_tmp$even # 検証用

# 各年度のデータを下記の並びに整理
# Station, Dist, Price, X, Y, Dist2Shijo, Year

dat1a <- data.frame(Station=dat1$Station, Dist=dat1$Dist, Price=dat1$p2010, 
X=dat1$X, Y=dat1$Y, Dist2Shijo=dat1$Dist2Shijo, Year=2010)

dat2a <- data.frame(Station=dat2$Station, Dist=dat2$Dist, Price=dat2$p2010, 
X=dat2$X, Y=dat2$Y, Dist2Shijo=dat2$Dist2Shijo, Year=2010)

year_list <- 2011:2021

for (year in year_list) {
	
	# "p2011", "p2012"のような列名を作成
	yearid <- paste0("p", year)

	dat1b <- data.frame(Station=dat1["Station"], Dist=dat1["Dist"], 
        Price=dat1[yearid], X=dat1["X"], Y=dat1["Y"], 
        Dist2Shijo=dat1["Dist2Shijo"], Year=year)
	
	dat2b <- data.frame(Station=dat2["Station"], Dist=dat2["Dist"], 
        Price=dat2[yearid], X=dat2["X"], Y=dat2["Y"], 
        Dist2Shijo=dat2["Dist2Shijo"], Year=year)
	

	# 上記のままでは列名が"Price"ではなく、"p2001"のように代入されるため、
	# 下記の処理で"Price"へ置換
	dat1b <- dplyr::rename(dat1b, Price=yearid)
	dat2b <- dplyr::rename(dat2b, Price=yearid)

	dat1a <- rbind(dat1a, dat1b)
	dat2a <- rbind(dat2a, dat2b)
}

# Priceを対数化
dat1a$Pricelog <- log(dat1a$Price)
dat2a$Pricelog <- log(dat2a$Price)

# Pricelogが無効な行を削除
dat1 <- dat1a[!(dat1a$Pricelog == "-Inf"),]
dat2 <- dat2a[!(dat2a$Pricelog == "-Inf"),]


# モデル用関数1
f1 <- Pricelog ~ te(X, Y, Year, k=c(50, 10), d=c(2, 1)) 
+ s(Year) + s(Dist) + s(Dist2Shijo)
model1   <- gam(f1, data = dat1)

# 結果の出力1
sink("output_f1.txt")
summary(model1)

jpeg("fig1-1.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model1, select=1)
dev.off()

jpeg("fig1-2.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model1, select=2,ylim=c(-1,1))
dev.off()

jpeg("fig1-3.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model1, select=3,ylim=c(-1,1))
dev.off()

jpeg("fig1-4.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model1, select=4)
dev.off()

# 検証1
pred01  <- predict(model1, dat2, se.fit = TRUE)
pred1   <- data.frame(dat2[,c("X","Y")], 
pred=pred01$fit, se.pred=pred01$se.fit)
pred
sink()

# jetカラーマップの作成
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", 
"cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))


cuts_val <- c(8, seq(10, 13, 0.25), 14)
cuts_err <- c(0, seq(0.03, 0.1, 0.01), 0.2)

coordinates(pred1)<-c("X","Y")

jpeg("fig1-val.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")

spplot(pred1, "pred", xlab="Lat (deg)", ylab="Lon (deg)", cex.lab=1.5, 
cex.axis=2.0, cuts=cuts_val, col.regions=jet.colors(7), cex=2.5, pch=20, 
main="Estimated land price (log): Year, Dist, Dist2Shijo")

dev.off()

jpeg("fig1-error.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")

spplot(pred1, "se.pred", xlab="Lat (deg)", ylab="Lon (deg)", cex.lab=1.5, 
cex.axis=2.0, cuts=cuts_err, col.regions=jet.colors(7), cex=2.5, 
pch=20, main="Stardard error (log) distribution: Year, Dist, Dist2Shijo")

dev.off()


# モデル用関数2：四条駅までの距離を説明変数から除外
f2 <- Pricelog ~ te(X, Y, Year, k=c(50, 10), d=c(2, 1)) + s(Year) + s(Dist)
model2 <- gam(f2, data = dat1)

# 結果の出力2
sink("output_f2.txt")
summary(model2)

jpeg("fig2-1.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model2, select=1)
dev.off()

jpeg("fig2-2.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model2, select=2,ylim=c(-1,1))
dev.off()

jpeg("fig2-3.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")
plot(model2, select=3, ylim=c(-1,1))
dev.off()


# 検証2
pred02  <- predict(model2, dat2, se.fit = TRUE)
pred2   <- data.frame(dat2[,c("X","Y")], pred=pred02$fit, se.pred=pred02$se.fit)
pred
sink()

cuts_val <- c(8, seq(10, 13, 0.25), 14)
cuts_err <- c(0, seq(0.03, 0.1, 0.01), 0.2)

coordinates(pred2)<-c("X","Y")

jpeg("fig2-val.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")

spplot(pred2, "pred", xlab="Lat (deg)", ylab="Lon (deg)", cex.lab=1.5, 
cex.axis=2.0, cuts=cuts_val, col.regions=jet.colors(7), cex=2.5, pch=20, 
main="Estimated land price (log): Year, Dist")

dev.off()

jpeg("fig2-error.jpg", units = "px", pointsize = 24, quality = 80, bg = "white")

spplot(pred2, "se.pred", xlab="Lat (deg)", ylab="Lon (deg)", cex.lab=1.5, 
cex.axis=2.0, cuts=cuts_err, col.regions=jet.colors(7), cex=2.5, pch=20, 
main="Stardard error (log) distribution: Year, Dist")

dev.off()

QGIS上での結果の表示

推定結果をQGIS上で表示する。上記の出力結果「L01-21_26_Kyoto_rst1.csv（モデル1）」「L01-21_26_Kyoto_rst2.csv（モデル2）」をGIS上に取り込でいく。その前に、各列の名前（フィールド名）が存在しないので、Excelでもテキストエディタでも結構なので、下記のように命名する。

次に、他のデータと重ね合わせて表示させるためにも、上記のファイルに座標系を与えて保存する。

この時点で表示されている「Kyoto10_Price_rst」は表示画面から削除する。最終結果を分かりやすく表示させるために、下記のデータも事前にダウンロードし表示させておく。（注意：別名で保存して、座標系の変換が必要）

最後に、上記で作成した「L01-21_26_Kyoto_rst1.shp」を表示させ、標準誤差に応じて表示させる色を変えていく。

課題

3人一組のグループで以下の課題に取り組んで、3人一組で一つのレポートを作成して提出する。


モデル1による推定の標準誤差	モデル2による推定の標準誤差