NetCDFについて
大気・海洋系の研究で用いられているNetCDFファイルは変数、次元数、属性情報を持ったファイル形式で、4次元データの保存に適している。CやFortranで扱われることが多いNetCDFを、Rで扱うことができるようにメモを残します。
Rでは、ncdf4やRNetCDFといった、NetCDFファイルを扱うことができるパッケージが存在します。今回はNetCDF4に対応したncdf4パッケージで説明します。
NetCDFについての詳しい説明は省略します。Fortran版NetCDFユーザーマニュアルに日本語版のテキストがありますので、参照してください。C版のマニュアルも存在します。
RでNetCDFを扱うために
まず、ncdf4パッケージをインストールします。
install.packages("ncdf4")
library(ncdf4)
ncdf4パッケージでよく使用する関数を説明します。ncdf4パッケージのマニュアルはCRANサイトを参照してください。
NetCDFファイルをRで読み込むときはnc_open()関数を使います。
試しにWRFの計算結果を読み込みます。
nc <- nc_open("wrfout_d01_2014-02-01_00_00_00")
このままprint(nc)とすると、ncdump -hとした時と同じ結果が返ってきます。
File ./wrfout_d01_2014-02-01_00_00_00 (NC_FORMAT_64BIT):
146 variables (excluding dimension variables):
char Times[DateStrLen,Time]
float LU_INDEX[west_east,south_north,Time]
FieldType: 104
MemoryOrder: XY
description: LAND USE CATEGORY
units:
stagger:
coordinates: XLONG XLAT
.
.
.
.
.
float LANDMASK[west_east,south_north,Time]
FieldType: 104
MemoryOrder: XY
description: LAND MASK (1 FOR LAND, 0 FOR WATER)
units:
stagger:
coordinates: XLONG XLAT
float SST[west_east,south_north,Time]
FieldType: 104
MemoryOrder: XY
description: SEA SURFACE TEMPERATURE
units: K
stagger:
coordinates: XLONG XLAT
10 dimensions:
Time Size:24 *** is unlimited ***
DateStrLen Size:19
west_east Size:94
south_north Size:74
bottom_top Size:37
bottom_top_stag Size:38
soil_layers_stag Size:4
west_east_stag Size:95
south_north_stag Size:75
land_cat_stag Size:24
98 global attributes:
TITLE: OUTPUT FROM WRF V3.4.1 MODEL
DX: 81000
DY: 81000
.
.
.
多すぎるので省略します。
NetCDFファイルからデータを抽出する際には、変数(variable)と次元数(dimension)の情報をもとに読み解いていきます。
変数は、データの名前であり、WRFの出力では気温(T)や地表面気圧(PSFC)、緯度（XLONG)、経度(XLAT)など多くの変数があります。
次元数は緯度、経度、高さ、時間などの情報に該当し、その内容はprint(nc)したときにdimentions:として出力されます。
例えば、print()で出力された、気温(T)の情報を見ていくと
float T[west_east,south_north,bottom_top,Time]
FieldType: 104
MemoryOrder: XYZ
description: perturbation potential temperature (theta-t0)
units: K
stagger:
coordinates: XLONG XLAT
このように出力されています。変数名（T)のあとに次元数(west_east,south_north,bottom_top,Time)がありますが、これは、T変数が経度(west_east)、緯度(south_north)、高さ(bottom_top)、時間(Time)の次元数を持っていることがわかります。また、単位(units)はケルビン(K)です。
この変数Tのデータを抽出してみます。
T <- ncvar_get(nc, "T")
抽出したTの次元数を確認すると...
> dim(T)
[1] 94 74 37 24
となって、先ほどの次元数(west_east,south_north,bottom_top,Time)の順番で出力されます。なので、東西4・南北12の位置の、最下層・10時間目のデータを抽出したい場合は次のように記述する。
> T[4, 12, 1, 10]
[1] -0.7268804
## ncvar_get(nc, "T")[4, 12, 1, 10]としても同じ
ある次元のデータをすべて書き出したい場合は、次元数を空白にする。試しに東西4・南北12の位置・10時間目のデータを最下層から最上層まですべて出力する。
> T[4, 12, , 10]
[1] -0.72688037 -0.74692106 -0.74950713 -0.74189156 -0.38697922
[6] 0.03983501 0.49886507 1.57127380 3.51204228 6.52097797
[11] 9.38071728 12.41840744 15.76916504 18.81344604 21.68756104
[16] 24.52053833 27.69459343 31.33041191 34.73131561 37.54943085
[21] 39.70961761 41.24932861 42.43940735 43.66003036 45.37441635
[26] 47.49922943 49.54561234 51.56838989 53.93046951 56.02065659
[31] 59.96704483 66.28604889 77.90301514 95.81775665 118.12465668
[36] 140.97705078 164.84016418
次に、同じ条件で、最下層から最上層までの平均値を抜き出したい場合、mean(T[4, 12, , 10]とすると求まる。しかし、これだと2つ以上の次元では使えない。例えば、mean(T[4, 12, , ])として、時間ごとに最下層から最上層までの平均値を求めることはできない。
このようなときはapply()関数を使う。
> x <- T[4, 12, , ]
> dim(x)
[1] 37 24 ## 高さ、時間
> mean(x)
[1] 38.79975
> apply(x, 2, mean)
[1] 38.72858 38.71227 38.72974 38.75706 38.78006 38.79540 38.81224 38.82204
[9] 38.84623 38.85404 38.86388 38.85925 38.84810 38.85500 38.85953 38.85500
[17] 38.84894 38.85622 38.86247 38.83385 38.78480 38.72373 38.67199 38.63350
apply(data, margin, function)関数では、dataで指定した配列に、functionで指定した関数をmarginで指定した次元に適応した結果を返す。
marginの数字は次元の順番であり、(west_east,south_north,bottom_top,Time)の次元を持つデータでは、west_east,south_north,bottom_top,Timeの順番で1〜4の数字が割り振られる。
上記の例では、2次元の配列xに平均値を求める関数meanを時間ごとに適応させた結果が返ってきた。marginを1にすると、高さごとの、時間平均値が返る。
> apply(x, 1, mean)
[1] -1.0751176 -1.1103524 -1.1106088 -1.0420452 -0.8099559 -0.3352173
[7] 0.2404019 1.7825204 4.3521008 7.0892539 9.5852735 12.4250741
[13] 15.6243184 18.6672581 21.5216343 24.3097252 27.5070387 31.0727283
[19] 34.4129815 37.2718484 39.4816875 41.0503216 42.3647316 43.7326900
[25] 45.3379526 47.2209280 49.2762682 51.5797598 54.2426057 56.3446231
[31] 60.2103895 66.4429398 78.1916310 96.3676236 118.6383750 140.8374748
[37] 163.8917669
もし、緯度・経度ユークリッド平面上での最下層の期間平均値を求める場合は、次のように記述する。
> head(apply(T[,,1,], c(1,2), mean))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01576137 0.04714638 0.05050302 0.02025020 -0.02503504 -0.03634928
[2,] -0.03251573 -0.07180138 -0.04693846 -0.04799919 -0.09376347 -0.19082553
[3,] -0.05194719 -0.06634700 -0.09137851 -0.11628789 -0.13117390 -0.16278988
[4,] -0.08131074 -0.03762976 -0.09773546 -0.34580672 -0.54329378 -0.48351659
[5,] -0.09038579 0.00158573 -0.08224722 -0.32997115 -0.91169928 -0.93532613
[6,] -0.13368743 -0.02925559 -0.05932397 -0.20304477 -0.76751008 -0.79081073
.
.
.
.
.
[,68] [,69] [,70] [,71] [,72] [,73] [,74]
[1,] -44.89541 -47.00260 -50.54811 -53.16967 -52.98195 -53.79609 -46.89111
[2,] -43.25560 -46.82128 -49.70825 -51.41909 -52.70751 -53.46209 -55.92016
[3,] -39.67345 -43.78516 -47.86383 -50.77978 -52.87378 -53.38339 -54.67168
[4,] -36.44224 -40.21573 -46.46169 -51.48956 -53.84814 -53.79808 -53.39638
[5,] -36.27189 -39.76556 -46.03089 -51.26846 -53.19774 -53.35141 -52.82303
[6,] -40.24731 -42.20425 -46.12080 -50.63491 -52.61677 -52.13189 -50.25226
これだけわかれば、NetCDFファイルから任意のデータを抽出し、平均値などに加工してグラフにすることが容易になります。
NetCDFファイルの作成
WRFやCMAQの出力はNetCDF形式であり、計算結果が変数に分割されて格納されている。WRFの結果でも湿度や海面更正気圧、風向、風速などはよく使用するが、そのままデータが格納されているわけではなく、いくつかの変数に別れて格納されているため、計算式に当てはめて変換する必要がある。
1地点の地上付近時系列データが欲しい場合は、先ほどのようにncvar_get()やapply()関数を使って時系列データを抽出し、csvファイルなどに書きだしておけば良い。しかし、緯度・経度ユークリッド平面空間でのデータを時間ごとにほしい場合は、csvファイルを時間数分作らなければならない(csvファイルは2次元までしかデータを格納できない)。CMAQは粒子の結果がいくつかの変数に別れて格納されているため、データを必要とする度に変換するのはなかなかめんどくさい。
そこで、使用するデータを、先にすべて変換し、まとめてNetCDFファイルに保管する。こうするとNetCDFファイルからデータを抜き出すだけで欲しいデータが手に入り、更に格納変数が減るためデータサイズが小さくなり、扱いやすくなるメリットがある。
ncdf4パッケージには一通りNetCDFファイルの編集に関する関数が準備されている。
NetCDFファイルを作る手順
簡単には以下の手順で作成する。
1. 次元の定義
2. 変数の定義
3. NetCDFファイルを作成
4. 変数・属性に値を入力
WRFの結果を編集し、必要なデータだけNetCDF形式にまとめる
まずは例のごとく、WRFのデータを読み込む。
library(ncdf4)
WRF <- nc_open("wrfout_d01_2014-02-01_00_00_00")
print(WRF)とすると、次元情報(dimention)に次のように出力されている。
10 dimensions:
Time Size:24 *** is unlimited ***
DateStrLen Size:19
west_east Size:94
south_north Size:74
bottom_top Size:37
bottom_top_stag Size:38
soil_layers_stag Size:4
west_east_stag Size:95
south_north_stag Size:75
land_cat_stag Size:24
この次元数を使って新しい次元数を指定する。
nt <- 24 ## Time Size
nx <- 94 ## west_east Size
ny <- 74 ## south_north Size
nz <- 37 ## bottom_top Size
date.strength <- 19 ## DateStrLen Size
データを確認していく。まずは時間情報を抜き出してみるとUTC(協定世界時)なので、JST(日本標準時)に変換する。
> ncvar_get(WRF, "Times")
[1] "2014-02-01_00:00:00" "2014-02-01_01:00:00" "2014-02-01_02:00:00"
[4] "2014-02-01_03:00:00" "2014-02-01_04:00:00" "2014-02-01_05:00:00"
[7] "2014-02-01_06:00:00" "2014-02-01_07:00:00" "2014-02-01_08:00:00"
[10] "2014-02-01_09:00:00" "2014-02-01_10:00:00" "2014-02-01_11:00:00"
[13] "2014-02-01_12:00:00" "2014-02-01_13:00:00" "2014-02-01_14:00:00"
[16] "2014-02-01_15:00:00" "2014-02-01_16:00:00" "2014-02-01_17:00:00"
[19] "2014-02-01_18:00:00" "2014-02-01_19:00:00" "2014-02-01_20:00:00"
[22] "2014-02-01_21:00:00" "2014-02-01_22:00:00" "2014-02-01_23:00:00"
str.T <- as.POSIXlt(paste(substring(ncvar_get(WRF, "Times")[1], 1, 10), "00:00:00", sep=" "), "%Y-%m-%d %H:%M:%OS") + 9*60*60
JST <- seq(str.T, by="hour", length=nt)
> JST
[1] "2014-02-01 09:00:00" "2014-02-01 10:00:00" "2014-02-01 11:00:00"
[4] "2014-02-01 12:00:00" "2014-02-01 13:00:00" "2014-02-01 14:00:00"
[7] "2014-02-01 15:00:00" "2014-02-01 16:00:00" "2014-02-01 17:00:00"
[10] "2014-02-01 18:00:00" "2014-02-01 19:00:00" "2014-02-01 20:00:00"
[13] "2014-02-01 21:00:00" "2014-02-01 22:00:00" "2014-02-01 23:00:00"
[16] "2014-02-02 00:00:00" "2014-02-02 01:00:00" "2014-02-02 02:00:00"
[19] "2014-02-02 03:00:00" "2014-02-02 04:00:00" "2014-02-02 05:00:00"
[22] "2014-02-02 06:00:00" "2014-02-02 07:00:00" "2014-02-02 08:00:00"
次に変数情報を抜き出していく。次元数が変数によって異なるが、今回は時間情報を1次元、位置情報を2次元、それ以外を3,4次元として作成します。
## 位置データ
lon <- array(ncvar_get(WRF, "XLONG")[,,1], dim=c(nx, ny)) # degree
lat <- array(ncvar_get(WRF, "XLAT")[,,1], dim=c(nx, ny)) # degree
hgt <- array(ncvar_get(WRF, "HGT")[,,1], dim=c(nx, ny, nt)) # terrain elevation
## 気温、湿度、降水量、気圧、風速、境界層高度データ
t2 <- array(ncvar_get(WRF, "T2"), dim=c(nx, ny, nt)) - 273.15 ## 地上2mの気温(℃)
Q2 <- array(ncvar_get(WRF, "Q2"), dim=c(nx, ny, nt)) ## 地上2mの混合比(水蒸気質量kg/乾燥空気質量kg)
PHUM <- 6.1078*10^(
7.5*(t2)/((t2)+237.3)
)#飽和水蒸気圧の計算。Tetensの式だが、Goff-Gratchの式と近似できる。
rh <- (100*psfc*Q2)/(PHUM*(0.622+Q2)) # 地上2mの湿度(%)
rain <- array(ncvar_get(WRF, "RAINC"), dim=c(nx, ny, nt)) ## 累積降水量
psfc <- array(ncvar_get(WRF, "PSFC"), dim=c(nx, ny, nt))/100 ## 地表面気圧(hPa)
pres <- psfc*(1-0.0065*hgt/(t2+273.15+0.0065*hgt))^(-5.257) ## 海面更正気圧
pblh <- array(ncvar_get(WRF, "PBLH"), dim=c(nx, ny, nt)) # 境界層高度
U10 <- array(ncvar_get(WRF, "U10"), dim=c(nx, ny, nt)) ## 地上10m東西方向風速(m/s)
V10 <- array(ncvar_get(WRF, "V10"), dim=c(nx, ny, nt)) ## 地上10m南北方向風速
## 4次元気温、風速データ
temp <- array(ncvar_get(WRF, "T"), dim=c(nx, ny, nz, nt)) - 273.15 ## 全層気温(℃)
U <- array(ncvar_get(WRF, "U"), dim=c(nx, ny, nz, nt)) ## 全層東西方向風速(m/s)
V <- array(ncvar_get(WRF, "V"), dim=c(nx, ny, nz, nt)) ## 全層東南北向風速(m/s)
次元を定義します。
tdim <- ncdim_def("Time", "", 1:nt, unlim=TRUE, create_dimvar=FALSE) ## 時間数
xdim <- ncdim_def("west_east","", 1:nx, create_dimvar=FALSE) ## 東西グリッド数
ydim <- ncdim_def("south_north","", 1:ny, create_dimvar=FALSE) ## 南北グリッド数
zdim <- ncdim_def("bottom_to_top","", 1:nz, create_dimvar=FALSE) ## 高さグリッド数
char.dim <- ncdim_def("date_strength","", 1:date.strength, create_dimvar=FALSE) ## 時間の文字列長さ指定("2014-01-01 00:00:00"で19文字)
変数を定義します。
time.var <- ncvar_def("TIME", "Date (JST)", list(char.dim, tdim), prec="char")
long.var <- ncvar_def("LON", "degreeE", list(xdim, ydim), 1.e30, longname="Longitude")
lat.var <- ncvar_def("LAT", "degreeN", list(xdim, ydim), 1.e30, longname="Latitude")
temp2.var <- ncvar_def("TEMP2", "degree C", list(xdim, ydim, tdim), 1.e30, longname="Temperature at 2m")
rh.var <- ncvar_def("RH", "percent (%)", list(xdim, ydim, tdim), 1.e30, longname="Humidity at 2m")
rain.var <- ncvar_def("RAIN", "mm", list(xdim, ydim, tdim), 1.e30, longname="Total Precipitation")
pres.var <- ncvar_def("SEA_PRES", "hPa", list(xdim, ydim, tdim), 1.e30, longname="Sea Level Correction Pressure")
prsfc.var <- ncvar_def("PRSFC", "hPa", list(xdim, ydim, tdim), 1.e30, longname="Surface pressure")
pblh.var <- ncvar_def("PBLH", "m", list(xdim, ydim, tdim), 1.e30, longname="Planetary boundary layer height")
u10.var <- ncvar_def("U10", "m/s", list(xdim, ydim, tdim), 1.e30, longname="west_east Wind speedat 10m")
v10.var <- ncvar_def("V10", "m/s", list(xdim, ydim, tdim), 1.e30, longname="south_north Wind speed at 10m")
temp.var <- ncvar_def("TEMP", "degree C", list(xdim, ydim, zdim, tdim), 1.e30, longname="Temperature")
u.var <- ncvar_def("U", "m/s", list(xdim, ydim, zdim, tdim), 1.e30, longname="west_east Wind")
v.var <- ncvar_def("V", "m/s", list(xdim, ydim, zdim, tdim), 1.e30, longname="south_north Wind speed")
NetCDFファイルを作成します。このとき、先ほど定義した変数名を指定します。
nc <- nc_create(filename="./data.nc",
list(time.var, long.var, lat.var,
temp2.var, rh.var, rain.var,
pres.var, prsfc.var, pblh.var,u10.var, v10.var,
temp.var, u.var, v.var
)
)
属性を定義ます。定義の方法はncatt_put(file_name, varid, attname, value)で、varidに、属性をつける変数名を指定します。varidを"0"にすると、グローバル属性になります。attnameで属性名、valueで属性値を指定します。
今回は簡単にグローバル属性だけ指定します。
ncatt_put(nc, 0, "CEN_LAT", 33.3389205932617) ## 計算中心位置の緯度
ncatt_put(nc, 0, "CEN_LON", 122.964904785156) ## 計算中心位置の経度
作成したNetCDFファイルは、次元、変数はありますがデータが空の状態なので、さっき定義した変数を入力します。定義の方法は、ncvar_put(file_name, varid, vals, start, count)で、varidで値を入れる変数名、valsで入力する値、startで値を入れる最初の位置指定、countで次元を指定する。startはすべて1で問題なく、countは変数の定義の時に指定した次元をそのまま指定すればよい。
ncvar_put(nc, long.var, lon, start=c(1,1), count=c(nx,ny))
ncvar_put(nc, lat.var, lat, start=c(1,1), count=c(nx,ny))
ncvar_put(nc, temp2.var, t2, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, rh.var, rh, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, rain.var, rain, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, pres.var, pres, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, prsfc.var, psfc, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, pblh.var, pblh, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, u10.var, U10, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, v10.var, V10, start=c(1,1,1), count=c(nx,ny,nt))
ncvar_put(nc, temp.var, temp, start=c(1,1,1,1), count=c(nx,ny,nz,nt))
ncvar_put(nc, u.var, U, start=c(1,1,1,1), count=c(nx,ny,nz,nt))
ncvar_put(nc, v.var, V, start=c(1,1,1,1), count=c(nx,ny,nz,nt))
ncvar_put(nc, time.var, JST)
nc_close(nc) ## 最後に作詞したファイルを閉じる。
何故か時間変数の入力は、一番最後でないとうまくいきませんでした。
作成されたファイルを、端末からncdumpしてみます。
$ ncdump -h ./data.nc
netcdf data {
dimensions:
date_strength = 19 ;
Time = UNLIMITED ; // (24 currently)
west_east = 94 ;
south_north = 74 ;
bottom_to_top = 37 ;
variables:
char TIME(Time, date_strength) ;
TIME:units = "Date (JST)" ;
float LON(west_east, south_north) ;
LON:units = "degreeE" ;
LON:_FillValue = 1.e+30f ;
LON:long_name = "Longitude" ;
float LAT(west_east, south_north) ;
LAT:units = "degreeN" ;
LAT:_FillValue = 1.e+30f ;
LAT:long_name = "Latitude" ;
float TEMP2(Time, west_east, south_north) ;
TEMP2:units = "degree C" ;
TEMP2:_FillValue = 1.e+30f ;
TEMP2:long_name = "Temperature at 2m" ;
float RH(Time, west_east, south_north) ;
RH:units = "percent (%)" ;
RH:_FillValue = 1.e+30f ;
RH:long_name = "Humidity at 2m" ;
float RAIN(Time, west_east, south_north) ;
RAIN:units = "mm" ;
RAIN:_FillValue = 1.e+30f ;
RAIN:long_name = "Total Precipitation" ;
float SEA_PRES(Time, west_east, south_north) ;
SEA_PRES:units = "hPa" ;
SEA_PRES:_FillValue = 1.e+30f ;
SEA_PRES:long_name = "Sea Level Correction Pressure" ;
float PRSFC(Time, west_east, south_north) ;
PRSFC:units = "hPa" ;
PRSFC:_FillValue = 1.e+30f ;
PRSFC:long_name = "Surface pressure" ;
float PBLH(Time, west_east, south_north) ;
PBLH:units = "m" ;
PBLH:_FillValue = 1.e+30f ;
PBLH:long_name = "Planetary boundary layer height" ;
float U10(Time, west_east, south_north) ;
U10:units = "m/s" ;
U10:_FillValue = 1.e+30f ;
U10:long_name = "west_east Wind speedat 10m" ;
float V10(Time, west_east, south_north) ;
V10:units = "m/s" ;
V10:_FillValue = 1.e+30f ;
V10:long_name = "south_north Wind speed at 10m" ;
float TEMP(Time, bottom_to_top, west_east, south_north) ;
TEMP:units = "degree C" ;
TEMP:_FillValue = 1.e+30f ;
TEMP:long_name = "Temperature" ;
float U(Time, bottom_to_top, west_east, south_north) ;
U:units = "m/s" ;
U:_FillValue = 1.e+30f ;
U:long_name = "west_east Wind" ;
float V(Time, bottom_to_top, west_east, south_north) ;
V:units = "m/s" ;
V:_FillValue = 1.e+30f ;
V:long_name = "south_north Wind speed" ;
// global attributes:
:CEN_LAT = 33.3389205932617 ;
:CEN_LON = 122.964904785156 ;
}
もちろん、ncviewで開くこともできます。
海面更正気圧や湿度はすべて計算してあり、また、元は500MB程度のファイルサイズが80MB程度になっていますので、次回から簡単に扱うことができます。
ただ、次元、属性をきちんと付けないと、WRF専用のviewerでは開けないみたいです。今回はデータの保存を目的としたNetCDFファイルの作成ですのでここまでにします。