我需要使用 terra 包读取 R 中的 NetCDF 文件。这是 NetCDF 的快照
(NC <- nc_open("Test11.nc"))
File Test11.nc (NC_FORMAT_NETCDF4):
1 variables (excluding dimension variables):
float Var[x,y,lambert_azimuthal_equal_area] (Contiguous storage)
units: UNITS
_FillValue: -99999
long_name: Map1 Description
grid_mapping: lambert_azimuthal_equal_area
3 dimensions:
x Size:390
units: m
standard_name: projection_x_coordinate
long_name: x coordinate of projection
y Size:404
units: m
standard_name: projection_y_coordinate
long_name: y coordinate of projection
lambert_azimuthal_equal_area Size:1
grid_mapping_name: lambert_azimuthal_equal_area
false_easting: 4321000
false_northing: 3210000
latitude_of_projection_origin: 52
longitude_of_projection_origin: 10
long_name: CRS definition
longitude_of_prime_meridian: 0
semi_major_axis: 6378137
inverse_flattening: 298.257222101
spatial_ref: PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]
crs_wkt: PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]
GeoTransform: 2630000 10000 0 5420000 0 -10000
使用光栅读取文件:
raster::raster("Test11.nc", var = "Var")
class : RasterLayer
dimensions : 404, 390, 157560 (nrow, ncol, ncell)
resolution : 10000, 10000 (x, y)
extent : 2630000, 6530000, 1380000, 5420000 (xmin, xmax, ymin, ymax)
crs : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
source : Test11.nc
names : Map1.Description
z-value : 1
zvar : Var
图层的名称是
long_name
变量的 Var
属性,而不是变量名称本身。是否可以直接强制 raster 将图层名称分配给感兴趣的变量(即不使用 names(R) = "Var"
)?
names(raster::raster("Test11.nc", var = "Var"))
# "Map1.Description"
通过 terra 读取时,图层名称不正确。
terra::rast("Test11.nc")
class : SpatRaster
dimensions : 404, 390, 1 (nrow, ncol, nlyr)
resolution : 10000, 10000 (x, y)
extent : 2630000, 6530000, 1380000, 5420000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
source : Test11.nc
varname : Var (Map1 Description)
name : Var_lambert_azimuthal_equal_area=1
unit : UNITS
> varnames(terra::rast("Test11.nc"))
[1] "Var"
> names(terra::rast("Test11.nc"))
[1] "Var_lambert_azimuthal_equal_area=1"
> longnames(terra::rast("Test11.nc"))
[1] "Map1 Description"
层的名称是与 CRS 维度
lambert_azimuthal_equal_area
连接的变量名称,后跟“=1”。
terra 和 raster 默认读取哪些属性? 我尝试显式地将
varname
或名称等属性分配给 Var
变量,而不影响 terra 读取图层名称的方式。
有什么解释或建议吗?
您的数据文件存在问题,因为维度“lambert_azimuthal_equal_area”根本不应该是维度,而应该是网格映射变量。变量“Var”有一个指向此处的属性
grid mapping
。现在它也被表示为一个维度,raster
和 terra
都将其视为 Z 维度。由于这两个包都不能处理真正的 3D 数据,因此任何 Z 或 T 维度都被解释为多个层。名称“Var_lambert_azimuthal_equal_area=1”代表变量 Z 维度的第一层。您只能通过修复数据文件来解决此问题。