使用 terra 和 raster R 包读取 NetCDF 文件时出现意外的图层名称

问题描述 投票:0回答:1

我需要使用 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 读取图层名称的方式。 有什么解释或建议吗?

r netcdf r-raster terra ncdf4
1个回答
0
投票

您的数据文件存在问题,因为维度“lambert_azimuthal_equal_area”根本不应该是维度,而应该是网格映射变量。变量“Var”有一个指向此处的属性

grid mapping
。现在它也被表示为一个维度,
raster
terra
都将其视为 Z 维度。由于这两个包都不能处理真正的 3D 数据,因此任何 Z 或 T 维度都被解释为多个层。名称“Var_lambert_azimuthal_equal_area=1”代表变量 Z 维度的第一层。您只能通过修复数据文件来解决此问题。

© www.soinside.com 2019 - 2024. All rights reserved.