2

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:

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?

  • 2
    Yeah, I have also had this problem, see this question: gis.stackexchange.com/questions/279949/…. A bug has been filed. My solution was always to use MAX(ST_NearestValue(... which worked for me. – John Powell Feb 14 at 17:07
  • @JohnPowell thanks for this link. I incorporated ST_NearestValue() as described in my answer, but mostly it seems I'm stuck with my work-around described in the question. – pdavis Feb 15 at 14:37
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.

  • Well, the hope, is that the bug will be fixed, but, we do at least, have a workaround. – John Powell Feb 15 at 16:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.