2

I'm having trouble using rioxarray for converting one time layer of the .nc file to a GeoTIFF. I get the following error:

MissingSpatialDimensionError: y dimension not found. 'rio.set_spatial_dims()' or using 'rename()' to change the dimension name to 'y' can address this. Data variable: spei

MissingSpatialDimensionError: y dimension not found. 'rio.set_spatial_dims()' or using 'rename()' to change the dimension name to 'y' can address this. Data variable: spei

The nc I'm using is from https://spei.csic.es/spei_database/#map_name=spei12#map_position=1439.

My code is:

import rioxarray
import xarray

nc_fo = r'./spei12.nc'
xds = xarray.open_dataset(nc_fo)
xds.rio.write_crs("epsg:4326", inplace=True)
xds.rio.set_spatial_dims("lon", "lat", inplace=True)
xds["spei"][0].rio.to_raster('test_spei.tif')

1 Answer 1

2

I don't know why, but each operation you make on xds seems to remove the spatial dimension information, which you can check through .rio._y_dim and .rio._x_dim:

import xarray

nc_fo = r'./spei12.nc'
xds = xarray.open_dataset(nc_fo)

print(f"1) y: {xds.rio._y_dim}, x: {xds.rio._x_dim}")

xds.rio.write_crs("epsg:4326", inplace=True)

print(f"2) y: {xds.rio._y_dim}, x: {xds.rio._x_dim}")

xds.rio.set_spatial_dims("lon", "lat", inplace=True)

print(f"3) y: {xds.rio._y_dim}, x: {xds.rio._x_dim}")

xda = xds["spei"]

print(f"4) y: {xda.rio._y_dim}, x: {xda.rio._x_dim}")

xda.rio.set_spatial_dims("lon", "lat", inplace=True)

print(f"5) y: {xda.rio._y_dim}, x: {xda.rio._x_dim}")

xda_1 = xda.isel(time=0)

print(f"6) y: {xda_1.rio._y_dim}, x: {xda_1.rio._x_dim}")

xda_1.rio.set_spatial_dims("lon", "lat", inplace=True)

print(f"7) y: {xda_1.rio._y_dim}, x: {xda_1.rio._x_dim}")

xda_1.rio.to_raster('test_spei.tif')

You can see that the spatial dims are lost after selecting the variable "spei" and selecting the first time:

1) y: None, x: None
2) y: None, x: None
3) y: lat, x: lon
4) y: None, x: None
5) y: lat, x: lon
6) y: None, x: None
7) y: lat, x: lon

So it's probably safest to just rename the dims:

import xarray

nc_fo = r'./spei12.nc'
xds = xarray.open_dataset(nc_fo)
xds.rio.write_crs("epsg:4326", inplace=True)
xds = xds.rename({"lon": "x", "lat": "y"})
xds["spei"][0].rio.to_raster('test_spei.tif')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.