Good practices for storing CRS in a dataset

[pierre-manchon]

Hello,

I'm encountering problems when I try to save (to_netcdf) and reopen a netCDF file with rasterio.open ; The geotransform is empty.

After some searching I found that:

• gdal looks for grid_mapping then spatial_ref and GeoTransform (link)
• rioxarray looks for grid_mapping/spatial_ref and grid_mapping/crs_wkt in the attributes but writes spatial_ref as a Coordinates (link)
• the cf-convention indicates that the CRS should be stored as a crs variable (data_var I guess) with a crs_wkt attribute which is what I do currently as you can see just after (link)

ds['crs'] = xr.DataArray(name='crs', attrs=pyproj.CRS(4326).to_cf())
ds.attrs['grid_mapping'] = 'crs'

• xarray can use grid_mapping and bounds to decode coordinates using decode_coords when opening files with open_rasterio(Add full support for xarray decode_coords corteva/rioxarray#282)

Cue to my questions:

• What is the proper way to store the CRS attributes in a netCDF file ?
• Do using spatial_ref should be prioritized against using a crs data_var (which does not work with rasterio.open ?
  Pierre

[pierre-manchon]

Reused #2288 (comment) for now