RでWRFのデータを地図上にプロット

#WRFの結果をggplot2を使ってマップ上に描画
　WRFとCMAQを使ったとき、投影法をlccにすると計算結果を地図上に描画するときに苦労したのでメモを残しておく。
　普段はGMTを使っているが、データの解析をしつつ、コンター図を描画したい時にRを使うと便利なので、今回はデータの読み込みからマッププロットまですべてRを使って行う。
　使用したRのバージョンはR-3.2.3で、気温などのカラーイメージ、等圧線、風向風速の3つを描画できるようにする。
　ちなみに、ポイントのデータをきちんと見たい場合以外（なんとなく全体の空間分布が把握できれば良い）は、filled.contour()、contour()、fieldsパッケージのimage.plot()、rasterVisパッケージのlevelplot()などの関数を使うほうが早い。

##  WRFの結果はNetCDF形式であるので、ncdf4パッケージを使用してデータの抽出を行う。
##  必要なパッケージの読み込み 
library(ncdf4)
library(maps)
library(mapdata)
library(ggplot2)
library(mapproj)
#library(grid)
library(colorRamps)
library(raster)

###NetCDF形式のファイルから必要なデータを読みだす。

## WRFファイルの読み込み
nc <- nc_open("~/wrfout_d01_2014-01-01_00:00:00")
t <- 1 ## 抜き出すデータの時間

## 位置データの抜き出し
LON <- ncvar_get(nc, "XLONG")[,,t]
LAT <- ncvar_get(nc, "XLAT")[,,t]

## 地上2mの気温、最下層の気圧、地上10mの風速成分を取り出す
TEMP <- ncvar_get(nc, "T2")[,,t]-273.15 ## [K] → [℃]
HGT <- ncvar_get(nc, "HGT")[,,time] ## 地表面高さ[m]
PSFC <- ncvar_get(nc, "PSFC")[,,time] ## [Pa]
PRES <- PSFC*(1-0.0065*HGT/(TEMP+273.15+0.0065*HGT))^(-5.257)/100  ## 海面更正気圧[hPa]
U <- ncvar_get(nc, "U10")[,,t] ## 西→東方向の風速
V <- ncvar_get(nc, "V10")[,,t] ## 南→北方向の風速
WRF <- data.frame(LON=matrix(LON), LAT=matrix(LAT), conc=matrix(TEMP), Col=matrix(TEMP), pres=matrix(PRES))

###コンター図作成のため、気圧データをラスタライズする
　地図投影法がlccだとgeom_contourでコンター図ができないため、気圧データをラスタライズする

##  raster()で形を作り、extent()で範囲を設定する
##  参考値で100列、100行のラスターデータを作成するが、ここはWRFの結果に合わせて変更する
ColNum=100
RowNum=100
r <- raster(ncols=ColNum, nrows=RowNum)  ## 空白のらラスタ―データ作成
extent(r) <- c(min(LON), max(LON), min(LAT), max(LAT))  ## データ範囲の設定
Pres <- rasterize(WRF[1:2], r, pres)  ## 気圧をz軸データとしたラスターデータを作成

##  作成したラスタデータから気圧データを抽出する
xval <- matrix(matrix(nrow=RowNum, ncol=ColNum, seq(min(LON), max(LAT), length=RowNum)))
yval <- matrix(t(matrix(nrow=RowNum, ncol=ColNum, rev(seq(min(LAT), max(LAT), length=ColNum)))))
ext <- extent(min(LON), max(LON), min(LAT), max(LAT))
presVal <- extract(Pres, ext)
r.PRES <- data.frame(xval, yval, presVal)

###風向・風速データの準備

##  そのままではデータ数が多すぎるため、適当に端折る
cutLev <- 4  ##  4つ置きのデータを作成する
x <- dim(U)[1]
y <- dim(U)[2]
wlon <- LON[seq(1, x, by=cutLev),seq(1, y, by=cutLev)]
wlat <- LAT[seq(1, x, by=cutLev),seq(1, y, by=cutLev)]
wu <- U[seq(1, x, by=cutLev),seq(1, y, by=cutLev)]
wv <- V[seq(1, x, by=cutLev),seq(1, y, by=cutLev)]
wind <- data.frame(LON=matrix(wlon), LAT=matrix(wlat), U=matrix(wu), V=matrix(wv))

