CRS is lost after writing nodata value

[cordmaur]

After a call to write_nodata method, the CRS info is lost.

Here is a sample code:

print('Printing variables for `vul` Array')
print(f'CRS={str(vul.rio.crs)}')
print(f'NO DATA={str(vul.rio.nodata)}')

print('After writing nodata')
urban_vul = vul.rio.write_nodata(0)
print(f'CRS={str(urban_vul.rio.crs)}')
print(f'NO DATA={str(urban_vul.rio.nodata)}')

Output:
Printing variables for `vul` Array
CRS=EPSG:4326
NO DATA=None
After writing nodata
CRS=None
NO DATA=0

Expected output would be CRS=EPSG:4326

Below is the environment information. Versions rioxarray (0.14.1) deps: rasterio: 1.3.7 xarray: 2023.5.0 GDAL: 3.6.4 GEOS: 3.11.1 PROJ: 9.0.1 PROJ DATA: /usr/local/lib/python3.10/dist-packages/rasterio/proj_data GDAL DATA: None

Other python deps: scipy: 1.11.4 pyproj: 3.6.0

System: python: 3.10.6 (main, May 29 2023, 11:10:38) [GCC 11.3.0] executable: /bin/python machine: Linux-5.15.133.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

[rbavery]

I'm also experiencing a loss of crs information when writing to a raster with rioxarray '0.15.1'

data_arrays[0].rio.crs
CRS.from_epsg(32612)

with xr.set_options(keep_attrs=True):
    for xarr in data_arrays:
        orig_crs = xarr.rio.crs
        xarr_attrs = xarr.attrs
        xarr = xarr.rename({'band': 'band_old'})
        xarr = xarr.stack(band=('band_old', 'time')).transpose('band', 'y', 'x')
        xarr = xarr.rio.update_attrs(xarr_attrs)
        xarr = xarr.rio.update_attrs(new_attrs = {"count": 36})
        xarr.rio.write_crs(orig_crs, inplace=True)
        xarr.rio.to_raster(f"{outdir}/{str(xarr.coords['s2:granule_id'][0].values)}.tif", driver="COG")

path = "solar_test_scenes/S2B_OPER_MSI_L2A_TL_ESRI_20201204T131109_A019558_T12SUC_N02.12.tif"
local_path = "./solar_test_scenes/S2A_OPER_MSI_L2A_TL_ESRI_20201129T020734_A028395_T12STB_N02.12.tif"
testxarr = rioxarray.open_rasterio(path)

testxarr.rio.crs == None
True

[snowman2]

Mind providing a simple, reproducible example: https://stackoverflow.com/help/minimal-reproducible-example

[cordmaur]

Here is a minimal example, using an empty array. I divided it in 3 steps to keep track of crs and no_data.

1- Create an empty array with xarray and checking CRS and NO_DATA:

empty_array = xr.DataArray()
print(f'CRS={empty_array.rio.crs}')
print(f'NO DATA={empty_array.rio.nodata}')
--------
empty_array:
CRS=None
NO DATA=None

2- Setting CRS='epsg:4326' :

new_array = empty_array.rio.set_crs('epsg:4326')
print(f'CRS={new_array.rio.crs}')
print(f'NO DATA={new_array.rio.nodata}')
--------
CRS=EPSG:4326
NO DATA=None

3- Setting NO_DATA=0:

nodata_array = new_array.rio.write_nodata(0)
print(f'NO DATA={nodata_array.rio.nodata}')
print(f'CRS={nodata_array.rio.crs}')
--------
NO DATA=0.0
CRS=None

Here is the environment provided by show_versions():

rioxarray (0.15.1) deps:
  rasterio: 1.3.9
    xarray: 2024.2.0
      GDAL: 3.6.4
      GEOS: 3.11.1
      PROJ: 9.0.1
 PROJ DATA: /usr/local/lib/python3.10/dist-packages/rasterio/proj_data
 GDAL DATA: /usr/local/lib/python3.10/dist-packages/rasterio/gdal_data

Other python deps:
     scipy: 1.11.4
    pyproj: 3.6.1

System:
    python: 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]
executable: /usr/bin/python3
   machine: Linux-6.1.58+-x86_64-with-glibc2.35

[snowman2]

Related #356

Use rio.write_crs not rio.set_crs.

[cordmaur]

It works now. Closing the issue. Thx.