Problems with coordinates in NetCDF using rioxarray.open_rasterio (NotGeoreferencedWarning)

[jamaa]

I am having problems processing NetCDF files.

This a sample file: https://opendata.dwd.de/climate_environment/CDC/grids_germany/monthly/hyras_de/precipitation/pr_hyras_1_2022_v3-0_de_monsum.nc

It uses EPSG:3034 as the coordinate system.

This is the code I am using to open the file:

import xarray as xr
import rasterio
import rioxarray
file = r"path\to\downloaded\file\linked\above"
xds = rioxarray.open_rasterio(file, parse_coordinates=True, variable="pr", masked=True)

This produces the following warning: "NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned."

The read x and y coordinates are not correct, as can be seen in the following screenshot (they are just cell indices offset by half a cell):

My goal is to use xds.rio.clip() but that always fails with the message "NoDataInBounds: No data found in bounds.", I presume because the coordinates were not read/interpreted correctly.

I am using

• rioxarray 0.12.4
• rasterio 1.3.3
• xarray 2022.11.0
• gdal 3.4.3

on Windows

This used to work just fine for me, but unfortunately after some software problems I recently had to reinstall everything from scratch. Unfortunately, I do not know which package versions I was working with when it worked.

So i presume something has changed in either rioxarray or some of its dependencies that is now causing this problem.

I saw here #454 (reply in thread) that someone suggested using xarray.open_dataset() instead of rioxarray.open_rasterio() and that actually solves my problem. But why does rioxarray.open_rasterio() not work as it used to before?

[snowman2]

I believe it is related to:

• BUG: Affine missing for netcdf file OSGeo/gdal#4075
• TST: Tests fail with GDAL 3.3.1 #371

If you use rio.write_coordinate_system after opening with xarray.open_dataset and write to a netCDF, it should work again.
See: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions

[jamaa]

Thanks for the reply, but I'm afraid it doesn't work for me. This is what I tried:

ds = xr.open_dataset("pr_hyras_1_2022_v3-0_de_monsum.nc")
ds.rio.write_coordinate_system(inplace=True)
ds.to_netcdf("temp.nc", mode="w")

xds = rioxarray.open_rasterio("temp.nc", parse_coordinates=True, variable="pr", masked=True)

This still results in the same warning "NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned."

If I add the following before saving to temp.nc:

ds.rio.write_crs("epsg:3034", inplace=True)

then I get the following error from .to_netcdf: "ValueError: failed to prevent overwriting existing key grid_mapping in attrs. This is probably an encoding field used by xarray to describe how a variable is serialized. To proceed, remove this key from the variable's attributes manually."

The bugs you linked above are both marked as fixed, does that mean this problem should not be occurring anymore in the most recent versions?

[snowman2]

It seems the root of the issue is rasterio/GDAL is unable to determine the transform:

>>> rds = rasterio.open("netcdf:pr_hyras_1_2022_v3-0_de_monsum.nc:pr")
>>> rds.crs
CRS.from_wkt('PROJCS["ETRS89 / LCC Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980 ",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",35],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_origin",52],PARAMETER["central_meridian",10],PARAMETER["false_easting",4000000],PARAMETER["false_northing",2800000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3034"]]')
>>> rds.transform
Affine(1.0, 0.0, 0.0,
       0.0, 1.0, 0.0)

[snowman2]

>>> from osgeo import gdal
>>> ds = gdal.Open("netcdf:pr_hyras_1_2022_v3-0_de_monsum.nc:pr")
>>> ds.GetGeoTransform()
(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

[snowman2]

Forwarded: OSGeo/gdal#6656

[jamaa]

thanks for looking into this

[RichardScottOZ]

Yes, this is pretty common with netcdf origin things in my experience - so a gdal fix for this in a new version of gdal Alan?

[snowman2]

The fix is tagged for 3.6.1: OSGeo/gdal#6663