ST_VALUE not resampling using nearest neighbor

[Herzsprung, Severin (LDBV)]

Hi,

thank you for the answer.

Sadly st_nearestvalue produces the exact same results:

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',
ST_NearestValue(rast, 'SRID=4326;POINT(0 2)'::geometry) as top_left,
ST_NearestValue(rast, 'SRID=4326;POINT(1 2)'::geometry) as top_right,
ST_NearestValue(rast, 'SRID=4326;POINT(0 1)'::geometry) as bottom_left,
ST_NearestValue(rast, 'SRID=4326;POINT(1 1)'::geometry) as bottom_right, 
ST_NearestValue(rast, 'SRID=4326;POINT(0.9 1)'::geometry) as should_be_bottom_right_but_is_bottom_left,
ST_NearestValue(rast, 'SRID=4326;POINT(1 1.1)'::geometry) as should_be_bottom_right_but_is_top_right,
ST_NearestValue(rast, 'SRID=4326;POINT(0.9 1.1)'::geometry) as should_be_bottom_right_but_is_top_left 
FROM r;
 ?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


The only difference between st_value and st_nearestvalue I could ever spot was that st_nearestvalue returns the nearest value if the given point lies outside the raster.


Looking at this segment of the code for st_nearestvalue it appears that the value is returned early if a value could be found:

if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
			rt_raster_destroy(raster);
			PG_FREE_IF_COPY(pgraster, 0);
			lwgeom_free(lwgeom);
			PG_FREE_IF_COPY(geom, 2);
			elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
			PG_RETURN_NULL();
		}

		/* value at point, return value */
		if (!exclude_nodata_value || !isnodata) {
			rt_raster_destroy(raster);
			PG_FREE_IF_COPY(pgraster, 0);
			lwgeom_free(lwgeom);
			PG_FREE_IF_COPY(geom, 2);

			PG_RETURN_FLOAT8(value);
		}

But maybe I'm reading this wrong, I am far from knowledgeable in C.


Am Mittwoch, dem 19.02.2025 um 22:28 -0500 schrieb lr at pcorp.us:
> I forget why we created ST_NearestValue if ST_Value supposedly gives
> nearest already.
> 
> @pramsey remember what the difference is?
> 
> Anyway I'd suggest trying ST_NearestValue
> https://postgis.net/docs/RT_ST_NearestValue.html and see if that give
> you what you are expecting.
> 
> Hope that helps,
> Regina
> 
> -----Original Message-----
> From: Herzsprung, Severin (LDBV) <Severin.Herzsprung at ldbv.bayern.de> 
> Sent: Wednesday, February 19, 2025 2:12 PM
> To: postgis-users at lists.osgeo.org
> Subject: ST_VALUE not resampling using nearest neighbor
> 
> 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
> 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 7552 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250220/e5fb1d9d/attachment.bin>