R 中 netCDF 文件的月平均值答案

【问题标题】：Monthly average from netCDF files in RR 中 netCDF 文件的月平均值
【发布时间】：2017-05-05 20:33:47
【问题描述】：

我有一个 netCDF 文件 (.nc)，其中包含 16 年（1998 - 2014）的每日降水量（5844 层）。 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

希望我把我的问题说清楚了，我们将不胜感激第一次朝着正确的方向前进。样本数据 here

【问题讨论】：

我发现示例nc文件的尺寸与您在问题中提到的不同。您能否上传更好的样本数据（与您在问题中提到的维度完全相同）？
你说得对，我认为只需几天时间，尺寸就会保持不变。由于原始文件不大，我提供了新的示例数据！谢谢@raymkchow
我建议将您的数据转换为xts 类，并使用函数apply.monthly 和apply.yearly 进行计算。但是@joberlin 的方法可能更好，因为它使用了 rasterstack。

标签： r netcdf4

【解决方案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 中调用这些函数。

【讨论】：

【解决方案2】：

这是使用zoo-package 的一种方法：

### 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')

在您的示例数据集中，您只有 18 层，全部来自 1998 年 1 月。但是，以下内容也应该适用于更多层（月）。首先，我们将构建一个函数，该函数操作一个值向量（即像素时间序列），以使用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-package 的calc-function 用于创建栅格输出（这将是一个栅格堆栈，其层数与您的数据中的月份一样多）

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

【讨论】：

您似乎也想以错误的方式访问堆栈的层。要获取堆栈对象的第 190 层，请执行以下操作：stack_object[[190]]
如果数据是每小时的观察值，而我们想要的是每日平均值，该怎么办？我尝试将 Monthly_mean_stack 的第 3 行更改为 out <- as.numeric(aggregate(pixel.ts, as.Date, mean, na.rm=TRUE)) ，但最终得到了一个栅格堆栈，其层数与我每小时观察到的一样多。

【解决方案3】：

我认为最容易转换为栅格砖然后转换为 data.frame。

然后可以很容易地使用通用代码 DF$weeklymean

【讨论】：