return only one value per point from ST_Value() - currently returning 2 values for some points

Question

I'm testing a point sample query with a subset of 5 points on a single band elevation raster. Using my usual approach (which I've previously not noticed this behavior), I get six records returned:

SELECT t1.pt_id, CEILING(ST_Value(rast, geom)) as elev_m
FROM (
    SELECT pt_id, geom
    FROM test.pt
    LIMIT 5
) t1
LEFT JOIN test.raster t2 ON ST_Intersects(t1.geom, rast);

Output:

It's clear which record is wrong. I understand that the offending point intersects with two pixels (of possibly two different tiles) and therefore, is returning two values.

I temporarily fixed this particular problem, by eliminating the duplicate pt_id and the value of 0 with:

SELECT t1.pt_id, CEILING(t2.elev_m)
FROM (
    SELECT pt_id, geom
    FROM test.pt
    LIMIT 5
) t1
JOIN LATERAL (
    SELECT t1.pt_id, MAX(ST_Value(rast, t1.geom)) as elev_m
    FROM test.raster 
    WHERE ST_Intersects(t1.geom, rast)
    GROUP BY t1.pt_id
) t2 ON TRUE;

But this solution will not be reusable necessarily, because other rasters may have outliers that are high, not low like in this example. Also, I'm not concerned about a few outliers in my dataset. I'm more concerned with returning only raster value per point.

So, how can I return one value per point, preferably the value of the pixel that it intersects 'most' with? I need a general solution that can be re-used consistently across different datasets. Also, efficiency is key, as this will be scaled up to millions of points once I get this worked out.

I found this question, but its answers remain unaccepted and did not help me for this particular problem.

Again, my question is: How can I create a reusable query that will ensure that only one raster value is returned per point?

Answer 1

Thanks to John Powell's question and this bug report, I've learned that I'm not the only one who has experienced similar behavior. My current solution for this problem is basically the same work-around described in my question, with the substitution of ST_NearestValue() instead of ST_Value():

SELECT t1.pt_id, CEILING(t2.elev_m)
FROM (
    SELECT pt_id, geom
    FROM test.pt
    LIMIT 5
) t1
JOIN LATERAL (
    SELECT t1.pt_id, MAX(ST_NearestValue(rast, t1.geom)) as elev_m
    FROM test.raster 
    WHERE ST_Intersects(t1.geom, rast)
    GROUP BY t1.pt_id
) t2 ON TRUE;

I made this substitution because I did not previously know about this function, and seems like a good choice instead of ST_Value() in my case because it will choose the closest non-NoData value instead of returning NoData values. However, that is not the solution to my problem, which is why I simply incorporated it into my work-around described in the question.

I should note that although I was able to obtain the desired result by incorporating MAX(ST_NearestValue()) into the LEFT JOIN with a GROUP BY clause, I chose to stick with the LATERAL JOIN approach (for now) because in my larger view definition the performance hit was negligible and the code is shorter and cleaner.

Hopefully in the future either I or someone else will be able to describe an actual solution, not some work-around.