WRFの結果からRでHubenyの式を使ってポイントデータを抽出する

#WRFデータから緯度・経度を指定してデータを抜き出す

　WRFの計算結果からあるグリッド（または、指定する地点に最も近いグリッド）のデータを取得して、時系列のデータ解析を行いたい時、Rでデータを取得する方法をメモする。
　WRFの結果はグリッドの位置情報が緯度、経度で出力されるので、指定した位置と最も近い距離にあるグリッドを求めるよう、計算式にHubenyの公式を用いた。
　Hubenyの公式は、カシミール3Dなどの距離計算に使用されている簡易式で、計算式の概要は[二地点の緯度・経度からその距離を計算する](http://yamadarake.jp/trdi/report000001.html)[1]や、[ヒュベニの式が2種類ありますがどこが違うのですか？](http://www.amano-tec.com/apps/paceruler.html)[2]のページに説明されています。[1]にR版スクリプトが掲載されていましたので、一部改変して使用させていただきました。
　正直、ここまで厳密じゃなくても、緯度、経度差が一番小さくなるポイントを抜き出せばすむのだが、別件でHubenyの式を使ったので、使い方を記録するつもりで書きました。
***

##WRFの結果から必要な情報を取得する

```r
##  ライブラリの読み込み
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]
```

次元数を確認してみる

```r
> dim(LON)
[1] 94 74
```

　WRFの次元数は、Rでは[経度, 緯度, z軸レイヤ, 時間]として認識される。ここでは経度、緯度のみ取得していたため、2つしか出力されない。
　このデータでは、東西方向に94のポイントがあり、南北方向に74のポイントがある。

```r
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))`などとして求める。

```r
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)
}

```


##任意の地点がどのグリッドに対応するか確認する
```r
## x:経度, y:緯度, grid:グリッド間距離(km; namelistで設定)
Select.Point(x=130.376501, y=33.582399, grid=81)
## 例として福岡管区気象台の位置を設定
## グリッドは81km間隔で計算してる
```

実行すると、以下のような出力と画像がプロットされる。

```r
> 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）

![sample.png](https://qiita-image-store.s3.amazonaws.com/0/113843/53c0cc6f-0ca8-69d9-a6ef-a82654d714b5.png)
　青い菱型のポイントが、x, yで指定した福岡管区気象台の位置。赤い領域がWRFの結果で対応する領域。この赤い領域のデータを抽出する。
　データを抜き出すポイントがわかったので、あとは`ncvar_get(WRF, "T2")[ROW, COL, ]`などとして、WRF計算結果から任意の地点の時系列データを抽出する。

##任意の地点のデータを取得する
　抜き出す位置がわかったので、それを元にWRFデータから、地上2mの気温[℃]、地上2mの湿度[%]、地表面気圧[hPa]、海面更生気圧[hPa]、地上10m風速[m/s]、累積降水量[mm]、境界層高度[m]、標高[m]を取得する。

　WRFの結果は使いやすいよう、1日ごとに分割して出力するようにしているので、まずはファイルのリストを作成する。

```r
## 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ファイルに出力する。

```r
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の出力と同じであることを確認する。

```r
> 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形式で出力する。

```r
write.table(WRF.OutData, "./WRF_data.csv", sep=",", row.names=FALSE)
```

風向風速を角度で求めたいときは、`atan2()`関数を使って求めると良い

```r
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を足します。

```r
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)