rio.ro_raster changes the coordinate

[hamiddashti]

rio.to_raster slightly changes the coordinate. Here is a small netcdf file which I used in the following code to reproduce the issue:

File (13kb)

xr.__version__
'0.17.0'
rioxarray.__version__
'0.3.2.dev0'

nc = xr.open_dataarray("nc.nc") # open netcdf
nc.rio.to_raster("test_tif.tif") # save it as raster
tif = xr.open_rasterio("test_tif.tif") # open the raster file
nc.y.values==tif.y.values  # compare y (lat) values
>>>array([ True, False,  True])  # as you can see they are different

# If you look at values they look very similar
nc.y.values
>>> array([3148842.63092262, 3147916.07989068, 3146989.52885875])
tif.y.values
>>> array([3148842.63092262, 3147916.07989068, 3146989.52885875])

# But a closer look at the "y" that is different between nc and tif shows a tiny change in the last digit
nc.y.values[1]
>>> 3147916.0798906838

tif.y.values[1]
3147916.0798906833

# same thing for x values: 
nc.x.values==tif.x.values
>>> array([ True, False, False,  True]

# An example problem that comes from this issue is we cannot do simple calculations:

tif-nc
<xarray.DataArray (band: 1, y: 2, x: 2)>
array([[[0., 0.],
        [0., 0.]]])
Coordinates:
  * y        (y) float64 3.149e+06 3.147e+06   
  * x        (x) float64 -1.878e+06 -1.875e+06
 
# Note coordinates dropped. 

This might lead to other issues. For example I recently asked a question on stackoverflow using xarray.where (here) which I think is related to this problem. Following suggestion in #298 I just reassign the coordinates to solve the problem...

I wonder is this a bug or I am doing something wrong here?

[snowman2]

I wonder is this a bug or I am doing something wrong here?

It has do do with rounding error in floating point aritmetic ref. In the geotiff, it does not store the raw coordinate values, instead it stores the transform. The transform is calculated based on the x/y coordinates in the file and written to the geotiff.

There are two places the rounding error could happen:

1. When calculating the transform to store in the geotiff
2. When generating the coordinates based on the transform in the geotiff

I think it would be nice for xarray to have a "close enough" option where a certain number of decimal places is enough to consider them equal. Not sure the best place to report this, but you could add your use case here.

[hamiddashti]

Thank you, that answers my question.