微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

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

如何解决将 netCDF 文件剪切为 shapefile 并在 R 中克隆元数据变量

我有 NetCDF 文件(例如 https://data.ceda.ac.uk/neodc/esacci/lakes/data/lake_products/L3S/v1.0/2019 全局域),我想根据 shapefile 边界(在本例中为湖 - 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 中较小的域/形状文件时,是否可以从原始 NetCDF 数据集中克隆/复制所有元数据变量,然后另存为 NetCDF?任何有关如何在 R 中执行此操作的指导将不胜感激。 (NetCDF 和 R 对我来说都是新手,所以我不确定我缺少什么或有深入的知识来对此进行排序)。

解决方法

您有一个包含许多 (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,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,y)
#extent      : -83.475,-78.85,41.38333,42.9  (xmin,ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       : water~datum,... 
#min values  :         NaN,NaN,271.170,0.283,... 
#max values  :         NaN,277.090,0.622,... 
#time        : 2019-01-01 
 

将其保存到这样的 GTiff 文件会很方便(或者甚至更好地在crop 中使用文件名参数)

gtf <- writeRaster(rc,"test.tif",overwrite=TRUE)
gtf
#class       : SpatRaster 
#dimensions  : 182,ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : test.tif 
#names       : water~datum,... 

改变的是数据现在在文件中,而不是在内存中。而且您仍然拥有图层(变量)名称。

要将层作为变量写入 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,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,... 
#time        : 2019-01-01 

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

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

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

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。

相关推荐


Selenium Web驱动程序和Java。元素在(x,y)点处不可单击。其他元素将获得点击?
Python-如何使用点“。” 访问字典成员?
Java 字符串是不可变的。到底是什么意思?
Java中的“ final”关键字如何工作?(我仍然可以修改对象。)
“loop:”在Java代码中。这是什么,为什么要编译?
java.lang.ClassNotFoundException:sun.jdbc.odbc.JdbcOdbcDriver发生异常。为什么?
这是用Java进行XML解析的最佳库。
Java的PriorityQueue的内置迭代器不会以任何特定顺序遍历数据结构。为什么?
如何在Java中聆听按键时移动图像。
Java“Program to an interface”。这是什么意思?