cog/point/ endpoint only doing 'nearest' resampling

[scottyhq]

Problem description

I'm only seeing 'nearest' resampling coming back for point sampling:

https://titiler.xyz/cog/point/-119.9019,34.3468?url=https%3A%2F%2Fcdn.proj.org%2Fus_nga_egm96_15.tif&resampling=bilinear

{
  "coordinates": [-119.9019, 34.3468],
  "values": [-36.6048774719238],
  "band_names": [
    "b1"
  ]
}

Expected Output

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdallocationinfo -geoloc -r bilinear /vsicurl/https://cdn.proj.org/us_nga_egm96_15.tif -119.9019 34.3468
Report:
  Location: (240.8924P,223.1128L)
  Band 1:
    Value: -35.9906852645825

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdallocationinfo -geoloc -r nearest /vsicurl/https://cdn.proj.org/us_nga_egm96_15.tif -119.9019 34.3468
Report:
  Location: (240P,223L)
  Band 1:
    Value: -36.6048774719238

Environment Information

titiler.xyz 0.21.1

[vincentsarago]

I think this is more a rasterio/rio-tiler design issue

👇 rio-tiler code

from rio_tiler.io import Reader

lon, lat = -119.9019, 34.3468
with Reader("https://cdn.proj.org/us_nga_egm96_15.tif") as src:
    print(src.point(lon, lat, coord_crs="epsg:4326").data)
    print(src.point(lon, lat, coord_crs="epsg:4326", resampling_method="bilinear").data)

[-36.604877]
[-36.604877]

which can be translated in rasterio code

import rasterio
from rasterio.warp import transform as transform_coords
from rasterio import windows

from rasterio.enums import Resampling

lon, lat = -119.9019, 34.3468
with rasterio.open("https://cdn.proj.org/us_nga_egm96_15.tif") as src:
    xs, ys = transform_coords("epsg:4326", src.crs, [lon], [lat])
    lon, lat = xs[0], ys[0]

    row, col = src.index(lon, lat)

    print(
        src.read(
            window=windows.Window(row_off=row, col_off=col, width=1, height=1),
            masked=True,
            resampling=Resampling["nearest"],
        )[0][0]
    )
    print(
        src.read(
            window=windows.Window(row_off=row, col_off=col, width=1, height=1),
            masked=True,
            resampling=Resampling["bilinear"],
        )[0][0]
    )

[-36.604877]
[-36.604877]

I think it comes down about how the pixel coordinates are rounded.

in the GDAL command, when using bilinear the location is (240.8924P,223.1128L) while it's 240,233 for nearest.

in rio-tiler, we will always use the rounded coordinates returned by rasterio's index() method https://github.com/rasterio/rasterio/blob/35f6fbe344b68d18f36d757531f816903cc80e28/rasterio/transform.py#L84-L135

To be honest I'm not quite sure how GDAL gdallocationinfo works but rasterio's sample selection seems different

[vincentsarago]

looking at https://github.com/OSGeo/gdal/blob/a599c264fc26c9f183c526b3a8dbba6a61e57c61/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py#L178-L200 it seems that GDAL will use a buffer to calculate value, but if nearest resampling it will use the pixel value directly https://github.com/OSGeo/gdal/blob/a599c264fc26c9f183c526b3a8dbba6a61e57c61/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py#L249-L276

This explain the difference between GDAL and Rasterio returned value. It's safe to assume that resampling will have no effect for rio-tiler/titiler point method 😓

[vincentsarago]

so GDAL does not use buffer but allow partial window read, while with rasterio/rio-tiler we use the index() method to get the pixel coordinates (as integer). So this is the main difference

[vincentsarago]

👇 I can replicate your result when using partial window

import rasterio
from rasterio.warp import transform as transform_coords
from rasterio import windows

from rasterio.enums import Resampling

row, col = 223.1128, 240.8924 # given in https://github.com/developmentseed/titiler/discussions/1104#discussion-8092649

with rasterio.open("https://cdn.proj.org/us_nga_egm96_15.tif") as src:
    print(
        src.read(
            window=windows.Window(row_off=row - 0.5, col_off=col - 0.5, width=1, height=1),
            masked=True,
            resampling=Resampling["nearest"],
        )[0][0]
    )
    print(
        src.read(
            window=windows.Window(row_off=row - 0.5, col_off=col - 0.5, width=1, height=1),
            masked=True,
            resampling=Resampling["bilinear"],
        )[0][0]
    )

[-36.604877]
[-35.990685]

[vincentsarago]

@scottyhq The comments ☝️ show there are differences between GDAL and rio-tiler/rasterio concept of point/location info. TBH I'm not sure which is the best one:

• GDAL: takes a partial pixel coordinate value and read a 1 pixel window
• Rasterio/rio-tiler: return the pixel value from the closest (round) pixel coordinates

rio-tiler implementation is similar to https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py#L280-L329 function.

Two choice here:

1. we make the rio-tiler/titiler documentation better by explaining that the point method return the pixel value at coordinates and we remove the resampling option
2. we update rio-tiler to replicate what GDAL does

[scottyhq]

Wow, thank you so much for taking the time to look into this in detail @vincentsarago ! Of the two options above, I'd advocate for number 2. The bilinear or cubic resampling really helps when trying to compare coarsely gridded data to higher resolution data.

[vincentsarago]

@scottyhq I've started an issue on rio-tiler cogeotiff/rio-tiler#793

feel free to comment there 🙏