* 第10回演習マニュアル [#q4b5460c]

*** 演習の目的 [#k45bd5b4]
- ポイントサンプリングされた環境データを使って,空間補間を行う.

** Rの起動とパッケージの読み込み [#f0c4898f]
- Rを起動し,spdepとgstatのパッケージを読み込みます.
 library(spdep)
 library(gstat)

- マニュアル:[[sp>http://cran.r-project.org/doc/packages/sp.pdf]],[[spdep>http://cran.r-project.org/doc/packages/spdep.pdf]],[[gstat>http://cran.r-project.org/doc/packages/gstat.pdf]]

** データの読み込み [#u9c8e7b1]
- データを[[ここ>http://web.sfc.keio.ac.jp/~maunz/SI07/yoko.zip]]からダウンロードして解凍してください。以下では、解凍したフォルダをyokoとしています。
- 「yoko」フォルダから,以下のファイルを読み込んでください.
 so2 <- read.table("YokoSO2j.csv",header=TRUE,sep=",")
 spm <- read.table("YokoSPMj.csv",header=TRUE,sep=",")
- 「so2」は 横浜市のSO2観測データを,「spm」はSPM(浮遊粒子状物質)の観測データです.
- 次に,点データを補完するメッシュデータの座標を読み込みます.
 mesh.grid <- read.table("mesh.csv",header=TRUE,sep=",")
 coordinates(mesh.grid) <- c("X","Y")
 mesh.grid <- as(mesh.grid, "SpatialPixelsDataFrame")

** シェープファイルを表示する [#aeaf26de]
- SPMの観測地点データ「YokoSPMj.shp」,区境界「Ward.shp」,およびクリギング結果を読み込みます.
 # spm.shp <- read.shape("YokoSPMj.shp")
 spm.shp <- readShapePoints("YokoSPMj.shp")
 # ward.shp <- read.shape("Ward.shp")
 ward.shp <- readShapePoly("Ward.shp")
 # ward.map <- Map2poly(ward.shp) 
 plot(spm.shp)
 plot(ward.shp,add=T)

** 逆距離加重法(IDW)による補完 [#i4f8f192]
- idw()関数を使う方法
 spm.idw1 <- idw(formula=ITEM060~1, locations=~X+Y, data=spm, mesh.grid, idp=2)
 plot(spm.shp) 
 spplot(spm.idw1["var1.pred"],add=T)
 plot(spm.shp,add=T)
 plot(ward.shp,add=T,border="gray")

- gstat()関数を使う方法
 spm.idw2 <- gstat(id="ID",formula=ITEM060~1, locations=~X+Y,data=spm, 
 set=list(idp=2))
 spm.idw.p <- predict(spm.idw2,mesh.grid)
 plot(spm.shp) 
 spplot(spm.idw.p["ID.pred"],add=T)
 plot(spm.shp,add=T)
 plot(ward.shp,add=T,border="gray")

** バリオグラムの計算 [#gd75f989]
- まず,全ての方向で等しい(isotropy)バリオグラムを作成する.
 spm.var<-variogram(object=ITEM060~1,locations=~X+Y,data=spm)
 plot(spm.var)
- バリオグラム雲を表示する
 spm.var.cl<-variogram(object=ITEM060~1,locations=~X+Y,data=spm,cloud=TRUE)
 plot(spm.var.cl)
- まず,見た目でシルsillとレンジrange,ナゲットnugget効果を決めて見ましょう.ここではガウス型モデルを適用する.
 model1=vgm(psill=0.000015,model="Gau",range=8000,nugget=0.000035)
 plot(spm.var,model=model1)
- 次に,定量的にフィットしたバリオグラムを作成しよう.
 model2=fit.variogram(object=spm.var,model=vgm(psill=0.000015,model="Gau",
 range=8000,nugget=0.000035))
 plot(spm.var,model=model2)
- このモデルだと,model2はmodel1のシルとレンジを固定して使っていますが,これを固定せずに使うには,var,model= vgm()に続いて以下のオプションを指定してください.
 fit.sills = FALSE, 
 fit.ranges = FALSE,
- また,vgm()関数の中で,以下のオプションを指定すると,誤差の二乗和が表示されるので,これが最も小さいヴァリオグラム=最も当てはまりがよいモデルを適用するとよいでしょう.
 print.SSE = TRUE,
- 異方性(anisotropy)バリオグラムは,次のようにalphaを指定して計算できます.
 spm.var3<-variogram(object=ITEM060~1,locations=~X+Y,data=spm,
 alpha=c(0,45,90,135))
 plot(spm.var3)

** クリギング [#o63e163d]
- 主要なクリギングとして,以下の4つのクリギングがあり,R言語でもこれらを扱えます. 
+ 単純型クリギング
-- 単純に,ある地点の属性値と標本全体の平均値を使って推定する.
+ 通常型クリギング
-- バリオグラムを取り入れて推定する.一般的なクリギング手法.
+ 普遍型クリギング
--回帰モデルのように,傾向面と残差を考えて推定する.
+ コクリギング
-- 既知の値Zとそれとは別の既知の値Yを用いて,未知のZ値を推定する.

以下では,1〜3のクリギングについて説明します.

** 通常型クリギング [#w2179f22]
- gstatライブラリでは,gstat関数を使って,データやバリオグラムモデルを管理することができます.さらに,predict.gstat関数を使うと,グリッドの座標位置を引数として予測を行える.以下の方法で,通常型クリギングを実行できる.
 spm.g <- gstat(id="ID",formula=ITEM060~1,locations=~X+Y,data=spm,model=model2)
 spm.p <- predict.gstat(spm.g,mesh.grid)

** 他のシェープファイルを表示する [#aeaf26de]
- SPMの観測地点データ「YokoSPMj.shp」,区境界「Ward.shp」,およびクリギング結果を表示します.
 spm.shp <- read.shape("YokoSPMj.shp")
 spm.shp <- readShapePoint("YokoSPMj.shp")
 plot(spm.shp)
 ward.shp <- readShapePoly("Ward.shp")
 # ward.map <- Map2poly(ward.shp)
 spplot(spm.p["ID.pred"],add=T)
 # library(lattice)
 # levelplot(exp(ID.pred)~X+Y,spm.p,add=TRUE,
 # main="SPM Kriging Predict in Yokohama")
 plot(spm.shp,add=T)
 plot(ward.shp,add=T,border="gray")

&ref(http://web.sfc.keio.ac.jp/~maunz/wiki/data/sdm06_07/image002.jpg, left);~
~

** 単純型クリギング [#ncc6ac92]
- gstatの引数にbeta(平均)を定義すれば単純型クリギングを推定できます.
-- exp()はなくてもいいです.(以下同じ)
 mean(spm$ITEM060)
 [1] 0.03960294
 spm.gs <- gstat(id="ID",formula=ITEM060~1,locations=~X+Y,data=spm,
 model=model2,beta=mean(spm$ITEM060))
 spm.ps <- predict.gstat(spm.gs,mesh.grid)
 # levelplot(exp(ID.pred)~X+Y,spm.ps,add=TRUE,main="SPM Simple Kriging Predict in Yokohama")
 plot(spm.shp)
 spplot(spm.ps["ID.pred"],add=T)
 plot(spm.shp,add=T)
 plot(ward.shp,add=T,border="gray")

&ref(http://web.sfc.keio.ac.jp/~maunz/wiki/data/sdm06_07/image005.jpg, left);~
~

** 普遍型クリギング [#r3b366c8]
- variogram関数のformulaの引数(説明変数)を~1ではなく,色々指定すると,普遍型クリギングを推定できる.以下の式ではX+Yが傾向面となる.
 spm.gu <- gstat(id="ID",formula=ITEM060~ X+Y,locations=~X+Y,
 data=spm,model=model2)
 spm.pu <- predict.gstat(spm.gu,mesh.grid)
 # levelplot(exp(ID.pred)~X+Y,spm.pu,add=TRUE,main="SPM Universal Kriging Predict in Yokohama")
 plot(spm.shp)
 spplot(spm.pu["ID.pred"],add=T)
 plot(spm.shp,add=T)
 plot(ward.shp,add=T,border="gray")

&ref(http://web.sfc.keio.ac.jp/~maunz/wiki/data/sdm06_07/image007.jpg, left);~
~

** 他のデータのクリギング補間 [#v299b26b]
- 地価やSO2の分布データを使ってクリギング補間を行って見ましょう.
- 地価データの場合,重回帰モデルで構築した際の説明変数を傾向面としてformulaを定式化し,普遍型クリギングを適用してみてください.

以上で今回の演習は終了です.

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS