第07回演習マニュアル

演習の目的

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

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

  • Rを起動し,「spdep」,「sp」,「gstat」のパッケージを読み込みます.
  • 「パッケージ」→「パッケージの読み込み」→「spdep」,「sp」及び「gstat」を選んでください.
  • マニュアル:spspdep,gstat
  • 今回は,Rコマンダーは使いません.

データの読み込み

  • 第05回演習で使った「yoko」フォルダが,Rフォルダの下にあることを確認してください.
  • 「yoko」フォルダから,以下のファイルを読み込んでください.
> so2 <- read.table("yoko/YokoSO2j.csv",header=TRUE,sep=",")
> spm <- read.table("yoko/YokoSPMj.csv",header=TRUE,sep=",")
  • 「so2」は 横浜市のSO2観測データを,「spm」はSPM(浮遊粒子状物質)の観測データです.

バリオグラムの計算

  • まず,全ての方向で等しい(isotropy)バリオグラムを作成する.
    > spm.var<-variogram(object=ITEM060~1,locations=~X+Y,data=spm)
    > plot(spm.var)
  • まず,見た目でシル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のクリギングについて説明します.

通常型クリギング

  • このmodel2を使って,まず通常型クリギングにより空間補間することにしましょう.そのためにまず,メッシュデータの座標を読み込みます.
    > mesh.grid <- read.table("yoko/mesh.csv",header=TRUE,sep=",")
  • 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 <- read.shape("yoko/YokoSPMj.shp")
Shapefile Type: Point   # of Shapes: 43
> plot(spm.shp)
> ward.shp <- read.shape("yoko/Ward.shp")
Shapefile Type: Polygon   # of Shapes: 24
> ward.map <- Map2poly(ward.shp)
> 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.map,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)
      [using simple kriging]
      > levelplot(exp(ID.pred)~X+Y,spm.ps,add=TRUE,main="SPM Simple Kriging Predict in Yokohama")
      > plot(spm.shp,add=T)
      > plot(ward.map,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)
[using universal kriging]
> levelplot(exp(ID.pred)~X+Y,spm.pu,add=TRUE,main="SPM Universal Kriging Predict in Yokohama")
> plot(spm.shp,add=T)
> plot(ward.map,add=T,border="gray")

 left

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

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

演習課題レポートの提出

  • 地価データを使ったクリギング補間結果についてまとめなさい.少なくとも,単純型,通常型,普遍型の3つのクリギング手法を適用し,更に普遍型クリギングについては,傾向面(回帰モデルの説明変数)に応じた推定結果を示しなさい.
    • gstat()関数やlevelplot()を適用する際に,exp(ID.pred)ではなく,単にID.predとしたほうが,よい場合があります.
  • 分量:A4 ページ数は特に指定しないが,読みやすくなるよう心がけること.
  • 提出期限・方法:12月11日(火)の授業時間中にレポートを提出すること.授業時間内の提出が困難な場合,提出期限である授業開始前までに担当者にメールで送付すること.

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


トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2008-05-22 (木) 08:35:38 (4077d)