从netCDF文件中在R中计算月平均值

Question

我有一个包含16年（1998-2014）日降水量（5844层）的netCDF文件(.nc)。三个维度是时间（大小为5844），纬度（大小为19）和经度（大小为20）。 在R中，是否有一种简单的方法来计算每个栅格单元格：

• 每月和每年平均值
• 累积比较（例如，jan-mar与所有jan-mar的平均值相比）

到目前为止，我有：

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

希望我的问题表述清楚了，初步的正确方向指引将不胜感激。 样本数据在这里。

- Jobbo90

---

我发现样本nc文件的维度与您在问题中提到的不同。您能否上传一个更好的样本数据（与您在问题中提到的完全相同的维度）？ - raymkchow

你是对的，我假设只需要几天，尺寸就会保持不变。由于原始文件不太大，我提供了新的样本数据！谢谢@raymkchow - Jobbo90

我建议将您的数据转换为 xts 类型，并使用函数 apply.monthly 和 apply.yearly 进行计算。但也许 @joberlin 的方法更好，因为它使用了 rasterstack。 - raymkchow

Answer 1

虽然问题要求使用 R 进行解决，但如果有人想要一个简单的命令行解决方案来完成此任务，这些统计数据是 CDO 的基础。

月平均值：

cdo monmean in.nc monmean.nc

年度平均值：

cdo yearmean in.nc yearmean.nc

计算一月、二月等所有月份的平均值：

cdo ymonmean in.nc ymonmean.nc

每月异常相对于长期年度周期：

cdo sub monmean.nc ymonmean.nc monanom.nc

如果您需要选择特定月份，只需使用 selmon 或 seldate 进行选择。

您可以从 R 中使用 system 命令调用这些函数。

- Adrian Tompkins

Answer 2

这里介绍一种使用zoo包的方法：

### 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个层。但是，以下方法也适用于具有更多层（月份）的情况。 首先，我们将构建一个函数，用于操作一个值向量（即像素时间序列），将输入转换为使用dates创建zoo对象，并使用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 检索单元格值后，在单元格值上应用函数或使用来自 raster 包的 calc 函数创建光栅输出（这将是一个具有与数据中月份数量相同层数的光栅堆）。

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

- joberlin

---

看起来你想以错误的方式访问堆栈的层。要获取一个堆栈对象的第190层，请执行以下操作： stack_object[[190]] - joberlin

如果数据是每小时观测一次，我们想要得到每日平均值怎么办？我尝试将monthly_mean_stack的第3行更改为out <- as.numeric(aggregate(pixel.ts, as.Date, mean, na.rm=TRUE))，但最终得到的栅格堆叠图层数与每小时观测次数相同。 - alaybourn

Answer 3

我认为最简单的方法是将其转换为栅格砖，然后再转换成数据框。

接着可以使用通用代码 DF$weeklymean <- rowMeans(DF[, ]) 很容易地提取统计信息。

- TJeff