如何在R中从每小时的NetCDF数据创建月度系列?

jaxagkaj  于 2023-09-27  发布在  Etcd
关注(0)|答案(1)|浏览(179)

我正在处理存储在NetCDF文件中的NEE(网络生态系统交换数据),我将其命名为 NEE_2022.nc。文件在我的wd里。空间覆盖范围是整个欧洲,分辨率约为1公里;时间范围是2022年,按小时时间分辨率。在这里,我提供了一些关于数据的更多信息:

**1 variables** (excluding dimension variables):
        float NEE[lon,lat,time]   
            long_name: net ecosystem exchange
            units: micromoles/m2/s
            _FillValue: 1.00000001504747e+30
            missing_value: 1.00000001504747e+30
**3 dimensions**:
        time  Size:8760   *** is unlimited *** 
            standard_name: time
            long_name: time endpoint
            units: hours since 2022-01-01 00:00:00
            calendar: standard
            axis: T
        lon  Size:400 
            standard_name: longitude
            long_name: longitude
            units: degrees_east
            axis: X
        lat  Size:480 
            standard_name: latitude
            long_name: latitude
            units: degrees_north
            axis: Y

**以下是我的目标:**我想平均每个像素每月;我还想为每个月的NEE光栅创建一个Map。不需要空间子集。
注意,R.version命令返回以下文本:

platform       x86_64-w64-mingw32               
arch           x86_64                           
os             mingw32                          
crt            ucrt                             
system         x86_64, mingw32                  
status                                          
major          4                                
minor          2.2                              
year           2022                             
month          10                               
day            31                               
svn rev        83211                            
language       R                                
version.string R version 4.2.2 (2022-10-31 ucrt)

我打开文件时遇到问题。当我运行以下代码时:

library(ncdf4)
library(raster) 
library(rgdal) 
library(ggplot2) 
library(terra)

NEE_2022 <- nc_open('NEE_2022.nc',readunlim=FALSE)
NEE.array <- ncvar_get(NEE_2022, "NEE")

我有一些记忆问题。我收到了以下错误消息:
错误:无法分配大小为12 gb的向量
然后,我通过在ncvar_get()函数中使用“start”和“count”参数成功地创建了数组。我只选择了第一天(24小时)。我做了以下操作:

NEE.array <- ncvar_get(NEE_2022, "NEE", start = c(1,1,1), count = c(-1,-1,24))

我该怎么办?

krugob8w

krugob8w1#

由于NEE文件太大,内存无法容纳,因此必须将其分段处理。

解决方案1:显而易见的方式

将每小时数据汇总为每日数据是最容易理解的方法,之后您必须将每日数据汇总为每月数据(在您的问题中,您提到“每月平均每个像素”,但我假设您希望汇总为“每月NEE”,而不是平均为“每月平均每小时NEE”)。执行此操作的基本循环如下所示:

library(ncdf4)
NEE_2022 <- nc_open('NEE_2022.nc')

# Number of days per month
days_in_month <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)

# First hour of every month in 2022
h1 <- c(1, 745, 1417, 2161, 2881, 3625, 4345, 5089, 5809, 6553, 7273, 8017)

# Array for the result
NEE_data <- array(0, dim = c(480, 400, 12))

for (m in 1:12) {
    mon <- array(0, dim = c(480, 400))
    for (d in 1:days_in_month[m]) {
        data <- ncvar_get(NEE_2022, "NEE", 
                          start = c(1, 1, h1[m]), 
                          count = c(-1, -1, 24))
        mon <- mon + rowSums(data, dims = 2) * 3600
    }
    NEE_data[, , m] <- mon
}
nc_close(NEE_2022)

你也可以将两个循环合并成一个,代价是从文件中读取更大的数据,这可能会耗尽你的内存:

for (m in 1:12) {
    data <- ncvar_get(NEE_2022, "NEE", 
                      start = c(1, 1, h1[m]), 
                      count = c(-1, -1, 24 * days_in_month[m]))
    NEE_data[, , m] <- rowSums(data, dims = 2) * 3600
}

解决方案2:R方式

R think 中,你应该用向量来构建你的问题,尽可能地避免循环。虽然由于文件的大小而无法避免循环,但仍然可以对处理逻辑进行向量化。在这种情况下,这意味着在time维度上创建一个 factor,然后使用它在一次操作中提取所有月度总和,由于文件大小,以逐像素为基础:

# Create the factor
f <- as.factor(rep(rep(1:12, days_in_month), each = 24))

# Loop through every pixel
for (lat in 1:480) {
    for (lon in 1:400) {
        data <- ncvar_get(NEE_2022, "NEE", 
                          start = c(lat, lon, 1), 
                          count = c(1, 1, -1))
        NEE_data[lat, lon, ] <- tapply(data, f, sum) * 3600
    }
}

然而,由于读取的数量非常大,这是非常缓慢的。一种更高效的方法是通过阅读和处理更大的块来减少循环次数,这是一次性处理整行的纬度:

for (lat in 1:480) {
    data <- ncvar_get(NEE_2022, "NEE", 
                      start = c(lat, 1, 1), 
                      count = c(1, -1, -1))
    NEE_data[lat, , ] <- apply(data, 1, tapply, f, sum) * 3600
}

潜在的陷阱

上面的代码通常应该放在一个函数中,这样它就可以很容易地存储和重用。问题是,上面的代码高度依赖于您正在处理的文件的细节。如果要处理另一个文件,而该数据覆盖地球仪的另一个区域,则尺寸可能会更改。如果你处理2020年的数据,有一个闰日是你必须考虑的。
如果数据文件符合CF元数据约定(很有可能),您可以查看CFtime package,因为它可以自动处理时间序列,并且它将通过阅读netCDF文件中的详细信息来生成第二个解决方案中使用的因子。

相关问题