Xarray Dataset to raster – ValueError: failed to prevent overwriting existing key grid_mapping in attrs

python 3rasteriorioxarrayxarray

I have a netcdf file located here http://thredds.northwestknowledge.net:8080/thredds/catalog/TERRACLIMATE_ALL/data/catalog.html?dataset=TERRACLIMATE_ALL_SCAN/data/TerraClimate_tmin_1958.nc.

I am trying to convert the first time slice of this netcdf file to a raster, or geotiff.

I am doing it like so:

import pandas as pd
import xarray as xr 

#import netcdf4 file
ds = xr.open_dataset"TerraClimate_tmin_1958.nc")

#get list of all times
times = ds['time'].values

#loop through times
for i, t in enumerate(times):

    #select first time slice (January 1st, 1958)        
    in_ras = ds.sel(time=t)

    #assign crs
    in_ras.rio.write_crs("epsg:4326", inplace=True)

    #write out dataset to a raster
    in_ras['tmin'].rio.to_raster('out.tif')

but this fails with:

Traceback (most recent call last):

  File "C:\Users\user\AppData\Local\Temp\3\ipykernel_11240\3274726430.py", line 1, in <cell line: 1>
    in_ras[var].rio.to_raster(os.path.join(out, str(month) + '.tif'))

  File "C:\Users\user\Anaconda3\envs\geospatial2\lib\site-packages\rioxarray\raster_array.py", line 1100, in to_raster
    return RasterioWriter(raster_path=raster_path).to_raster(

  File "C:\Users\user\Anaconda3\envs\geospatial2\lib\site-packages\rioxarray\raster_writer.py", line 252, in to_raster
    data = encode_cf_variable(out_data).values.astype(numpy_dtype)

  File "C:\Users\user\Anaconda3\envs\geospatial2\lib\site-packages\xarray\conventions.py", line 282, in encode_cf_variable
    pop_to(var.encoding, var.attrs, attr_name)

  File "C:\Users\user\Anaconda3\envs\geospatial2\lib\site-packages\xarray\coding\variables.py", line 129, in pop_to
    safe_setitem(dest, key, value, name=name)

  File "C:\Users\user\Anaconda3\envs\geospatial2\lib\site-packages\xarray\coding\variables.py", line 112, in safe_setitem
    raise ValueError(

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.

If I do not assign the crs like this:

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

the output raster does not have a spatial reference. Before assigning it this in_ras.crs returns:

<xarray.DataArray 'crs' (crs: 1)>
array([3], dtype=int16)
Coordinates:
    time     datetime64[ns] 1958-01-01
  * crs      (crs) int16 3
Attributes:
    grid_mapping_name:            latitude_longitude
    longitude_of_prime_meridian:  0.0
    semi_major_axis:              6378137.0
    inverse_flattening:           298.257223563
    long_name:                    crs

after assigning EPSG 4326 it returns:

<xarray.DataArray 'crs' ()>
array(0)
Coordinates:
    time     datetime64[ns] 1958-01-01
    crs      int32 0
Attributes:
    crs_wkt:                      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...
    semi_major_axis:              6378137.0
    semi_minor_axis:              6356752.314245179
    inverse_flattening:           298.257223563
    reference_ellipsoid_name:     WGS 84
    longitude_of_prime_meridian:  0.0
    prime_meridian_name:          Greenwich
    geographic_crs_name:          WGS 84
    grid_mapping_name:            latitude_longitude
    spatial_ref:                  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...

I also tried this which should remove the duplicate attributes but it fails the same way:

keys = list(in_ras.crs.attrs.keys())

for k in keys:
        
     del in_ras.crs.attrs[k]
        
        
      in_ras.rio.write_crs("epsg:4326", inplace=True)
      in_ras[var].rio.to_raster('output.tif')

Best Answer

See: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html

ds = xr.open_dataset("TerraClimate_tmin_1958.nc", decode_coords="all")

With decode_coords="all", there is no need to write the CRS:

ds.rio.crs

Output:

CRS.from_wkt('GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6378137,298.257223563]],PRIMEM["undefined",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')
Related Question