A suspicious bug on raster.xy2ij

[liuh886]

Describe the bug

I have a short test on the raster's function of xy2ij, and ij2xy below:

x_1,y_1 = raster.ij2xy(i_1, j_1)
i_2, j_2 = raster.xy2ij(x_1, y_1)

The results show i_1, j_1 != i_2, j_2, which means the two functions are not reversible.

To Reproduce

Here is a minimal test function:

import numpy as np
import dem
dem = xdem.DEM('test_dem.tif')

# i,j = 0,0 on both cases

# case 1

x_0,y_0 = dtm.ij2xy(0, 0, offset='center')
dem.xy2ij(x_0,y_0,op=np.float32, area_or_point='Point')
-> (array([1.], dtype=float32), array([1.], dtype=float32))

# case 2

x_0,y_0 = dtm.ij2xy(0, 0, offset='ul')
dem.xy2ij(x_0,y_0,op=np.float32,area_or_point='Area')
-> (array([0], dtype=float32), array([0], dtype=float32)

Expected behavior

I expected the two functions to be reversible on both case 1 and 2

However, the 'xy2ij' has a problem of turning xy back to ij when offset='center'.

When offset='ul' and area_or_point='Area', it works as expected.

System (please complete the following information):

• OS: Window 10
• Environment: geoutils 0.0.9, xdem 0.0.7

Additional context

The main issue is how to turn the index of rows and columns into ij when area_or_point='Point'.

It seems like lines 2202/2203 at https://github.com/GlacioHack/geoutils/blob/main/geoutils/georaster/raster.py should be:

    i -= 0.5
    j -= 0.5

[rhugonnet]

Hi @liuh886,
Good catch, thanks for reporting it. We'll correct this!

I think we only recently realized that Area and Point are actually solely about interpretation (see https://gdal.org/user/raster_data_model.html#metadata) and not georeferencing. We were misguided in our first implementation by shifting the coordinate depending on the interpretation. Maybe this is because we do have to shift coordinates by half a pixel sometimes: some data do come with such a shift, and it exists also when converting to netCDF, for instance.

I'm fully back to work now. Hoping to get on that next week! (and the datum shift in xDEM too)