如何让ERA5的netCDF文件与WATCH数据结构完全匹配？

我有两个NetCDF文件，目标是让ERA5的数据结构完全匹配WATCH文件，经过空间和时间聚合预处理后，两者变量仍不一致（WATCH含Tair+timestp，处理后ERA5含Tair+crs），以下是具体数据详情及解决步骤：

WATCH数据结构 #

> watch_Tair:

     2 variables (excluding dimension variables):
        int timestp[tstep]   
            title: time steps
            units: time steps (days) since 2018-01-01 00:00:00
            long_name: time steps (days) since start of month
        float Tair[lon,lat,tstep]   
            title: Tair average
            units: K
            long_name: Average Near surface air temperature 
 at 2 m at time stamp
            actual_max: 314.173614501953
            actual_min: 212.381393432617
            _FillValue: 1.00000002004088e+20

     4 dimensions:
        lon  Size:720 
            title: Longitude
            units: grid box centre degrees_east
            actual_max: 179.75
            actual_min: -179.75
        lat  Size:360 
            title: Latitude
            units: grid box centre degrees_north
            actual_max: 89.75
            actual_min: -89.75
        tstep  Size:31   *** is unlimited *** (no dimvar)
        day  Size:31 
            title: day
            units: integer
            long_name: day of the month

原始ERA5数据结构 #

> era5_Tair:

     1 variables (excluding dimension variables):
        short t2m[longitude,latitude,time]   
            scale_factor: 0.00169511169020999
            add_offset: 265.715689919741
            _FillValue: -32767
            missing_value: -32767
            units: K
            long_name: 2 metre temperature

     3 dimensions:
        longitude  Size:3600 
            units: degrees_east
            long_name: longitude
        latitude  Size:1801 
            units: degrees_north
            long_name: latitude
        time  Size:744 
            units: hours since 1900-01-01 00:00:00.0
            long_name: time
            calendar: gregorian

已完成的预处理步骤 #

# 日尺度聚合
era5_Tair_daily <- tapp(era5_Tair, "days", mean) 

# 空间重采样匹配WATCH分辨率
era5_temp_resampled <- terra::resample(era5_Tair_daily, watch_Tair, method = "bilinear")

保存后的处理后ERA5数据结构：

> era5_temp_resample:

     2 variables (excluding dimension variables):
        float Tair[longitude,latitude,time]   (Contiguous storage)  
            units: K
            _FillValue: -1.17549402418441e+38
            long_name: Average Near surface air temperature at 2 m at time stamp
            grid_mapping: crs
        int crs[]   (Contiguous storage)  
            crs_wkt: GEOGCRS["unknown",    DATUM["World Geodetic System 1984",        ELLIPSOID["WGS 84",6378137,298.257223563,            LENGTHUNIT["metre",1]],        ID["EPSG",6326]],    PRIMEM["Greenwich",0,        ANGLEUNIT["degree",0.0174532925199433],        ID["EPSG",8901]],    CS[ellipsoidal,2],        AXIS["longitude",east,            ORDER[1],            ANGLEUNIT["degree",0.0174532925199433,                ID["EPSG",9122]]],        AXIS["latitude",north,            ORDER[2],            ANGLEUNIT["degree",0.0174532925199433,                ID["EPSG",9122]]]]
            spatial_ref: GEOGCRS["unknown",    DATUM["World Geodetic System 1984",        ELLIPSOID["WGS 84",6378137,298.257223563,            LENGTHUNIT["metre",1]],        ID["EPSG",6326]],    PRIMEM["Greenwich",0,        ANGLEUNIT["degree",0.0174532925199433],        ID["EPSG",8901]],    CS[ellipsoidal,2],        AXIS["longitude",east,            ORDER[1],            ANGLEUNIT["degree",0.0174532925199433,                ID["EPSG",9122]]],        AXIS["latitude",north,            ORDER[2],            ANGLEUNIT["degree",0.0174532925199433,                ID["EPSG",9122]]]]
            proj4: +proj=longlat +datum=WGS84 +no_defs
            geotransform: -80 0.5 0 15 0 -0.5

     3 dimensions:
        longitude  Size:100 
            units: degrees_east
            long_name: longitude
        latitude  Size:70 
            units: degrees_north
            long_name: latitude
        time  Size:31 
            units: days since 1970-1-1
            long_name: time
            calendar: standard

实现结构完全匹配的步骤 #

1. 删除自动生成的crs变量 #

用terra包直接删除era5_temp_resampled中的crs变量：

# 删除crs变量
era5_temp_resampled <- del(era5_temp_resampled, "crs")

2. 创建并添加timestp变量 #

根据WATCH中timestp的定义（从2018-01-01起的天数步长），生成对应数据并添加：

# 生成timestp数据：1到31的整数（对应1月31天）
timestp_data <- 1:31
# 创建timestp变量，关联时间维度
timestp_var <- terra::rast(nrows=1, ncols=1, nlyrs=31, vals=timestp_data)
# 设置变量属性
terra::set.names(timestp_var, "timestp")
terra::units(timestp_var) <- "time steps (days) since 2018-01-01 00:00:00"
terra::longname(timestp_var) <- "time steps (days) since start of month"
terra::title(timestp_var) <- "time steps"

# 将timestp变量合并到ERA5数据中
era5_temp_resampled <- c(era5_temp_resampled, timestp_var)

3. 调整维度名称匹配WATCH #

把ERA5的longitude/latitude/time维度名改为lon/lat/tstep：

# 修改空间维度名称
names(terra::ext(era5_temp_resampled)) <- c("lon", "lat")
# 修改时间维度名称
names(dimensions(era5_temp_resampled))[3] <- "tstep"

4. 补充缺失的元数据 #

添加WATCH中的day维度，以及Tair变量的缺失属性：

# 添加day维度（大小31，对应1月的天数）
era5_temp_resampled <- add_dim(era5_temp_resampled, "day", 31)
# 设置day维度属性
terra::dim_attr(era5_temp_resampled, "day", "title") <- "day"
terra::dim_attr(era5_temp_resampled, "day", "units") <- "integer"
terra::dim_attr(era5_temp_resampled, "day", "long_name") <- "day of the month"

# 补充Tair变量的属性
terra::title(era5_temp_resampled$Tair) <- "Tair average"
terra::set.values(era5_temp_resampled$Tair, "actual_max", max(values(era5_temp_resampled$Tair), na.rm=TRUE))
terra::set.values(era5_temp_resampled$Tair, "actual_min", min(values(era5_temp_resampled$Tair), na.rm=TRUE))
# 统一FillValue为WATCH的数值
terra::set.values(era5_temp_resampled$Tair, "_FillValue", 1.00000002004088e+20)

5. 保存最终匹配的NetCDF文件 #

writeCDF(era5_temp_resampled, filename = "./era5_temp_matched.nc", varname = c("Tair", "timestp"), 
         unit=c("K", "time steps (days) since 2018-01-01 00:00:00"), 
         zname="tstep", overwrite = TRUE)

内容的提问来源于stack exchange，提问作者herakles_1950