How to write CRS info to NetCDF in a way QGIS can read? (Python/Xarray)

Question

I have a NetCDF file I have subset by geometry and time, and now I need to save it as a new NetCDF file. The original file does not have the CRS set, thus I am setting it in this one. The problem is that when I load the layer in QGIS, as either a mesh layer or a raster, it doesn't find the CRS information attached to the file.

Workflow:

from xarray import open_dataset
import rioxarray from rio

# Open the dataset and add coordinate reference system
ds = open_dataset(file, decode_coords="all")
ds.rio.write_crs("epsg:4326", inplace="True")

# Clip the netCDF data to a polygon (defined elsewhere)
clipped = ds.rio.clip(geometries=gdf.geometry.values, crs=4326)

#Subset the data by time
subset = clipped.sel(time=slice('1980-01-01', '1980-12-31'))

#Calculate mean
summarized = subset.mean(dim="time", skipna=False, keep_attrs=True)

If I print summarized it looks like this:

<xarray.Dataset>
Dimensions:      (lon: 574, lat: 288)
Coordinates:
  * lon          (lon) float64 -114.0 -114.0 -113.9 ... -90.21 -90.17 -90.13
  * lat          (lat) float64 37.04 37.08 37.12 37.17 ... 48.92 48.96 49.0
    spatial_ref  int64 0
Data variables:
    prcp         (lat, lon) float32 nan nan nan nan nan ... nan nan nan nan nan
Attributes:
    title:                    SECURE Water Act 9505V3 Assessment
    creation_date:            19-Jul-2023 14:35:47

I save it as a NetCDF file:

summarized.to_netcdf(path=os.path.join(outputDir, outputName + ".nc"), format="NETCDF4", engine="netcdf4")

If I reload it, I can see the crs value is correct:

import rioxarray as rio
import xarray as xarr

ds = xarr.open_dataset(file, decode_coords="all")


print(ds.rio.crs)
//EPSG:4326


print(ds)
//<xarray.Dataset>
//Dimensions:      (lon: 574, lat: 288)
//Coordinates:
//  * lon          (lon) float64 -114.0 -114.0 -113.9 ... -90.21 -90.17 -90.13
//  * lat          (lat) float64 37.04 37.08 37.12 37.17 ... 48.92 48.96 49.0
//    spatial_ref  int64 ...
//Data variables:
//    prcp         (lat, lon) float32 ...
//Attributes:
//    title:                    SECURE Water Act 9505V3 Assessment
//    creation_date:            19-Jul-2023 14:35:47


print(ds.prcp)
//<xarray.DataArray 'prcp' (lat: 288, lon: 574)>
//[165312 values with dtype=float32]
//Coordinates:
//  * lon          (lon) float64 -114.0 -114.0 -113.9 ... -90.21 -90.17 -90.13
//  * lat          (lat) float64 37.04 37.08 37.12 37.17 ... 48.92 48.96 49.0
//    spatial_ref  int64 ...
//Attributes:
//    long_name:  annual total precipitation
//    units:      mm/year

If I load it in QGIS, it can't seem to understand the CRS:

How do I get the CRS to map correctly so QGIS can read it? One other note, I am not having this problem when I output to a geotiff:

summarized.rio.to_raster(raster_path=os.path.join(outputDir, outputName + ".tif"), driver="GTiff")

Answer 1

For a NetCDF file, you have to call rio.write_coordinate_system(..) for it to be understood by GDAL 3.3.1+.

See:

• https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions
• https://github.com/corteva/rioxarray/issues/371