Unexpected layer name when reading NetCDF file using terra and raster R packages

27 Views Asked by At

I need to read a NetCDF file in R using terra package. Here is a snapshot of the 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

Reading the file using raster:

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 

The name of the layer is the long_name attribute of the Var variable, not the variable name itself. Is it possible to directly force raster to assign the name of the layer to the variable of interest (i.e., without using names(R) = "Var")?

names(raster::raster("Test11.nc", var = "Var"))
# "Map1.Description"

When reading it by terra the name of the layer is incorrect.

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"

The name of the layer is the name of the variable concatenated with the CRS dimension lambert_azimuthal_equal_area followed by "=1".

Which attributes are read by default by terra and raster? I tried explicitly assigning attributes like varname or name to the Var variable without affecting how terra reads layer name. Any explanation or suggestions?

0

There are 0 best solutions below