使用 R 的 netcdf 栅格堆栈或栅格砖的时间和地理子集
Posted
技术标签:
【中文标题】使用 R 的 netcdf 栅格堆栈或栅格砖的时间和地理子集【英文标题】:time and geographical subset of netcdf raster stack or raster brick using R 【发布时间】:2018-06-09 06:12:06 【问题描述】:对于以下包含 2016 年每日全球海面温度的 netcdf 文件,我尝试 (i) 时间上的子集,(ii) 地理上的子集,(iii) 然后对每个像素采用长期平均值并创建一个基本的情节。
文件链接:here
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))
最后,我想取所有数据天数中每个像素的平均值,将其分配给单个栅格,然后绘制一个简单的绘图...
感谢您的帮助和指导。
【问题讨论】:
【参考方案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)
【讨论】:
【参考方案2】:不是在 R 中,只是指出子集和平均任务在 CDO 中很容易从命令行完成:
cdo timmean -sellonlatbox,lon1,lon2,lat1,lat2 -seldate,date1,date2 in.nc out.nc
其中 lon1,lon2 等定义要剪切的 lon-lat 区域,date1,date2 是日期范围。
然后您可以将生成的文件读入 R 进行绘图,或使用 ncview 快速查看。
【讨论】:
一些文件,如 CMIP5 气候模型 (esgf-node.llnl.gov/search/cmip5),有一些时间变量,如 average_T1 和 average_T2。当我运行您的示例时,我收到以下消息。有没有办法强制时间变量?有没有办法在这段代码中选择变量?警告(find_time_vars):发现多个时间变量,跳过变量average_T1!警告(find_time_vars):发现多个时间变量,跳过变量average_T2! 60 cdontime:处理 2 个变量超过 60 个时间步。 @fvfaleiro 抱歉,现在才看到这个。当存在多个时间变量时,CDO 可能会遇到困难,这也是 S2S 数据库的一个问题,我通常会尝试构建我的检索来避免这种情况,或者在这些情况下使用来自 ECWMF 的 gribex 或 eccodes 进行操作。在这种情况下,CDO 是否去除了其中一个时间变量,或者尽管有警告,操作仍然有效? Adrian Tompkins 尽管有警告,该操作仍被执行。输出是有道理的,但我不确定它是否正确。 输出是否具有相同的尺寸,或者您是否丢失了一个时间轴?你对 ncdump -h 有什么看法? 尺寸相同,但正如预期的那样,时间从 time = UNLIMITED 改变; //(目前为 240)到时间 = UNLIMITED ; //(当前为 1 个)。我主要关心的是每个单元格中的值,我不确定这些警告是否正确执行了计算。这是 ncdum -h 输出的第一行:netcdf pr_mean_204101-206012 dimensions: time = UNLIMITED ; // (当前为 1) bnds = 2 ;隆 = 144 ;纬度 = 90;变量:双倍时间(时间);时间:标准名称=“时间”;时间:长名称=“时间”;时间:界限=“time_bnds”; time:units = "自 2006-01-01 00:00:00 以来的天数"以上是关于使用 R 的 netcdf 栅格堆栈或栅格砖的时间和地理子集的主要内容,如果未能解决你的问题,请参考以下文章