issues: 1516234283

What is your issue?

When decode_coords="all" is used to load a NetCDF file, references to grid mapping variables are stored in the data variables encoding. When it is saved to disk again, the grid mapping name should be added again to the data variables attributes to be CF-conform. This seems not to be the case when using to_netcdf() with extra encoding info for the respective data variable.

I will use the following NetCDF file to illustrate the issue: netcdf TL_2021_in { dimensions: x = 654 ; y = 866 ; time = 24 ; variables: float transverse_mercator ; transverse_mercator:_FillValue = NaNf ; string transverse_mercator:crs_wkt = "PROJCRS[\"DHDN / 3-degree Gauss-Kruger zone 3\",BASEGEOGCRS[\"DHDN\",DATUM[\"Deutsches Hauptdreiecksnetz\",ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4314]],CONVERSION[\"3-degree Gauss-Kruger zone 3\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",9,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",1,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",3500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"northing (X)\",north,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"easting (Y)\",east,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"unknown\"],AREA[\"Germany - West-Germany - 7.5°E to 10.5°E\"],BBOX[47.27,7.5,55.09,10.51]],ID[\"EPSG\",31467]]" ; transverse_mercator:semi_major_axis = 6377397.f ; transverse_mercator:semi_minor_axis = 6356079.f ; transverse_mercator:inverse_flattening = 299.1528f ; transverse_mercator:reference_ellipsoid_name = "Bessel 1841" ; transverse_mercator:longitude_of_prime_meridian = 0.f ; transverse_mercator:prime_meridian_name = "Greenwich" ; transverse_mercator:geographic_crs_name = "DHDN" ; transverse_mercator:horizontal_datum_name = "Deutsches Hauptdreiecksnetz" ; transverse_mercator:projected_crs_name = "DHDN / 3-degree Gauss-Kruger zone 3" ; transverse_mercator:grid_mapping_name = "transverse_mercator" ; transverse_mercator:latitude_of_projection_origin = 0.f ; transverse_mercator:longitude_of_central_meridian = 9.f ; transverse_mercator:false_easting = 3500000.f ; transverse_mercator:false_northing = 0.f ; transverse_mercator:scale_factor_at_central_meridian = 1.f ; transverse_mercator:spatial_ref = "PROJCS[\"DHDN / 3-degree Gauss-Kruger zone 3\",GEOGCS[\"DHDN\",DATUM[\"Deutsches_Hauptdreiecksnetz\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],AUTHORITY[\"EPSG\",\"6314\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4314\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"31467\"]]" ; int x(x) ; x:standard_name = "projection_x_coordinate" ; x:long_name = "x coordinate of projection" ; x:units = "m" ; int y(y) ; y:standard_name = "projection_y_coordinate" ; y:long_name = "y coordinate of projection" ; y:units = "m" ; int time(time) ; time:timezone = "UTC" ; time:units = "hours since 2020-12-31T23:00:00" ; time:calendar = "proleptic_gregorian" ; short tas(time, y, x) ; tas:_FillValue = -9999s ; tas:standard_name = "air_temperature" ; tas:long_name = "Near-Surface Air Temperature" ; tas:units = "degree_C" ; tas:cell_methods = "time: mean" ; tas:scale_factor = 0.1f ; tas:grid_mapping = "transverse_mercator" ; }

With the code ds = xr.open_dataset("TL_2021_in.nc", decode_coords="all") ds.to_netcdf("TL_2021_out.nc") the file will be written correctly to disk, the output of ncdump is the same for TL_2021_in.nc and TL_2021_out.nc.

But when I save with extra encoding information fot the data variable tas, f.e. with ds.to_netcdf("TL_2021_out.nc", encoding={"tas": {"dtype": "float32"}}) the grid_mapping attribute is missing in the file for this variable. Edit: scale_factor is also missing and _FillValue is now NaNf instead of -9999s.

I suspect this is because of the possibility that the previous encoding might become invalid when the dataset is changed or encoded differently? But this makes it quite hard not to lose the grid mapping info, which should normally not be affected by data type, fill value or similar.

I use python 3.10.5 and xarray 2022.3.0.