4

16 年 (1998 年 - 2014 年) 相当の日降水量 (5844 層) を含む 1 つの netCDF ファイル (.nc) があります。3 つの次元は、時間 (サイズ 5844)、緯度 (サイズ 19)、および経度 (サイズ 20) です。R で各ラスターセルを計算する簡単な方法はありますか。

  • 月間・年平均
  • 累積比較 (たとえば、すべての 1 月から 3 月の平均と 1 月から 3 月を比較)

これまでのところ、私は持っています:

library(ncdf4)
library(raster)

Rname <- 'F:/extracted_rain.nc'
rainfall <- nc_open(Rname)
readRainfall <- ncvar_get(rainfall, "rain") #"rain" is float name
raster_rainfall <- raster(Rname, varname = "rain") # also tried brick()
asdatadates <- as.Date(rainfall$dim$time$vals/24, origin='1998-01-01') #The time interval is per 24 hours

私の最初の課題は、各ラスター セルの月平均の計算です。究極の目標（累積比較）を念頭に置いたまま、どのように進めるのが最善かわかりません。特定の月の日付のみに簡単にアクセスするにはどうすればよいですか?

raster(readRainfall[,,500])) # doesn't seem like a straightforward approach

うまくいけば、質問が明確になりました。正しい方向への最初のプッシュは大歓迎です。サンプルデータはこちら

4

3 に答える 3

3

zoo-packageを使用したアプローチの 1 つを次に示します。

### first read the data
library(ncdf4)
library(raster)
library(zoo)

### use stack() instead of raster
stack_rainfall <- stack(Rname, varname = "rain")

### i renamed your "asdatadates" object for simplicity
dates <- as.Date(rainfall$dim$time$vals/24, origin='1998-01-01')

サンプル データセットでは、すべて 1998 年 1 月からの 18 レイヤーしかありません。しかし、以下はより多くのレイヤー (月) でも機能するはずです。zooまず、 を使用して入力をオブジェクトに変換し、 を使用datesして平均を計算するために、値の 1 つのベクトル (つまり、ピクセル時系列) を操作する関数を作成しますaggregate。この関数は、 の月数に等しい長さのベクトルを返しますdates。

monthly_mean_stack <- function(x) {
    require(zoo)
    pixel.ts <- zoo(x, dates)
    out <- as.numeric(aggregate(pixel.ts, as.yearmon, mean, na.rm=TRUE))
    out[is.nan(out)] <- NA     
    return(out)
}

次に、出力をベクトル/マトリックス/データ フレームにするか、ラスター形式のままにするかに応じて、セル値を で取得した後にセル値に関数を適用するか、 - 関数から-関数をgetValues使用できます。パッケージを使用してラスター出力を作成します (これは、データ内の月と同じ数のレイヤーを持つラスター スタックになります)calcraster

v <- getValues(stack_rainfall) # every row displays one pixel (-time series)


# this should give you a matrix with ncol = number of months and nrow = number of pixel
means_matrix <- t(apply(v, 1, monthly_mean_stack))

means_stack <- calc(stack_rainfall, monthly_mean_stack)

大規模なラスター データセットを操作している場合、 関数を使用して関数を並列に適用することもできますclusterR。見る ？clusterR

于 2016-12-20T22:55:06.773 に答える