9 大気再解析

再解析データとは，比較的新しい予報モデルに対して過去の観測データを同化したものです。数値天気予報においては，予報モデルに対して観測を同化した解析値が作成され，初期値として用いられています。予報モデルは改良が重ねられていくため，解析値の品質はモデルの変遷と共に変化していきます。データ同化手法も大きく進歩し，数十年前とは全く異なる手法が用いられています。モデルやデータ同化手法の改良の影響を取り除くことにより，時間的により品質が安定したデータセットになります。再解析とは，現業予報で一度行った解析を再びやり直すということを意味しています。

9.1 ラスタデータ

terraで格子点値はSpatRasterクラスで扱います。

SpatRasterにデータが格納される順序を確認します。

library(terra)

terra 1.8.5

r <- rast(ncol=12, nrow=6, xmin=0, xmax=360, ymin=-90, ymax=90)
values(r) <- 1:ncell(r)
plot(r)

列方向に連続する行優先で、列は西から東、行は北から南であることが分かります。セルに色が付いているので半格子ずれるが、ひとまず目をつぶることにします。

9.2 気候値

京都大学生存圏研究所生存圏データベースグローバル大気観測データでは、NCEP再解析を提供しています。月平均/その他の統計処理済データ、surfaceとたどりましょう。

NCEP/NCAR再解析データの海面気圧の月別気候値slp.mon.ltm.ncを取得します。

NetCDFは、RNetCDFやncdf4で読むことができます。 terraもrast()でNetCDFを読めることになっているが、上記ファイルは読めませんでした。

library(RNetCDF)

nc <- open.nc("slp.mon.ltm.nc")
slp <- var.get.nc(nc, "slp")
dim(slp)

[1] 144  73  12

変数はvar.get.nc()で取得します。次元は経度、緯度の順で144x73です。 Rの配列はFortranやMATLABと同じ列優先（左の添え字が先に変わる）なので、rast()に渡すときは、転置t()する必要があります。 SpatRasterも北から南向きなので、南北は反転しません。データと経度を合わせて日付変更線を図の中心にするため、Chapter 8 で定義したshift.lon()を使っています。 1月の海面気圧を描きましょう。

#|echo: false
shift.lon <- function(v, dlon=0) {
  e <- crop(v, ext(dlon+0, 180, -90, 90))
  w <- crop(v, ext(dlon-180, dlon, -90, 90))
  w <- shift(w, 360)
  rbind(e, w)
}

slp.ras <- rast(xmin=0, xmax=360, ymin=-90, ymax=90, ncols=144, nrows=73)
values(slp.ras) <- t(slp[,,1])
cshp <- "ne_50m_coastline/ne_50m_coastline.shp"
c50 <- vect(cshp)
c50 <- shift.lon(c50)
plot(slp.ras)
plot(c50, add=TRUE)

カラーパレット

軸のフォントサイズなどのパラメタはpaxにリストとして指定します。同様にカラーバーのパラメタはplgに指定します。

plot(slp.ras, pax=list(cex.axis=2), plg=list(cex=2))
plot(c50, add=TRUE)

slp.ras.c <- crop(slp.ras, ext(0, 360, -30, 90))

newcrs <- "+proj=stere +lon_0=135e +lat_0=90n"
c50p <- project(c50, newcrs)

slp.ras.p <- project(slp.ras.c, newcrs)
g <- graticule(30, 30, crs=newcrs)
plot(slp.ras.p, axes=FALSE, ext=ext(-1e+7, 1e7, -1e7, 1e7))
plot(c50p, add=TRUE)
plot(g, add=TRUE)

カラーパレット

色を変えるには、color.paletteにカラーパレットを生成する函数を与えるか、colに明示的に色を与えます。colはcolor.paletteに優先し、色の数はレベルより一つ少なくします。

RColorBrewerなど、様々なパレットがパッケージとして提供されています。
NCLのカラーテーブルはrcolorsに収められています。

9.3 月平均

NCEP再解析の月平均/その他の統計処理済データ、surfaceから、今度は月平均データを取得します。ここでは層厚（1000 hPaと500 hPaの厚さ）thickness.mon.mean.ncを選びました。

特定の年月に対応する番号は、1948年1月が最初に入っているので、 (y - 1948) * 12 + mで計算できます。

あるいは、日時に関する函数を使って次のように求めることができます。

library(RNetCDF)

nc <- open.nc("thickness.mon.mean.nc")
1time <- var.get.nc(nc, "time")
2tunit <- att.get.nc(nc, "time", "units")
3time.posixct <- utcal.nc(tunit, time, type="c")
4ymd <- as.POSIXct(paste(c(1970, 1, 1), collapse="-"))
5t <- which.min(abs(time.posixct - ymd))

1: 変数timeを取得します。
2: 時刻の単位を取得します。
3: RNetCDFのutcal.nc()でPOSIXct型に変換します。
4: c()で作った年月日のベクトルを-で繋ぎ、POSIXct型に変換します。
5: which.min()で最も小さい番号を求めます。

練習

同じ月について、最近の層厚と数十年前の層厚を比べてみよう。