将 netCDF 文件裁剪为 shapefile 并在 R 中克隆元数据变量

Question

我有 NetCDF 文件（例如 https://data.ceda.ac.uk/neodc/esacci/lakes/data/lake_products/L3S/v1.0/2019 全局域），我想根据 shapefile 边界提取数据（在本例中为 Lake - https://www.sciencebase.gov/catalog/item/530f8a0ee4b0e7e46bd300dd），然后将剪裁的数据保存为 NetCDF 文件，但是保留剪辑文件中的所有原始元数据和变量名称。这是我到目前为止所做的

library(rgdal)
library(sf)
library(ncdf4)
library(terra)

#Read in the shapefile of Lake 
Lake_shape <- readOGR("C:/Users/CEDA/hydro_p_LakeA/hydro_p_A.shp")
# Reading the netcdf file using Terra Package function rast
test <- rast("ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190705-fv1.0.nc")
# List of some of variables names for orginal dataset 
      head(names(test))
[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"           "lake_surface_water_extent"                 
[4] "lake_surface_water_extent_uncertainty"      "lake_surface_water_temperature"             "lswt_uncertainty"   
                                 
#Clipping data to smaller Lake domain using the crop function in Terra Package
test3 <- crop(test, Lake_shape)
#Listing the some variables names for clipped data
head(names(test3))
[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"           "lake_surface_water_extent"                 
[4] "lake_surface_water_extent_uncertainty"      "lake_surface_water_temperature"             "lswt_uncertainty" 


# Writing the crop dataset as netcdf or Raster Layer using the WriteCDF function 

filepath<-"Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0"
fname <- paste0( "C:/Users/CEDA/",filepath,".nc")
rnc <- writeCDF(test3, filename =fname, overwrite=T)”

我的主要问题是，当我读入剪辑的 netCDF 文件时，我似乎无法保留原始 NetCDF 的数据变量的名称。当我使用 writeCDF 函数将剪辑的数据集保存为新的 netCDF 时，它们都会自动重命名。

#Reading in the new clipped file
 LakeA<-rast("Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0.nc")
> head(names(LakeA))
[1] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_1" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_2"
[3] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_3" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_4"
[5] "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_5" "Lake_A_ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20020501-fv1.0_6"

那么当剪辑到 R 中较小的域/shapefile 时，是否可以从原始 NetCDF 数据集中克隆/复制所有元数据变量，然后另存为 NetCDF？任何有关如何在 R 中执行此操作的指导将不胜感激。 （NetCDF 和 R 对我来说都是新的，所以我不确定我缺少什么或有深入的知识来排序）。

Answer 1

您有一个包含许多 (52) 个变量（子数据集）的 NetCDF 文件。当您使用rast 打开文件时，这些将成为“层”。或者，您可以使用 sds 打开文件以保留子数据集结构，但这对您没有帮助（您需要跳过前两个，见下文）。

library(terra)
f <- "ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc"
r <- rast(f)
r
#class       : SpatRaster 
#dimensions  : 21600, 43200, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:water_surface_height_above_reference_datum  
              ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:water_surface_height_uncertainty  
              ESACCI-LAKES-L3S-LK_PRODUCTS-MERGED-20190101-fv1.0.nc:lake_surface_water_extent  
              ... and 49 more source(s)
#varnames    : water_surface_height_above_reference_datum (water surface height above geoid) 
              water_surface_height_uncertainty (water surface height uncertainty) 
              lake_surface_water_extent (Lake Water Extent) 
              ...
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#unit        :           m,           m,         km2,         km2,      Kelvin,      Kelvin, ... 
#time        : 2019-01-01 

请注意，有 52 个层和源（子数据集）。有名字

head(names(r))
#[1] "water_surface_height_above_reference_datum" "water_surface_height_uncertainty"          
#[3] "lake_surface_water_extent"                  "lake_surface_water_extent_uncertainty"     
#[5] "lake_surface_water_temperature"             "lswt_uncertainty"                          

还有“长名”（它们通常比变量名长得多，在这种情况下不是）

head(longnames(r))
# [1] "water surface height above geoid" "water surface height uncertainty" "Lake Water Extent"               
# [4] "Water extent uncertainty"         "lake surface skin temperature"    "Total uncertainty"               

你也可以用sds打开文件，但需要跳过“lon_bounds”和“lat_bounds”变量（维度）

s <- sds(f, 3:52)

现在读取一个矢量数据集（在本例中为 shapefile）并裁剪

lake <- vect("hydro_p_LakeErie.shp")
rc <- crop(r, lake)
rc 

#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#min values  :         NaN,         NaN,         NaN,         NaN,     271.170,       0.283, ... 
#max values  :         NaN,         NaN,         NaN,         NaN,     277.090,       0.622, ... 
#time        : 2019-01-01 
 

这样保存到GTiff文件会很方便（或者在crop中使用文件名参数更好）

gtf <- writeRaster(rc, "test.tif", overwrite=TRUE)
gtf
#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : test.tif 
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#min values  :         NaN,         NaN,         NaN,         NaN,     271.170,       0.283, ... 
#max values  :         NaN,         NaN,         NaN,         NaN,     277.090,       0.622, ... 

改变的是数据现在在文件中，而不是在内存中。而且您仍然拥有层（变量）名称。

要将图层作为变量写入 NetCDF 文件，您需要创建一个SpatRasterDataset。你可以这样做：

x <- as.list(rc)
s <- sds(x)
names(s) <- names(rc)
longnames(s) <- longnames(r)
units(s) <- units(r)

注意longnames(r) 和units(r) 的使用（不是rc）。这是因为r 有子数据集（每个都有一个长名称和一个单位），而rc 没有。

现在使用writeCDF

z <- writeCDF(s, "test.nc", overwrite=TRUE)
 
rc2 <- rast("test.nc")
rc2

#class       : SpatRaster 
#dimensions  : 182, 555, 52  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -83.475, -78.85, 41.38333, 42.9  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : test.nc:water_surface_height_above_reference_datum  
              test.nc:water_surface_height_uncertainty  
              test.nc:lake_surface_water_extent  
              ... and 49 more source(s)
#varnames    : water_surface_height_above_reference_datum (water surface height above geoid) 
              water_surface_height_uncertainty (water surface height uncertainty) 
              lake_surface_water_extent (Lake Water Extent) 
              ...
#names       : water~datum, water~ainty, lake_~xtent, lake_~ainty, lake_~ature, lswt_~ainty, ... 
#unit        :           m,           m,         km2,         km2,      Kelvin,      Kelvin, ... 
#time        : 2019-01-01 

所以看起来我们有一个具有相同结构的 NetCDF。

请注意，如果只有一个时间步，terra 的当前 CRAN 版本会删除 time 变量。开发版（1.3-11）保持时间维度，即使只有一步。

您可以安装开发版本 install.packages('terra', repos='https://rspatial.r-universe.dev')