使用 R 中的 netcdf 数据集计算网格单元的加权空间全球年平均值

Question

我目前正在从事一个涉及存储在 NetCDF 文件中的气候模型数据的项目。我目前正在尝试计算 "weighted" 空间年度 "global" 降水平均值。我需要为我拥有的 95 年的全球降水数据中的每一年执行此操作。这个想法是通过使用其纬度的余弦以某种方式将权重应用于每个网格单元（这意味着赤道处的纬度网格单元的权重为 1（即 0 度的余弦为 1），并且两极将具有值为 1（因为 90 的余弦为 1））。然后，我将能够根据每个网格单元的平均值来计算年度加权平均值。

我知道如何从概念上做到这一点，但我不确定从哪里开始在 R 中编写脚本以将权重应用于所有网格单元格，然后对 95 年中的每一年进行平均。我将不胜感激任何对此的帮助，或任何可能有用的资源！！！

最起码我打开了.nc文件，读入了NetCDF变量，如下图：

ncfname<-"MaxPrecCCCMACanESM2rcp45.nc"
Prec<-raster(ncfname)
print(Prec)
Model<-nc_open(ncfname)
get<-ncvar_get(Model,"onedaymax")
longitude<-ncvar_get(Model, "lon")
latitude<-ncvar_get(Model, "lat")
Year<-ncvar_get(Model, "Year")

此外，假设我想为特定位置或区域创建这些新派生的加权平均值的时间序列，我之前使用以下代码显示 95 年来一天最大降水量的趋势可行，但我只需要稍作更改即可使用年度加权平均值？ :

r_brick<-brick(get, xmn=min(latitude), xmx=max(latitude), ymn=min(longitude),                               
ymx=max(longitude), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84   
+no_defs+ towgs84=0,0,0"))
r_brick<-flip(t(r_brick), direction='y')
randompointlon<-13.178
randompointlat<--59.548
Random<-extract(r_brick, 
SpatialPoints(cbind(randompointlon,randompointlat)),method='simple')
df<-data.frame(year=seq(from=1, to=95, by=1), Precipitation=t(Hope))
ggplot(data=df, aes(x=Year, y=Precipitation,   
group=1))+geom_line()+ggtitle("One-day maximum precipitation (mm/day) trend  
for Barbados for CanESM2 RCP4.5")

此外，如果有帮助，以下是 .nc 文件包含的内容：

 3 variables (excluding dimension variables):
    double onedaymax[lon,lat,time]   (Contiguous storage)  
        units: mm/day
    double fivedaymax[lon,lat,time]   (Contiguous storage)  
        units: mm/day
    short Year[time]   (Contiguous storage)  

 3 dimensions:
    time  Size:95
    lat  Size:64
        units: degree North
    lon  Size:128
        units: degree East

再次强调，任何帮助都将非常有价值！期待您的回复！

Answer 1

请一次问一个清楚的问题，并提供示例数据（通过代码）。

我认为您没有以正确的方式阅读 ncdf 数据。我觉得你应该做

library(raster)
ncfname <- "MaxPrecCCCMACanESM2rcp45.nc"
Prec <- brick(ncfname, var="onedaymax")

（不使用nc_open等）

获取全局加权平均值

示例数据

library(raster)
r <- abs(init(raster(), 'y'))
s <- stack(r, r, r)

s 是一个 RasterStack，极点值为 90，赤道值为 0

未加权的全局平均值。首先对层进行平均，然后对单元进行平均（在这种情况下倒序也适用）

sm <- mean(s, na.rm=TRUE)
cellStats(sm, mean, na.rm=TRUE)
[1] 45

现在使用权重（获得较低的数字是高纬度获得较少的权重）

# raster with latitude cell values 
w <- init(s, 'y')
# cosine after transforming to radians
w <- cos(w  * (pi/180))
# multiply weights with values
x <- sm * w
# compute weighted average
cellStats(x, sum) / cellStats(w, sum)
#[1] 32.70567

另一种可能更简单的解决方案是使用每个单元格的面积（与 cos(lat) 成正比）。结果可能更精确一点（因为面积不只考虑小区中心）。

a <- area(s) / 10000
y <- sm * a
cellStats(y, sum) / cellStats(a, sum)
#[1] 32.72697

以后:

对于时间序列，只需使用s。

未加权

cellStats(s, mean) 
#layer.1 layer.2 layer.3 
# 45      45      45

加权

a <- area(s) / 10000
y <- s * a
cellStats(y, sum) / cellStats(a, sum)
# layer.1  layer.2  layer.3 
#32.72697 32.72697 32.72697

Answer 2

不是我想让你远离 R，而是这种计算是直接从命令行进行的 cdo（气候数据运算符）的绝对面包和黄油！

计算空间加权平均值（这说明了纬度的罪恶，也可以处理减少的高斯网格等）：

cdo fldmean input.nc fldmean.nc

计算年均值

cdo yearmean input.nc yearmean.nc

通过将两者结合起来计算年度全球均值的时间序列（即使用 fldmean.nc 作为第二个命令的输入），或者您可以通过管道在一行上执行：

cdo yearmean -fldmean input.nc yearglobal.nc

那是什么？你想为你说的经纬度盒子区域计算它，而不是全球平均值？没问题，先用sellonlatbox切出一块区域

cdo sellonlatbox,lon1,lon2,lat1,lat2 in.nc out.nc

所以管道这个：

cdo  yearmean -fldmean -sellonlatbox,lon1,lon2,lat1,lat2 in.nc yearregion.nc

但是等等！您现在想要一个特定的位置，而不是区域平均值？好吧，你可以选择离

使用 R 中的 netcdf 数据集计算网格单元的加权空间全球年平均值

Calculating weighted spatial global annual averages across grid cells using netcdf dataset in R

global

r

spatial

weighted-average

netcdf