使用R语言对netCDF栅格数据集进行时间和地理子集的处理

Question

对于以下2016年全球每日海表面温度的netcdf文件，我试图 (i) 在时间上进行子集，(ii) 在地理上进行子集，(iii) 然后对每个像素取长期平均值并创建一个基本图。

文件链接: 这里

library(raster)
library(ncdf4)

设置工作目录后打开netcdf文件。

nc_data <- nc_open('sst.day.mean.2016.v2.nc')

更改时间变量，使其易于解释。

time <- ncdf4::ncvar_get(nc_data, varid="time")
head(time)

将日期更改为我可以解释的格式

time_d <- as.Date(time, format="%j", origin=as.Date("1800-01-01"))

现在我想只选取9月1日到10月15日的数据子集，但不知道如何实现...

在进行时间子集操作之后，创建栅格砖块（或堆栈）和地理子集。

b <- brick('sst.day.mean.2016.v2.nc') # I would change this name to my file with time subest

地理子集

b <- crop(b, extent(144, 146, 14, 16))

最后，我希望对所有数据中的每个像素进行平均，并将其分配给单个光栅，然后制作一个简单的图表...感谢任何帮助和指导。

- Peter Houk

Answer 1

在执行b <- brick('sst.day.mean.2016.v2.nc')命令后，我们可以输入b查看栅格数据的信息。

b
# class       : RasterBrick 
# dimensions  : 720, 1440, 1036800, 366  (nrow, ncol, ncell, nlayers)
# resolution  : 0.25, 0.25  (x, y)
# extent      : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : C:\Users\basaw\Downloads\sst.day.mean.2016.v2.nc 
# names       : X2016.01.01, X2016.01.02, X2016.01.03, X2016.01.04, X2016.01.05, X2016.01.06, X2016.01.07, X2016.01.08, X2016.01.09, X2016.01.10, X2016.01.11, X2016.01.12, X2016.01.13, X2016.01.14, X2016.01.15, ... 
# Date        : 2016-01-01, 2016-12-31 (min, max)
# varname     : sst 

请注意Date插槽包含了从2016-01-01到2016-12-31的信息，这意味着Z值已经包含了日期信息，我们可以使用它来对光栅砖块进行子集操作。

我们可以使用getZ函数来访问存储在Z值中的数值。输入getZ(b)，我们可以看到一系列日期。

head(getZ(b))
# [1] "2016-01-01" "2016-01-02" "2016-01-03" "2016-01-04" "2016-01-05" "2016-01-06"

class(getZ(b))
# [1] "Date"

我们可以使用以下代码来对栅格砖块进行子集划分。

b2 <- b[[which(getZ(b) >= as.Date("2016-09-01") & getZ(b) <= as.Date("2016-10-15"))]]

我们可以根据您提供的代码裁剪图像。

b3 <- crop(b2, extent(144, 146, 14, 16))

要计算平均值，只需使用mean函数。

b4 <- mean(b3, na.rm = TRUE)

最后，我们可以绘制平均值。

plot(b4)

- www

Answer 2

子集和平均任务在CDO中很容易完成：

cdo timmean -sellonlatbox,lon1,lon2,lat1,lat2 -seldate,date1,date2 in.nc out.nc

其中lon1、lon2等定义了要切割的经纬度区域，date1、date2是日期范围。

您可以使用气候运算符包直接从R中调用此命令，就像这个问题一样。

例如，不使用管道，在R中将被拆分成3行：

cdo("seldate,date1,date2",in.fname,out1.fname,debug=TRUE)
cdo("sellonlatbox,lon1,lon2,lat1,lat", out1.fname,out2.fname,debug=TRUE) 
cdo("timmean",out2.fname,out.fname,debug=TRUE) 

- Adrian Tompkins

---

有些文件，例如CMIP5气候模型（https://esgf-node.llnl.gov/search/cmip5/）中有一些时间变量，如average_T1和average_T2。当我运行您的示例时，我收到了以下消息。是否有一种方法可以强制使用时间变量？是否有一种方法可以在此代码中选择变量？警告（find_time_vars）：找到多个时间变量，跳过变量average_T1！警告（find_time_vars）：找到多个时间变量，跳过变量average_T2！60 cdo ntime：在60个时间步骤上处理了2个变量。 - fvfaleiro

@fvfaleiro 抱歉现在才看到。当存在多个时间变量时，CDO可能会出现问题，这也是S2S数据库的问题，我通常尝试构建我的检索以避免这种情况，或者在这些情况下使用ECWMF的gribex或eccodes进行操作。在这种情况下，CDO是否剥离了一个时间变量，还是尽管有警告，操作仍然有效？ - ClimateUnboxed

Adrian Tompkins执行了该操作，尽管收到了警告。输出结果是有意义的，但我不确定它是否正确。 - fvfaleiro

输出的维度是否相同，还是失去了一个时间轴？使用ncdump -h命令看到了什么？ - ClimateUnboxed

尺寸相同，但是预期的时间从time = UNLIMITED; //（目前为240）变为time = UNLIMITED; //（目前为1）。我最关心的是每个单元格中的值，我不确定计算是否正确执行了这些警告。这是ncdum -h输出的第一行：netcdf pr_mean_204101-206012 {dimensions：time = UNLIMITED; //（目前为1）bnds = 2; lon = 144; lat = 90;variables：double time（time）;time：standard_name =“time”;time：long_name =“time”;time：bounds =“time_bnds”;time：units =“days since 2006-01-01 00:00:00” - fvfaleiro

我认为从错误信息和您在此处的回复来看，事情进展得很顺利。它已经对时间变量进行了平均处理，并忽略了其他两个变量。新文件中唯一可能“不正确”的是，average_T1和average_T2这两个变量可能不再正确。我不熟悉CMIP中的这些变量，但我认为它们是某个时期内的平均时间步长。CDO没有触及它们，只是将它们复制到了新文件中，因此它们在新的平均文件中无效。如果您在分析此文件时仅使用“时间”，并忽略这两个变量，则可以正常使用。 - ClimateUnboxed