Physrisk pixel is area

[EglantineGiraud]

Hi Joe,

“_linear_interp_frac_coordinates” is correct given that it is a basic bilinear interpolation on a bidimensional array (so it should be agnostic to the GPS coordinates).

I think the issue is rather in “_get_coordinates” which should use by default the “PixelIsArea” mode. Indeed, when looking at the coordinates in the data set, they seems to always follows this convention, that is, the latitudes/longitudes are built as in the code snippet below extracted from xarray.open_rasterio:

x, _ = riods.transform * (np.arange(nx) + 0.5, np.zeros(nx) + 0.5)
_, y = riods.transform * (np.zeros(ny) + 0.5, np.arange(ny) + 0.5)

Hence, the algorithm finding the pixel from latitude/longitude should correspond to the inverse of this function, that is:
~transform * coordinates - 0.5

The interpolation will still be on the edges but the pixel to interpolate will be the green point instead of the red point:

Note that it will also impact the “floor” mode. The convention applied to any Zarr store mocker also needs to be updated to apply the 0.5-pixel offset.

[joemoorhouse]

Hi @EglantineGiraud,

Thanks for looking at that! I like the explicit use of pixel_is_area and agree (of course) that we are using xarray conventions:

x, _ = riods.transform * (np.arange(nx) + 0.5, np.zeros(nx) + 0.5)
_, y = riods.transform * (np.zeros(ny) + 0.5, np.arange(ny) + 0.5)

with inverse:
~transform * coordinates - 0.5

This is indeed in the AffineTransformer class of rasterio.
https://rasterio.readthedocs.io/en/stable/topics/transforms.html

One thing I suggest we change. I still think that floor behaviour is unchanged, but that's because 'floor' has a particular meaning by convention which appears to be the intuitive one: 'pick the pixel that the lat/lon falls in when pixel is projected into lat/lon space'. This means that as a pixel value it is not simply taking (~transform * coordinates - 0.5) and applying the mathematical floor.

I found that surprising following your logic above, so checked the AffineTransformer

https://github.com/rasterio/rasterio/blob/a5e8fd02f11b6ea7ea6fc428e59403c2cd31b3cf/rasterio/transform.py#L476

def xy(self, rows, cols, zs=None, offset='center'):
at line 396 it applies the 0.5, 0.5 offset - all makes sense

But for reversing this to get the pixel under the floor convention at line 347:
def rowcol(self, xs, ys, zs=None, op=math.floor, precision=None):
there is no 0.5, 0.5 offset.

So I think we change _get_coordinates as you suggest, but in the floor case we add back on the 0.5 pixel offset before truncating to keep it in alignment with that convention?

Thanks,
Joe