Hi,
I have been using PostGIS to build and query a DEM, while testing I noticed 
that I was not getting the expected results from st_value using resample => 
'nearest'. Instead of returning the value of the nearest pixel the x-coordinate 
is always simply rounded down, while the y-coordinate ist always rounded up. 
The following query showcases this behaviour:

WITH r AS ( 
        SELECT ST_SetValues( 
                ST_AddBand( 
                        ST_MakeEmptyRaster(
                        width => 2, height => 2, 
                        upperleftx => 0, upperlefty => 2, 
                        scalex => 1.0, scaley => -1.0, 
                        skewx => 0, skewy => 0, srid => 4326), 
                index => 1, pixeltype => '16BSI', 
                initialvalue => 0, nodataval => -999), 
        1,1,1, newvalueset =>ARRAY[ARRAY[10.0::float8, 50.0::float8], 
ARRAY[40.0::float8, 20.0::float8]]) AS rast 
)
SELECT
'Test 5',
round(ST_Value(rast, 1, 'SRID=4326;POINT(0 2)'::geometry, resample => 
'nearest')) as top_left,
round(ST_Value(rast, 1, 'SRID=4326;POINT(1 2)'::geometry, resample => 
'nearest')) as top_right,
round(ST_Value(rast, 1, 'SRID=4326;POINT(0 1)'::geometry, resample => 
'nearest')) as bottom_left,
round(ST_Value(rast, 1, 'SRID=4326;POINT(1 1)'::geometry, resample => 
'nearest')) as bottom_right, 
round(ST_Value(rast, 1, 'SRID=4326;POINT(0.9 1)'::geometry, resample => 
'nearest')) as should_be_bottom_right_but_is_bottom_left,
round(ST_Value(rast, 1, 'SRID=4326;POINT(1 1.1)'::geometry, resample => 
'nearest')) as should_be_bottom_right_but_is_top_right,
round(ST_Value(rast, 1, 'SRID=4326;POINT(0.9 1.1)'::geometry, resample => 
'nearest')) as should_be_bottom_right_but_is_top_left 
FROM r;

Resulting in:

?column? | top_left | top_right | bottom_left | bottom_right | 
should_be_bottom_right_but_is_bottom_left | 
should_be_bottom_right_but_is_top_right | 
should_be_bottom_right_but_is_top_left 
----------+----------+-----------+-------------+--------------+-------------------------------------------+-----------------------------------------+----------------------------------------
Test 5 | 10 | 50 | 40 | 20 | 40 | 50 | 10

When looking through the source code on git I found that the method 
'rt_band_get_pixel_resample' that is eventually called when using st_value 
simply rounds down the given raster coordinates instead of actually using 
nearest-neighbor: 

else if (resample == RT_NEAREST) { 
        return rt_band_get_pixel( 
                band, floor(xr), floor(yr), 
                r_value, r_nodata 
                ); 
        }

So my question is, is this the intended behaviour? 
Is there maybe another way to get a value from a raster using actual 
nearest-neighbor resampling?

I'm using Postgres 17.2 and postgis 3.6.0dev if that is relevant.

Thanks in advance

Attachment: smime.p7s
Description: S/MIME cryptographic signature

Reply via email to