使用 R 中的特定条件聚合 nc 文件

Question

我再次需要你的帮助。 我有 .nc 文件，元数据： 文件 minty.nc (NC_FORMAT_64BIT):

 1 variables (excluding dimension variables):
    short mn2t[longitude,latitude,time]   
        scale_factor: 0.000940940342005054
        add_offset: 259.294916797895
        _FillValue: -32767
        missing_value: -32767
        units: K
        long_name: Minimum temperature at 2 metres since previous post-processing

 3 dimensions:
    longitude  Size:57
        units: degrees_east
        long_name: longitude
    latitude  Size:49
        units: degrees_north
        long_name: latitude
    time  Size:90240
        units: hours since 1900-01-01 00:00:00.0
        long_name: time
        calendar: gregorian

2 global attributes:
    Conventions: CF-1.6
    

我有一个代码，适用于较小的 .nc 文件：

 library(raster)
 library(rgdal)
 library(ggplot2)
 nc_data = nc_open('file.nc')
 lon = ncvar_get(nc_data, "longitude")
 lat = ncvar_get(nc_data, "latitude", verbose = F)
 t = ncvar_get(nc_data, "time")
 head(lon)
 head(t)
 head(lat)
 mint.array = ncvar_get(nc_data, "mn2t")
 dim(mint.array)
 fillvalue = ncatt_get(nc_data, "mn2t", "_FillValue")
 fillvalue
 mint.array[mint.array == fillvalue$value] <- NA
 r_brick <- brick(mint.array, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
 r_brick = flip(t(r_brick), direction = 'y')

由于文件太大，我收到错误：“无法分配大小为 1.4 Mb 的向量” 我还使用 gc() 来清除未使用的内存。这没有帮助。 我不需要 file.nc 中的所有数据。在这种情况下，我需要以某种方式聚合它。为了进一步计算，我只需要每日最小值。在本例中，对于 df 我使用：

df(ff) <- aggregate(df, list(rep(1:(nrow(df)%(%n+1), each=24, len=nrow(df))), min)

不幸的是，我发现很难将此代码改编为 .nc 文件。也许有人可以帮助我。预先感谢您！

Answer 1

为了避免内存问题，您可以这样做：

library(raster)
r_brick <- brick('file.nc', "mn2t")

它还可以防止错误。例如，在您的代码中，这在两个方面是错误的：

xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon)

因为

x

应该是

lon

并且

y

应该是

lat

并且因为 ncdf 坐标指的是单元格的中心，而

xmn

、

xmx

、

ymn

和

ymx

指的是单元格的中心到边界。

您也可以使用现代的等价物

library(terra)
r <- rast('file.nc')

Answer 2

使用标准数组逐步引导您完成的答案。

数据文件

您的数据文件符合 CF 元数据约定，如全局属性所示。大多数气候观测和预测都是按照该公约发布的（有一些重要的例外）。这意味着您可以假设有关该文件的许多信息，例如是否存在

longitude

、

latitude

和

time

维度（以及可能的其他维度），它们定义了以 3- 形式存储的数据。一个或多个变量的一维数组。数据也按该维度顺序存储在文件中。

以块的形式读取数据

如果数据文件太大而无法全部加载到内存中，您可以分块加载并在读取更多数据之前对其进行处理。您应该始终按照尺寸顺序执行此操作，因为这可以确保快速读取文件。 （粗略的基准测试表明，沿着子午线（经纬度）读取全时间维度的速度比沿着平行线（经纬度）读取速度慢约 7 倍，后者本身比读取时间片（经纬度）慢约 7 倍）。在这种情况下，这意味着您读取

time

维度的切片。

日历

该文件以

gregorian

ly 时间步长使用

hour

日历。如果您确定

time

维度中没有间隙（每天有24个切片并且没有缺失的日子），您可以使用上面使用的简单算术，但使用可能会更好包

CFtime

可以处理所有 CF 日历（共有 9 个）和时间序列中的间隙。

内存限制和聚合到每日最小值的解决方案

此解决方案故意冗长，以容纳任何缺失的日期或观察结果。

library(ncdf4)
library(CFtime)

nc_data <- nc_open("file.nc")

# Read in the details of the `time` dimension
cf <- CFtime(nc_data$dim$time$units, nc_data$dim$time$calendar, nc_data$dim$time$vals)

# Make a factor of all days in the file, then find the starting index and the
# number of observations per day in the time dimension
days <- CFfactor(cf, "day")         # Missing days not a problem 
nobs <- tabulate(days)              # Missing observations not a problem
h1 <- cumsum(nobs) - nobs[1] + 1

# Create the array of daily minimum values, all NA values
tmin <- array(dim = c(nc_data$dim$longitude$length, nc_data$dim$latitude$length, nlevels(days)))

# Loop through the time dimension, day by day
for (d in 1:nlevels(days)) {
  hrs <- ncvar_get(nc_data, "mn2t", start = c(1, 1, h1[d]), count = c(-1, -1, nobs[d]))
  tmin[, , d] <- apply(hrs, 1:2, min)
}

# Set the dimension names
dimnames(tmin) <- list(nc_data$dim$longitude$vals, nc_data$dim$latitude$vals, levels(days))

nc_close(nc_data)

tmin

数组包含至少有一次观测的每一天的每日最低温度。使用相同的代码模式，您可以计算每日温度的其他有趣属性，例如一天中温度最低的时间。也有可能，如果观察数量低于一定数量，则不会计算每日最小值，因此您应该在结果中添加

NA

。

您可以使用

terra::SpatRaster

函数将此数组简单地转换为

terra::rast()

，但您必须确定数据文件是否具有“点”或“单元”数据才能设置正确的范围。您显示的元数据没有提供这方面的任何线索。

更多内存

如果您有足够的 RAM 来完成内存中的所有操作，则上述内容可简化为：

nc_data <- nc_open("file.nc")
cf <- CFtime(nc_data$dim$time$units, nc_data$dim$time$calendar, nc_data$dim$time$vals)
days <- CFfactor(cf, "day")
tmin <- aperm(apply(ncvar_get(nc_data, "mn2t"), 1:2, tapply, days, min), c(2, 3, 1))
dimnames(tmin) <- list(nc_data$dim$longitude$vals, nc_data$dim$latitude$vals, levels(days))
nc_close(nc_data)

也许升级你的硬件？