第10回演習マニュアル

演習の目的

  • ポイントサンプリングされた環境データを使って,空間補間を行う.

Rの起動とパッケージの読み込み

  • Rを起動し,spdepとgstatのパッケージを読み込みます.
    library(spdep)
    library(gstat)

データの読み込み

  • データをここからダウンロードして解凍してください。以下では、解凍したフォルダを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")

シェープファイルを表示する

  • 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)による補完

  • 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")

バリオグラムの計算

  • まず,全ての方向で等しい(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)

クリギング

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

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

通常型クリギング

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

他のシェープファイルを表示する

  • SPMの観測地点データ「YokoSPMj.shp」,区境界「Ward.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")

 left

単純型クリギング

  • 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")

 left

普遍型クリギング

  • 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")

 left

他のデータのクリギング補間

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

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


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-06-22 (火) 12:12:09 (3095d)