緯度経度情報を含む平面的なデータがあるとき、やはり真っ先に「俯瞰的に情報を把握したい」と思うのではないでしょうか。
そんなときに便利なのが等値線図 (コンター) ですよね。
でもこの等値線図、どうやって描きますか？
Surfer に GMTこの他にも様々なソフトがありますよね
もちろん R にも contour() という関数があり、等値線図を描くことができます。しかし、contour() を利用する前には補間をしておかなければなりません。補間の種類は色々とあるのですが、今回はクリギングができるパッケージを紹介したいと思います。
そして等値線図を R で表示しただけでは、海岸線が表示されていないので、あまり意味がありません。
そこで、今回は、Google Earth に等値線図を貼り付けるコードを紹介したいと思います。
出力イメージはこんな感じです。
なにはともあれデータが必要なのでこちらのテストデータをご利用下さい。
このデータは前回 Google Earth にプロットした九州西部の AMeDAS の緯度経度データと気温のデータです
データは 2014 年 3 月 30 日 12:00 のものです。
必要となるパッケージは automap です。
それと表示用の Google Earth です
使い方 ---------------------------
はじめにプログラム上部に、変更可能な変数があるので、任意の値を設定して下さい。
xgn = 300 はテストデータの西端と東端の間を 300 等分したときの点という意味です。
つまり、この値を増やすほど詳細に計算されますが、計算速度は当然落ちますし、結果のファイルも重たくなります。
これに関連して、reso で解像度を設定できます。reso = 1 のときは 480×480 ピクセルです。
これもあまりに大きくすると、綺麗ですが重たくなりますので、使用している PC のスペックに合わせて設定して下さい。
色の設定ではデフォルトで、白 → オレンジ に変わるようになっています。
R で指定できる色に変更可能ですし、カンマで切って繋げれば、複数設定することもできます。
例えば c
プログラムを実行すると、ファイル選択ウィンドウが開くので、テストデータあるいはご自分で作成したファイルを選んで下さい。
以上です。
読み込んだファイルがあるフォルダに「ファイル名_等値線図.kml」というファイルができているはずです。
これを開くと等値線図が見れます。
プログラムについて ---------------
補間に関しては、パッケージ automap の関数 autoKrige() に任せっきりです。
どのようなモデルで補間するかとかも指定できるので automap.pdf で確認して下さい。
デフォルトでは model="Gau" です。
途中でバリオグラムが表示されるので、正しく補間されているか確認して下さい。
特に自作関数は使ってないので、気になるところは一行ずつ実行してもらうと中身が分かると思います。
#作成日時 : 9:18 2011/07/06
#作成者名：T. Nakano
rm(list=ls())
gc();gc()
library(automap)
# //////////////// Google Earth で等値線図作成 ///////////////////
# 変更可 -----------------------------------------------------------
# 等値線の本数 ----------
LN = 10
# 等値線の値の桁数 -----
dgt = 1
# 外挿計算 --------------
extrapolate = "off" # on or off
# 補間点数 --------------
xgn = 300 # x 軸方向に xgn 分割する
ygn = 300 # y 軸方向に ygn 分割する
# 解像度 ----------------
reso = 2
# 色の設定 --------------
cols = colorRamp(c("white", "orange"))
# データの読込み -------------------------------------------------
library(tcltk)
fileName <- tclvalue(tkgetOpenFile())
act1 <- unlist(strsplit(fileName,"/"))
act1 <- act1[-length(act1)]
directory <- NULL
if(length(act1)!=1){
for(i in 1:(length(act1)-1)){
directory <- paste(directory,act1[i],"/",sep="")
}
}
directory <- paste(directory,act1[length(act1)],sep="")
setwd(directory)
SRX <- read.csv(fileName, header=T)
DT <- data.frame(x=SRX$LON, y=SRX$LAT, data=SRX$DATA)
coordinates(DT) =~ x+y
# 補間地点の緯度経度整理 -----------------------------------------
Xst <- seq(min(SRX$LON), max(SRX$LON), length=xgn)
Yst <- seq(min(SRX$LAT), max(SRX$LAT), length=ygn)
Xip <- rep(Xst, each=length(Yst))
Yip <- rep(Yst, length(Xst))
Ncode <- data.frame(x=Xip, y=Yip)
# クリギング -------------------------------------------------------
gridded(Ncode)=~x+y
kriging_result <- autoKrige(data~1, DT, Ncode, model="Gau")
WR <- data.frame(LON=Xip, LAT=Yip, PRED=kriging_result$krige_output$var1.pred)
plot(kriging_result)
# 凸閉包で外挿除外 -------------------------------------------------
if(extrapolate=="off"){
act1 <- data.frame(x=SRX$LON, y=SRX$LAT)
hullpnt <- chull(act1)
hullpnt <- c(hullpnt,hullpnt[1])
poly <- act1[hullpnt,]
pip <- point.in.polygon(Xip,Yip,poly$x,poly$y)
hullDT <- data.frame(x=Xip,y=Yip,pip)
WR[which(hullDT$pip == 0), 3] <- NA
}
# コンター作成 ----------------------------------------------------
con.x <- unique(WR$LON)
con.y <- unique(WR$LAT)
con.z <- matrix(WR$PRED, ncol=length(con.x), byrow=T)
# 出力設定 ----------
xlength <- max(SRX$LON)-min(SRX$LON) # デバイスサイズの計算
ylength <- max(SRX$LAT)-min(SRX$LAT)
xratio <- xlength/xlength
yratio <- ylength/xlength
png("filledcontour_kirg.png", xratio*480*reso, yratio*480*reso)
par(xaxs = "i", yaxs = "i", mar=c(0,0,0,0), omi=c(0,0,0,0), bg="transparent", xaxt="n", yaxt="n")
# コンター作成 ------
image(con.x, con.y, con.z, col=rgb(cols(0:(LN-1)/(LN-1))/255), bg="transparent", bty="n")
maxn <- max(con.z[which(!is.na(con.z))])
minn <- min(con.z[which(!is.na(con.z))])
LV <- round(seq(minn, maxn, length=LN+1), dgt)
contour(con.x, con.y, con.z, levels=LV, add=T, labcex=reso)
points(SRX$LON, SRX$LAT, pch=16, col=4, cex=reso)
dev.off()
# KML 作成 ----------------------------------------------------------
GEfn <- paste(unlist(strsplit(fileName, ".csv")), "_等値線図.kml", sep="")
sink(GEfn)
cat('<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<GroundOverlay>
<name>Contour</name>
<color >5cffffff</color>
<Icon>
<href >filledcontour_kirg.png</href >
<viewRefreshTime >0</viewRefreshTime >
</Icon>
<altitude>300</altitude>
<altitudeMode >clampToGround </altitudeMode >
<LatLonBox id="contour">
<north>', max( SRX$LAT),'</north>
<south>', min( SRX$LAT),'</south>
<east>', max( SRX$LON),'</east>
<west>', min( SRX$LON),'</west>
</LatLonBox>
</GroundOverlay>
</kml >', sep='')
sink( )
