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

Question

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?

Any chance to make your examples above reproducible by uploading this dataset somewhere? — dimfalk
– dimfalk, Commented Mar 18, 2024 at 22:40

Patrick · Accepted Answer · 2024-09-11 19:30:31Z

0

Your data file has a problem in that the dimension "lambert_azimuthal_equal_area" should not be a dimension at all but a grid mapping variable. The variable "Var" has an attribute grid mapping that refers there. Now that it is also indicated as a dimension both raster and terra consider this a Z dimension. Since neither package works with true 3D data, any Z or T dimension is interpreted as a number of layers. The name "Var_lambert_azimuthal_equal_area=1" stands for the first layer of the Z dimension of the variable. You can only fix this by fixing your data file.

answered Sep 11, 2024 at 19:30

Patrick

33k7 gold badges73 silver badges102 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

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

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related