4

I am developing a PostGIS query to sample rasters using a point table to detect land cover change. I have read the documentation on ST_value and I am able to get the raster value of a single point using the following query:

SELECT  rid, ST_Value(rast, 1, ST_SetSRID(ST_Point(484375.820,4742079.979), 32617)) as b1
FROM imagery.l8_2015_09_12 
WHERE ST_Intersects(rast, ST_SetSRID(ST_Point(484375.820,4742079.979), 32617)::geometry, 1);

However, when sampling using a table of points, I cannot get ST_Value to work. When running the following query in the SQL window of QGIS:

-- sample a raster from a table of points 
SELECT ST_Value(r.rast, 3, p.geom) As band3
FROM analysis_results.sample_points AS p, imagery.l8_2015_09_12 AS r
WHERE ST_Intersects(r.rast,p.geom);

I get the following error:

Attempting to get the value of a pixel with a non-point geometry

The table has a multipoint geometry data type, so why is ST_value not playing nicely?

2 Answers 2

4

ST_Value is strict about wanting a Point rather than a MultiPoint; one way to bypass this is to use ST_Dump. Taking your first query and converting the Point to a MultiPoint:

SELECT 
    rid,
    ST_Value(rast, 1, 
        (ST_Dump(
            ST_Multi(
                ST_SetSRID(ST_Point(484375.820, 4742079.979), 32617)
            )
        )).geom
    ) AS b1 
FROM imagery.l8_2015_09_12 
WHERE ST_Intersects(rast, ST_SetSRID(ST_Point(484375.820, 4742079.979), 32617)::geometry, 1);
3
  • So I need to convert the mulltipoint to point for the second query to work?
    – Cliff
    Commented Nov 25, 2015 at 17:51
  • 1
    Yep, ST_Value(r.rast, 3, (ST_Dump(p.geom)).geom).
    – sokimsu
    Commented Nov 25, 2015 at 17:58
  • Wow! Thanks so much for your help! I will post the solution below.
    – Cliff
    Commented Nov 25, 2015 at 18:01
2

The following solution works to sample a raster using a table of points.

SELECT rid, ST_Value(r.rast, 3, (ST_Dump(p.geom)).geom) As band3
FROM analysis_results.sample_points AS p, imagery.l8_2015_09_12 AS r
WHERE ST_Intersects(r.rast, ST_SetSRID(p.geom, 32617)::geometry, 1);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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