#WRFデータから緯度・経度を指定してデータを抜き出す
WRFの計算結果からあるグリッド（または、指定する地点に最も近いグリッド）のデータを取得して、時系列のデータ解析を行いたい時、Rでデータを取得する方法をメモする。
WRFの結果はグリッドの位置情報が緯度、経度で出力されるので、指定した位置と最も近い距離にあるグリッドを求めるよう、計算式にHubenyの公式を用いた。
Hubenyの公式は、カシミール3Dなどの距離計算に使用されている簡易式で、計算式の概要は二地点の緯度・経度からその距離を計算する[1]や、ヒュベニの式が2種類ありますがどこが違うのですか？[2]のページに説明されています。[1]にR版スクリプトが掲載されていましたので、一部改変して使用させていただきました。
正直、ここまで厳密じゃなくても、緯度、経度差が一番小さくなるポイントを抜き出せばすむのだが、別件でHubenyの式を使ったので、使い方を記録するつもりで書きました。
##WRFの結果から必要な情報を取得する
## ライブラリの読み込み
library(ncdf4)
library(maps)
library(mapdata)
library(fields)
## WRFデータを読み込む
WRF <- nc_open("./wrfout_d01_2014-01-01_00:00:00")
## グリッド中心の緯度、経度を取得する
LON <- ncvar_get(WRF, "XLONG")[,,1]
LAT <- ncvar_get(WRF, "XLAT")[,,1]
次元数を確認してみる
> dim(LON)
[1] 94 74
WRFの次元数は、Rでは[経度, 緯度, z軸レイヤ, 時間]として認識される。ここでは経度、緯度のみ取得していたため、2つしか出力されない。
このデータでは、東西方向に94のポイントがあり、南北方向に74のポイントがある。
n.lon <- dim(LON)[1] ## WRFの東西グリッド数を取得
n.lat <- dim(LON)[2] ## WRFの南北グリッド数を取得
## 空白のラスターデータを作成
d.field <- matrix(nrow=n.lon, ncol=n.lat, 0)
##Hubenyの公式を用いて指定した地点に最も近いグリッドを計算する
簡単には、(x1-x2)+(y1-y2)が最小値になる位置を計算すれば良いのだが、グリッド中心位置との距離を求めるため、Hubenyの式を使った計算式を書きます。
最短距離から簡単に指定した地点を探す場合は、which.min関数を使って、which.min(abs(LON - 130.000) + abs(LAT - 32.000))などとして求める。
Select.Point <- function(x, y, grid){
## x:経度, y:緯度, grid:グリッド間距離(km; namelistで設定)
## 日本域のデータ取得を前提として、1°≒100kmとして距離を計算
## グリッド間距離を1/2する
Rx <- grid/100/2
a = 6378137.0
b = 6356752.314140
## xで指定した経度に近いデータに絞る
## x±Rxの範囲で絞る
LONGITUDE <- which(LON<=(x+Rx)&LON>=(x-Rx), arr.ind=TRUE)
Len.LON <- length(LONGITUDE[,1])
min.Distance <- grid
## LONGITUDEで絞った中から、更に距離が最も近いポイントを選ぶ
## for文内でHubenyの式を使用
for (i in 1:Len.LON){
row.data <- LONGITUDE[i,2]
col.data <- LONGITUDE[i,1]
x2 <- LON[col.data, row.data]
y2 <- LAT[col.data, row.data]
dy <- (y - y2)*pi/180
dx <- (x - x2)*pi/180
my <- ((y + y2) / 2)*pi/180
e2 <- (a^2 - b^2) / a^2
Mnum <- a * (1 - e2)
W <- sqrt(1 - e2 * sin(my)^2)
M <- Mnum / W^3
N <- a / W
Distance <- sqrt((dy * M)^2 + (dx * N * cos(my))^2)/10^3
if (Distance <= min.Distance){
min.Distance <- Distance
min.row <- row.data
min.col <- col.data
}
}
lon <- LON[min.col, min.row]
lat <- LAT[min.col, min.row]
## 求めたグリッド位置にデータを入力
d.field[min.col, min.row] <- 1
## データ取得位置と入力した位置情報をマップ上にプロット
image.plot(LON, LAT, d.field, col=c(NA, "red"), border=TRUE,
main="データ抽出地点",xlab="Longitude(deg)",ylab="Latitude(deg)",
xlim=c(lon-grid/20, lon+grid/20), ylim=c(lat-grid/20, lat+grid/20))
map("worldHires", add=T)
map("japan", add=T)
points(x, y, col="blue", pch=18, cex=2)
legend("topleft", legend=c("MODEL", "REAL"), col=c("red", "blue"), pch=c(15,18), bg = "pink", cex=1.2)
text(lon, lat + grid/50, x, col="blue", cex=1.5)
text(lon, lat + grid/100, y, col="blue", cex=1.5)
text(lon, lat - grid/100, lon, col="red", cex=1.5)
text(lon, lat - grid/50, lat, col="red", cex=1.5)
grid(col="black")
select.point <- data.frame(ROW = min.row,
COL = min.col,
Latitude = lat,
Longitude = lon,
Distance=min.Distance)
return(select.point)
}
##任意の地点がどのグリッドに対応するか確認する
## x:経度, y:緯度, grid:グリッド間距離(km; namelistで設定)
Select.Point(x=130.376501, y=33.582399, grid=81)
## 例として福岡管区気象台の位置を設定
## グリッドは81km間隔で計算してる
実行すると、以下のような出力と画像がプロットされる。
> Select.Point(x=130.376501, y=33.582399, grid=81)
ROW COL Latitude Longitude Distance
col 44 65 33.51446 130.601 22.16525
ROW：xに最も近いグリッド経度地点
COL：yに最も近いグリッド緯度地点
Latitude：ROWに対応するグリッド値（緯度）
Longitude：COLに対応するグリッド値（経度）
Distance：任意の地点とグリッド中心位置との距離（km）
青い菱型のポイントが、x, yで指定した福岡管区気象台の位置。赤い領域がWRFの結果で対応する領域。この赤い領域のデータを抽出する。
データを抜き出すポイントがわかったので、あとはncvar_get(WRF, "T2")[ROW, COL, ]などとして、WRF計算結果から任意の地点の時系列データを抽出する。
##任意の地点のデータを取得する
抜き出す位置がわかったので、それを元にWRFデータから、地上2mの気温[℃]、地上2mの湿度[%]、地表面気圧[hPa]、海面更生気圧[hPa]、地上10m風速[m/s]、累積降水量[mm]、境界層高度[m]、標高[m]を取得する。
WRFの結果は使いやすいよう、1日ごとに分割して出力するようにしているので、まずはファイルのリストを作成する。
## Select.Pointで求めたCOL, ROWを設定する
COL <- 65
ROW <- 44
## WRFファイルのリストを作成
## データを抽出するファイルだけまとめてフォルダを作っておくと楽
wrf.dir <- "./WRF_outdata/" ## 最後に必ずスラッシュ"/"をつけること
files <- list.files(path=wrf.dir, pattern = "wrfout_d01")
wrf.len <- length(files)
注意：WRFのデータは、経度、緯度、z軸レイヤ、時間の4次元データを基準にしているため、1時間しかデータが格納されていないNetCDFファイルだとエラーが出る。そんなファイルは除いてください。
WRFファイル数だけ処理を繰り返し、まとめてcsvファイルに出力する。
WRF.OutData <- data.frame()
for (i in 1:wrf.len){
WRF <- nc_open(paste(wrf.dir, files[i], sep=""))
TIME <- ncvar_get(WRF, "Times")
## JSTに変換
str.T <- as.POSIXlt(paste(substring(TIME[1], 1, 10), "00:00:00", sep=" "), "%Y-%m-%d %H:%M:%OS") + 9*60*60
JST <- seq(str.T, by="hour", length=24)
TEMP <- ncvar_get(WRF, "T2")[COL, ROW, ]-273.15 #地上2mの気温
RAIN <- ncvar_get(WRF, "RAINC")[COL, ROW, ] #累積降水量(mm)
PRES <- ncvar_get(WRF, "PSFC")[COL, ROW, ]/100 #地表付近の気圧(hPa)
HGT <- ncvar_get(WRF, "HGT")[COL, ROW, 1] ## 地表面高さ[m]
S.PRES <- PRES*(1-0.0065*HGT/(TEMP+273.15+0.0065*HGT))^(-5.257) ## 海面更正気圧[hPa]
XWS <- ncvar_get(WRF, "U10")[COL, ROW, ] #地上10mのx軸風速(西→東)
YWS <- ncvar_get(WRF, "V10")[COL, ROW, ] #地上10mのy軸風速(南→北)
Lon <- ncvar_get(WRF, "XLONG")[COL, ROW, ] #経度
Lat <- ncvar_get(WRF, "XLAT")[COL, ROW, ] #緯度
QV <- ncvar_get(WRF, "Q2")[COL, ROW, ] #地上2mの混合比(水蒸気質量kg/乾燥空気質量kg)
## 飽和水蒸気圧の計算。Tetensの式だが、Goff-Gratchの式と近似できる。
PHUM <- 6.1078*10^(
7.5*TEMP/(TEMP+237.3)
)
RH <- 100*PRES*QV/(PHUM*(0.622+QV))
PBLH <- ncvar_get(WRF, "PBLH")[COL, ROW, ] #境界層高度
Mid.frame <- data.frame(TIME_UTC=TIME,
TIME_JST=JST,
LON_deg=Lon,
LAT_deg=Lat,
TEMP_degC=TEMP,
RH_percent=RH,
RAIN_mm=RAIN,
SFC_PRES_hPa=PRES,
SEA_LEVEL_PRES_hPa=S.PRES,
U10=XWS,
V10=YWS,
PBLH_m=PBLH,
HEIGHT_m=HGT)
WRF.OutData<- rbind(WRF.OutData, Mid.frame)
}
head(WRF.OutData)で緯度(LAT)、経度(LON)がSelect.Pointの出力と同じであることを確認する。
> head(WRF.OutData)
TIME_UTC TIME_JST LON_deg LAT_deg TEMP_degC RH_percent
1 2014-01-15_00:00:00 2014-01-15 09:00:00 130.601 33.51446 3.787439 60.37696
2 2014-01-15_01:00:00 2014-01-15 10:00:00 130.601 33.51446 5.069299 62.86162
3 2014-01-15_02:00:00 2014-01-15 11:00:00 130.601 33.51446 5.671594 60.55730
4 2014-01-15_03:00:00 2014-01-15 12:00:00 130.601 33.51446 6.196771 60.33883
5 2014-01-15_04:00:00 2014-01-15 13:00:00 130.601 33.51446 6.581842 61.80480
6 2014-01-15_05:00:00 2014-01-15 14:00:00 130.601 33.51446 6.771753 64.21387
RAIN_mm SFC_PRES_hPa SEA_LEVEL_PRES_hPa U10 V10 PBLH_m HEIGHT_m
1 0 1005.317 1028.986 1.713518 -1.446464 0.000 189.0172
2 0 1005.254 1028.811 2.567818 -4.762895 1257.336 189.0172
3 0 1004.958 1028.457 2.458457 -5.227210 1258.585 189.0172
4 0 1004.848 1028.300 2.323239 -5.595024 1260.606 189.0172
5 0 1004.771 1028.188 2.463428 -5.933176 1262.666 189.0172
6 0 1004.595 1027.992 2.521139 -6.201655 1264.699 189.0172
抽出したデータをcsv形式で出力する。
write.table(WRF.OutData, "./WRF_data.csv", sep=",", row.names=FALSE)
風向風速を角度で求めたいときは、atan2()関数を使って求めると良い
x <- c(1, sqrt(3)/2, 1/2, 0, -1/2, -sqrt(3)/2, -1, -sqrt(3)/2, -1/2, 0, 1/2, sqrt(3)/2, 1)
y <- c(0, 1/2, sqrt(3)/2, 1, sqrt(3)/2, 1/2, 0, -1/2, -sqrt(3)/2, -1, -sqrt(3)/2, -1/2, 0)
> atan2(y, x)/pi*180
[1] 0 30 60 90 120 150 180 -150 -120 -90 -60 -30 0
東が0°で、反時計回りに角度が決まる、三角関数で風向が出力される。ただし、180°以上ではマイナス値になるため、その場合は360を足します。
WD <- atan2(y, x)/pi*180
> WD
[1] 0 30 60 90 120 150 180 -150 -120 -90 -60 -30 0
WD <- ifelse(WD<0, WD+360, WD)
> WD
[1] 0 30 60 90 120 150 180 210 240 270 300 330 0
風速はsqrt(x^2 + y^2)で求まりますね。
抽出したデータはopenairパッケージを使うと綺麗なグラフが出来上がります。
参考：
[1] 日本は山だらけ〜 技術研究本部 報告書 第一号 (2009).（http://yamadarake.jp/trdi/report000001.html)
[2] アマノ技研. (http://www.amano-tec.com/apps/paceruler.html)
[3] 師匠の散歩 VBAと測地学 Hubenyの式を考察. (http://tancro.e-central.tv/grandmaster/excel/hubenystandard.html)