###ggplot2を使ってデータを描画

##  何度も条件を変えてマップを作成するため、関数化しておく
plot_model <- function(maxV, minV, title, legend, fillColor, xlim, ylim, xbk, ybk){
  ###  描画する濃度範囲の設定
  WRF$Col <- ifelse(WRF$conc > maxV, maxV, WRF$Col)
  WRF$Col <- ifelse(WRF$conc < minV, minV, WRF$Col)
  
  ###  世界地図のデータ
  outlines <- as.data.frame(map("worldHires", plot = FALSE,
                                xlim = c(50, 180),
                                ylim = c(-10, 90)
  )[c("x","y")])
  worldmap <- geom_path(aes(x, y), inherit.aes = FALSE,
                        data = outlines, show.legend = FALSE, col="black")
  
  ###  地図投影法をlccに設定する
  ###  WRFの計算条件に合わせてlat0, lat1, orientationは変更する
  projection <- coord_map(projection="lambert", lat0=30, lat1=60,
                            orientation=c(87.5,30,82.8))
  
  ###  カラースケールを最大値と最小値の間で10分割したラベルを表示
  breaks <- seq(minV, maxV, by=(maxV-minV)/10)

  ###  ■のポイントでプロットする。拡大に応じてsizeを変更する
  DATA <- geom_point(data=WRF, aes(x=LON, y=LAT, z=conc, col=Col), size=1.5, shape=15)

  ggplot() + DATA + worldmap + projection +
      ###  geom_contour()で等圧線の描画
      ###  4hPaごとに描画設定
      geom_contour(data=PRES, aes(x=xval, y=yval, z=presVal), binwidth=4, colour="black", size=0.7) +    
      
      ###  geom_segment()で風向・風速に対応した矢印を描画
      ###  見やすいように適当な数字でU, Vを割る
      geom_segment(data=wind, aes(x=LON-U/5, y=LAT-V/5, xend=LON, yend=LAT), arrow=arrow(length=unit(0.1, "cm")), size=0.25) +
      
      scale_colour_gradientn(colours=fillColor, breaks=breaks, guide="colorbar") +　
      scale_x_continuous(breaks=xbk, expand = c(-0.23, 0), limits=xlim) +
      scale_y_continuous(breaks=ybk, expand = c(0, -2.0), limits=ylim) +
      labs(x="Longitude(deg)", y="Latitude(deg)", title=title, color=legend) +
      theme(legend.key.width = unit(2, "cm")) +
      theme(legend.key.height = unit(0.3, "cm")) +
      theme(legend.position="bottom") +
      theme(plot.title=element_text(size=15)) +
      theme(axis.title=element_text(size=15)) +
      theme(axis.text=element_text(size=12, colour="black")) +
      theme(legend.title=element_text(size=15)) +
      theme(legend.text=element_text(size=12)) +
      theme(panel.background=element_rect(colour="black"))

}

####カラーパレットの作成

##  colorbrewerでも綺麗な図ができるが、個人的に濃度変化を見やすい色で設定
colmap <- colorRampPalette(c(NA, "#4682b4", "#1e90ff", "#20b2aa", "#00ff00", "#ffff00", "#ffa500", "#FF4500", "#ff0000"))

####地図の描画

plot_model(maxV=30,                           ## 濃度プロットの最高濃度
           minV=-30,                          ## 濃度プロットの最低濃度
           title="TEMP[degreeC] color image and PRESSURE[hPa] contour\ndata time; 2014-01-01 09:00 JST\nContours; lines for every 4 hPa",  ## グラフのタイトル
           legend="degreeC",                  ## 凡例のタイトル
           fillColor=colmap(9),               ## 使用するカラーパレット(9色)
           xlim=c(89.5,153.4),                ## 描画範囲（経度）,
           ylim=c(14.5,51.8),                 ## 描画範囲（緯度）
           xbk=seq(60,170,by=10),             ## x軸ラベル間隔
           ybk=seq(0,60,by=10))               ## y軸ラベル間隔

　等圧線に数値を印字したい場合は、directlabelsパッケージを使うと簡単なようだが、だいたいの気圧配置がわかればいいので、印字していない。
　プロットすると、以下のような図ができる。