データは JODC の 500m メッシュ水深を利用させていただこうと思います。
早速ですが、データの取得です。
取得方法は JODC の説明を見てもらえば分かると思いますが、簡単に説明します。
[ダウンロードするメッシュを選択] にチェックを入れて、選択し、[データ取得] をクリックします。
すると、zip ファイルがダウンロードされるので、それを解凍して下さい。
zip の中には [jodc-depth500mesh-...] という名前の txt ファイルが入っていると思います。
今回はこれをこのまま使います。
次に、以下のコードを実行し、海底地形を可視化します。
完成イメージはこんな感じです。
今回は沖縄周辺のデータを使いました。
プログラムについて ------------------------------------------
JODC の 500m メッシュ水深データは [緯度・経度・水深] というセットになっていて、厳密には格子状ではありません。
これを R が扱いやすい格子状にするために、直角平面座標系に変換して、格子状に整形します。
直角平面座標系への変換については測量計算サイトを参考にしました。
格子状に整形した後、平滑化を行って、パッケージ rgl で可視化します。
#作成日時 : 18:56 2011/08/05
#作成者名：T.Nakano
rm(list=ls())
gc();gc()
library(rgl)
library(sp)
# ///////////////////////////// JODC 500m メッシュ水深データの可視化 ////////////////////////////////////////////////////////////////////////////////////////////
# 整形するメッシュサイズ
ms <- 500 # 単位：m
# 子午線長 S を求める関数 ---------------------------------------------------------------------------------------------------------------------------
fun_calS <- function(φ){
a <- 6378137
f <- 1/298.257222101
e <- sqrt(2*f-f^2)
rad <- pi/180
A <- 1+(3/4*e^2)+(45/64*e^4)+(175/256*e^6)+(11025/16384*e^8)+
(43659/65536*e^10)+(693693/1048576*e^12)+
(19324305/29360128*e^14)+(4927697775/7516192768*e^16)
B <- (3/4*e^2)+(15/16*e^4)+(525/512*e^6)+(2205/2048*e^8)+(72765/65536*e^10)+
(297297/262144*e^12)+(135270135/117440512*e^14)+(547521975/469762048*e^16)
C <- (15/64*e^4)+(105/256*e^6)+(2205/4096*e^8)+(10395/16384*e^10)+
(1486485/2097152*e^12)+(45090045/58720256*e^14)+(766530765/939524096*e^16)
D <- (35/512*e^6)+(315/2048*e^8)+(31185/131072*e^10)+(165165/524288*e^12)+
(45090045/117440512*e^14)+(209053845/469762048*e^16)
E <- (315/16384*e^8)+(3465/65536*e^10)+(99099/1048576*e^12)+
(4099095/29360128*e^14)+(348423075/1879048192*e^16)
F <- (693/131072*e^10)+(9009/524288*e^12)+(4099095/117440512*e^14)+
(26801775/469762048*e^16)
G <- (3003/2097152*e^12)+(315315/58720256*e^14)+(11486475/939524096*e^16)
H <- (45045/117440512*e^14)+(765765/469762048*e^16)
I <- (765765/7516192768*e^16)
B1 <- a*(1-e^2)*A
B2 <- a*(1-e^2)*(-B/2)
B3 <- a*(1-e^2)*(C/4)
B4 <- a*(1-e^2)*(-D/6)
B5 <- a*(1-e^2)*(E/8)
B6 <- a*(1-e^2)*(-F/10)
B7 <- a*(1-e^2)*(G/12)
B8 <- a*(1-e^2)*(-H/14)
B9 <- a*(1-e^2)*(I/16)
S <- B1*(φ*rad)+B2*sin(2*φ*rad)+B3*sin(4*φ*rad)+B4*sin(6*φ*rad)+
B5*sin(8*φ*rad)+B6*sin(10*φ*rad)+B7*sin(12*φ*rad)+
B8*sin(14*φ*rad)+B9*sin(16*φ*rad)
return(S)
}
# ----------------------------------------------------------------------------
# 地理座標から直交座標に変換する関数 ---------------------------------------------------------------------------------------------------------------
fun_GC2RC <- function(φ, λ, No=0){ # 変換したい緯度(φ), 経度(λ), 当該都道府県の系番号(No) を引数に与える
# ベクトル型で与えることができる
if(No!=0){
oo <- data.frame(Range=1:19,
LON=c(129.5, 131, 132.166666, 133.5, 134.333333, 136, 137.166666,
138.5, 139.833333, 140.8333333, 140.25, 142.25, 144.25,
142, 127.5, 124, 131, 136, 154),
LAT=c(33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44,
44, 26, 26, 26, 26, 20, 26))
φ0 <- oo$LAT[which(oo$Range==No)]
λ0 <- oo$LON[which(oo$Range==No)]
}else if(No==0){
φ0 <- min(φ)
λ0 <- min(λ)
}
Δλ <- λ-λ0
a <- 6378137
f <- 1/298.257222101
e <- sqrt(2*f-f^2)
ed <- sqrt(2/f-1)/((1/f)-1)
rad <- pi/180
η2 <- ed^2*cos(φ*rad)^2
t <- tan(φ*rad)
m0 <- 0.9999
W <- sqrt(1-e^2*sin(φ*rad)^2)
N <- a/W
M <- a*(1-e^2)/W^3
S <- fun_calS(φ)
S0 <- fun_calS(φ0)
x <- ((S-S0)+(1/2)*N*cos(φ*rad)^2*t*(Δλ*rad)^2+
(1/24)*N*cos(φ*rad)^4*t*(5-t^2+9*η2+4*η2^2)*(Δλ*rad)^4-
(1/720)*N*cos(φ*rad)^6*t*(-61+58*t^2-t^4-270*η2+330*t^2*η2)*(Δλ*rad)^6-
(1/40320)*N*cos(φ*rad)^8*t*(-1385+3111*t^2-543*t^4+t^6)*(Δλ*rad)^8)*m0
y <- (N*cos(φ*rad)*(Δλ*rad)-
(1/6)*N*cos(φ*rad)^3*(-1+t^2-η2)*(Δλ*rad)^3-
(1/120)*N*cos(φ*rad)^5*(-5+18*t^2-t^4-14*η2+58*t^2*η2)*(Δλ*rad)^5-
(1/5040)*N*cos(φ*rad)^7*(-61+479*t^2-179*t^4+t^6)*(Δλ*rad)^7)*m0
m <- m0*(1+(y^2)/(2*M*N*m0^2)+(y^4)/(24*M^2*N^2*m0^4))
y <- y*m
return(data.frame(X=x, Y=y))
}
# 直交座標から地理座標に変換する関数 ------------------------------------------------------------------------------------------------------------------------
fun_RC2GCφ1 <- function(x, y, No=0){ # 変換したい x 座標, y 座標, 任意の原点 (φ0, λ0), 当該都道府県の系番号(No) を引数に与える
# No=0 のときは φ0, λ0 を必ず入力すること
# ベクトル型で与えることができる
if(No!=0){
oo <- data.frame(Range=1:19,
LON=c(129.5, 131, 132.166666, 133.5, 134.333333, 136, 137.166666,
138.5, 139.833333, 140.8333333, 140.25, 142.25, 144.25,
142, 127.5, 124, 131, 136, 154),
LAT=c(33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44,
44, 26, 26, 26, 26, 20, 26))
φ0 <- oo$LAT[which(oo$Range==No)]
λ0 <- oo$LON[which(oo$Range==No)]
}else if(No==0){
φ0 <- min(φ)
λ0 <- min(λ)
}
a <- 6378137
f <- 1/298.257222101
e <- sqrt(2*f-f^2)
ed <- sqrt(2/f-1)/((1/f)-1)
m0 <- 0.9999
rad <- pi/180
actφn1 <- rep(NA,length(x))
for(i in 1:length(x)){
φn <- φn1 <- M <- Sφn <- NULL
φn <- φ0
M <- fun_calS(φ0)+(x[i]/m0)
Sφn <- fun_calS(φn)
φn1 <- φn*rad+(2*(Sφn-M)*(1-e^2*sin(φn*rad)^2)^(3/2))/(3*e^2*(Sφn-M)*sin(φn*rad)*cos(φn*rad)*(1-e^2*sin(φn*rad)^2)^(1/2)-2*a*(1-e^2))
while(abs(φn1-φn) >= 2e-10){
φn <- φn1
Sφn <- fun_calS(φn)
φn1 <- φn+(2*(Sφn-M)*(1-e^2*sin(φn)^2)^(3/2))/(3*e^2*(Sφn-M)*sin(φn)*cos(φn)*(1-e^2*sin(φn)^2)^(1/2)-2*a*(1-e^2))
}
actφn1[i] <- φn1
cat(i, " ", length(x), "\n")
}
φ1 <- actφn1*rad
return(φ1)
}
# ------------------------------------------------------------------------------------------------------------------
fun_RC2GC <- function(φ1, y, No=0){
if(No!=0){
oo <- data.frame(Range=1:19,
LON=c(129.5, 131, 132.166666, 133.5, 134.333333, 136, 137.166666,
138.5, 139.833333, 140.8333333, 140.25, 142.25, 144.25,
142, 127.5, 124, 131, 136, 154),
LAT=c(33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44,
44, 26, 26, 26, 26, 20, 26))
φ0 <- oo$LAT[which(oo$Range==No)]
λ0 <- oo$LON[which(oo$Range==No)]
}else if(No==0){
φ0 <- min(φ)
λ0 <- min(λ)
}
a <- 6378137
f <- 1/298.257222101
e <- sqrt(2*f-f^2)
ed <- sqrt(2/f-1)/((1/f)-1)
m0 <- 0.9999
rad <- pi/180
η12 <- ed^2*cos(φ1)^2
t1 <- tan(φ1)
W <- sqrt(1-e^2*sin(φ1)^2)
N1 <- a/W
φ <- φ1-(1/2)*(1/N1^2)*t1*(1+η12)*(y/m0)^2+
(1/24)*(1/N1^4)*t1*(5+3*t1^2+6*η12-6*t1^2*η12-3*η12^2-9*t1^2*η12^2)*(y/m0)^4-
(1/720)*(1/N1^6)*t1*(61+90*t1^2+45*t1^4+107*η12-162*t1^2*η12-45*t1^4*η12)*(y/m0)^6+
(1/40320)*(1/N1^8)*t1*(1385+3633*t1^2+4095*t1^4+1575*t1^6)*(y/m0)^8
φ <- φ/rad
Δλ <- (1/(N1*cos(φ1)))*(y/m0)-
(1/6)*(1/(N1^3*cos(φ1)))*(1+2*t1^2+η12)*(y/m0)^3+
(1/120)*(1/(N1^5*cos(φ1)))*(5+28*t1^2+24*t1^4+6*η12+8*t1^2*η12)*(y/m0)^5-
(1/5040)*(1/(N1^7*cos(φ1)))*(61+622*t1^2+1320*t1^4+720*t1^6)*(y/m0)^7
Δλ <- Δλ/rad
W <- sqrt(1-e^2*sin(φ*rad)^2)
N <- a/W
M <- a*(1-e^2)/W^3
m <- m0*(1+(y^2)/(2*M*N*m0^2)+(y^4)/(24*M^2*N^2*m0^4))
λ <- λ0+(Δλ)/m
return(data.frame(緯度=φ, 経度=λ))
}
# 海底地形を平滑化する関数 ---------------------------------------------------------------------------------------
fun_SmoothBF <- function(depth){ # 海底地形の平滑化
# depth はマトリクス型の水深データ (陸地は0)
depth0 <- as.vector(depth) # データとその
depth1 <- as.vector(rbind(depth[-1,],rep(NA,ncol(depth)))) # 南側
depth2 <- as.vector(cbind(rep(NA,nrow(depth)),depth[,-ncol(depth)])) # 西側
depth3 <- as.vector(rbind(rep(NA,ncol(depth)),depth[-nrow(depth),])) # 北側
depth4 <- as.vector(cbind(depth[,-1],rep(NA,nrow(depth)))) # 東側のデータを
depthA <- cbind(depth0,depth1,depth2,depth3,depth4) # cbind で 5 列のデータにする。
depthcond <- depth0
depthcond[which(depthcond!=0)] <- 1
depthIP <- apply(depthA,1,function(da) mean(da[!is.na(da)]))*depthcond # 行方向で NA 以外の平均を算出し
depth <- matrix(depthIP,nrow(depth),ncol(depth)) # マトリクス型に変換して depth に返す
return(depth)
}
# 格子にしたことにより生じた幻の島を海底に沈めるする関数 ---------------------------------------------------------------------------------------
fun_sinking <- function(depth){ # 海底地形の平滑化とほとんど同じ
# depth はマトリクス型の水深データ (陸地は0)
depth0 <- as.vector(depth) # データとその
depth1 <- as.vector(rbind(depth[-1,],rep(NA,ncol(depth)))) # 南側
depth2 <- as.vector(cbind(rep(NA,nrow(depth)),depth[,-ncol(depth)])) # 西側
depth3 <- as.vector(rbind(rep(NA,ncol(depth)),depth[-nrow(depth),])) # 北側
depth4 <- as.vector(cbind(depth[,-1],rep(NA,nrow(depth)))) # 東側のデータを
depthA <- cbind(depth1,depth2,depth3,depth4) # cbind で 5 列のデータにする。
depthcond <- depth0
depthcond[which(depthcond==0)] <- NA
depthcond[which(!is.na(depthcond))] <- 0
depthcond[which(is.na(depthcond))] <- 1
depthIP <- apply(depthA,1,function(da) mean(da[!is.na(da)]))*depthcond # 行方向で NA 以外の平均を算出し
depthIP <- depthIP+depth0
depth <- matrix(depthIP,nrow(depth),ncol(depth)) # マトリクス型に変換して depth に返す
return(depth)
}
# ------------------------------------------------------------------------------------------------------------------------------------------------
# データの読込み ------------------------------------------------------------------------------------------
library(tcltk)
cat("JODC から取得したファイルを選択する","\n")
fileName <- tclvalue(tkgetOpenFile())
act1 <- unlist(strsplit(fileName,"/"))
filename2 <- act1[length(act1)]
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)
if(length(dir(pattern=paste(filename2,".xdr",sep="")))==0){
cat("データを格子状に整列する (約 10 分)","\n")
topoDT <- read.table(fileName, header=FALSE)
topoDT <- topoDT[,-1]
colnames(topoDT) <- c("LAT","LON","H")
φ <- topoDT$LAT
λ <- topoDT$LON
GCdt <- fun_GC2RC(φ, λ)
GCdt <- data.frame(GCdt,H=topoDT$H)
DT <- GCdt
plot(DT[,2:1], pch=16, cex=0.5)
# 格子状に整形 ----------------------------------
# このとき、理論上最大250mのズレが生じる可能性があります。
# 悪しからずご了承ください (ただし、今のところそのような状態には陥っていません)。
DT <- DT/ms
DT$X <- round(DT$X,0)
DT$Y <- round(DT$Y,0)
DT <- DT*ms
DT <- DT[rev(order(DT$X)),]
DT <- DT[order(DT$Y),]
plot(DT[,2:1], pch=16, cex=0.5)
# 同じ地点のデータを 1 つにする ------------------------------------------------------------------
actX <- unique(DT$X)
actX <- actX[rev(order(actX))]
actY <- unique(DT$Y)
actY <- actY[order(actY)]
MeshX <- rep(actX, length(actY))
MeshY <- as.vector(sapply(actY, function(da) rep(da, length(actX))))
MeshDT <- data.frame(No=1:length(MeshX), X=MeshX, Y=MeshY, H=NA)
MeshDT$H <- mapply(function(da,db,dc) {st <- Sys.time()
act <- mean(DT$H[which(DT$X==da & DT$Y==db)])
ed <- Sys.time()
endtime <- paste(Sys.time()+(as.numeric(difftime(ed,st,units="sec"))+0.01)*(length(MeshX)-dc))
cat(round(dc/length(MeshX)*100,2), "% ", endtime,"\n")
return(act)}, MeshDT$X, MeshDT$Y, MeshDT$No)
MeshDT$H[which(is.na(MeshDT$H))] <- 0
topoH <- matrix(MeshDT$H, length(actX), length(actY))
persp3d(t(-topoH[nrow(topoH):1,]),col="lightblue",aspect=c(1,1,0.5))
topoH <- fun_sinking(topoH)
topoH[1,] <- topoH[2,]
topoH[nrow(topoH),] <- topoH[nrow(topoH)-1,]
topoH[,1] <- topoH[,2]
topoH[,ncol(topoH)] <- topoH[,ncol(topoH)-1]
persp3d(t(-topoH[nrow(topoH):1,]),col="lightblue",aspect=c(1,1,0.5))
# 平滑化 -----------------------------------------------------------------------------
for(i in 1:1){
topoH <- fun_SmoothBF(topoH)
persp3d(t(-topoH[nrow(topoH):1,]),col="lightblue",aspect=c(1,1,0.5))
}
MeshDT$H <- as.vector(topoH)
uniX <- unique(MeshDT$X)
GCφ1dt <- fun_RC2GCφ1(uniX, MeshDT$Y)
φ1DTfm <- data.frame(X=uniX, φ1=GCφ1dt)
MeshDT <- data.frame(ID=1:nrow(MeshDT), MeshDT)
MeshDT <- merge(MeshDT,φ1DTfm)
MeshDT <- MeshDT[order(MeshDT$ID),]
GCdtfm <- fun_RC2GC(MeshDT$φ1, MeshDT$Y)
MeshDT <- data.frame(GCdtfm, H=MeshDT$H)
save(MeshDT, file=paste(fileName,".xdr",sep=""))
} # end if
# xdr がある場合はここからスタート ---------------------------------------------------------------------------------------------------------
topoDT <- read.table(fileName, header=FALSE)
topoDT <- topoDT[,-1]
colnames(topoDT) <- c("LAT","LON","H")
φ <- topoDT$LAT
λ <- topoDT$LON
GCdt <- fun_GC2RC(φ, λ)
GCdt <- data.frame(GCdt,H=topoDT$H)
DT <- GCdt
plot(DT[,2:1], pch=16, cex=0.5, asp=T)
# 格子状に整形 ----------------------------------
DT <- DT/ms
DT$X <- round(DT$X,0)
DT$Y <- round(DT$Y,0)
DT <- DT*ms
DT <- DT[rev(order(DT$X)),]
DT <- DT[order(DT$Y),]
plot(DT[,2:1], pch=16, cex=0.5, asp=T)
actX <- unique(DT$X)
actX <- actX[rev(order(actX))]
actY <- unique(DT$Y)
actY <- actY[order(actY)]
load(paste(filename2,".xdr",sep="")) # xdr の読込み
X <- rep(actX,length(actY))
Y <- as.vector(unlist(sapply(actY, function(da) rep(da,length(actX)))))
MeshDT <- data.frame(MeshDT,X=X,Y=Y)
topoH <- matrix(MeshDT$H, length(actX), length(actY))
#persp3d(t(-topoH[nrow(topoH):1,]),col="lightblue",aspect=c(1,1,0.5))
#contour(t(-topoH[nrow(topoH):1,]))
#act <- edit(MeshDT)
# 抽出範囲の選択用のコンターマップ作成 --------------------------------------------------------------
#pdf("ContourMap.pdf", family="Japan1GothicBBB")
XN <- ncol(topoH)
YN <- nrow(topoH)
LX <- seq(min(MeshDT$Y), max(MeshDT$Y), (max(MeshDT$Y)-min(MeshDT$Y))/(XN-1))
LY <- seq(min(MeshDT$X), max(MeshDT$X), (max(MeshDT$X)-min(MeshDT$X))/(YN-1))
cols = colorRamp(c("white", "royalblue1"))
n=length(seq(0,max(topoH), length=10))
image(LX, LY,t(topoH[nrow(topoH):1,]), col=rgb(cols(0:(n-1)/(n-1))/255), xlab="", ylab="", bty="n", xaxt="n" , yaxt="n", asp=T)
contour(LX, LY,t(topoH[nrow(topoH):1,]), add=TRUE, levels=seq(0,max(topoH), length=n+1))
topoH2 <- topoH
topoH2[which(topoH2 > 0)] <- NA
image(LX, LY,t(-topoH2[nrow(topoH2):1,]),col=1, xlab="", ylab="",bg="transparent",axes=F, asp=T, add=T)
#dev.off()
# 抽出範囲の選択 ---------------------------------------------------------------------------------------
cat("抽出範囲の選択","\n","左下端 → 右上端 の順でクリックする","\n")
for(i in 1:100){
range <- locator(2)
polygon(c(range$x[1],range$x[2],range$x[2],range$x[1],range$x[1]),
c(range$y[1],range$y[1],range$y[2],range$y[2],range$y[1]),col=i+1, density=0)
cond <- tkmessageBox(message="抽出範囲を正しく選択できましたか？",icon="question",type="yesno",default="yes")
if(as.character(cond)=="yes") break
}
poly <- data.frame(x=c(range$x[1],range$x[2],range$x[2],range$x[1],range$x[1]), y= c(range$y[1],range$y[1],range$y[2],range$y[2],range$y[1]))
pip <- point.in.polygon(MeshDT$Y,MeshDT$X, poly$x, poly$y)
MeshDT <- MeshDT[which(pip==1),]
actX <- unique(MeshDT$X)
actX <- actX[rev(order(actX))]
actY <- unique(MeshDT$Y)
actY <- actY[order(actY)]
topoH <- matrix(MeshDT$H, length(actX), length(actY))
XN <- ncol(topoH)
YN <- nrow(topoH)
LX <- seq(min(MeshDT$経度), max(MeshDT$経度), (max(MeshDT$経度)-min(MeshDT$経度))/(XN-1))
LY <- seq(min(MeshDT$緯度), max(MeshDT$緯度), (max(MeshDT$緯度)-min(MeshDT$緯度))/(YN-1))
ratio1 <- (max(LX)-min(LX))/(max(LY)-min(LY))
topoH2 <- topoH
topoH2[which(topoH2!=0)] <- NA
topoH2[which(topoH2==0)] <- 0.5
persp3d(LX, LY, -t(topoH[nrow(topoH):1,]), col="lightblue", aspect=c(ratio1, 1, 0.5), xlab="", ylab="", zlab="", axes=F, box=FALSE)
persp3d(LX, LY, t(topoH2[nrow(topoH2):1,]), col="darkgreen", add=T)
# コンターマップ作成 --------------------------------------------------------------
cols = colorRamp(c("white", "royalblue1"))
n=length(seq(0,max(topoH), length=10))
image(LX, LY,t(topoH[nrow(topoH):1,]), col=rgb(cols(0:(n-1)/(n-1))/255), xlab="", ylab="", bty="n", xaxt="n" , yaxt="n", asp=T)
contour(LX, LY,t(topoH[nrow(topoH):1,]), add=TRUE, levels=seq(0,max(topoH), length=n+1))
topoH2 <- topoH
topoH2[which(topoH2 > 0)] <- NA
image(LX, LY,t(-topoH2[nrow(topoH2):1,]),col=1, xlab="", ylab="",bg="transparent",axes=F, add=T)
使い方 ------------------------------------------
プログラムを実行すると、ファイル選択ウィンドウが開くので、解答した txt データを選択して下さい。
すると、格子状への整形が始まります。
これは少々時間がかかります。
JODC でダウンロードする際に選択した区画一つにつき、だいたい 10 分程度かかります (もちろん PC のスペックにもよります)。
整形作業の進捗をパーセンテージと終了時刻で表示しています。
こんな感じです。
整形が完了すると xdr ファイルで保存されます。
これで次から整形作業を待つ必要がなくなります。
次に、ダウンロードしたデータの範囲で等値線図が描かれます。
こんな感じです。
ここで 3D で可視化したい範囲を [左下 → 右上] の順にクリックします。
こんな感じです。
すると、この範囲で海底地形が 3D で可視化されます。
これは rgl で描かれているので、マウスで図をドラッグするとグリグリ動かすことができます。
ここで、さらにもう一つ、
本コードに続けて以下のコードを実行すると、海底地形図内にプロットを描くことも可能です。
plot3d( その他には、persp3d() という関数を使うと等値面を描くことも可能です。
例えば、CTD とかで大量にデータを取った場合なんかは、水温等の等値面を描くことができます。
データの可視化という意味では色々と使えそうだなと思います。
