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 l...@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.herzspr...@ldbv.bayern.de> > Sent: Wednesday, February 19, 2025 2:12 PM > To: postgis-users@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 >
smime.p7s
Description: S/MIME cryptographic signature