Point sampling raster with PostGIS

Question

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?

score 4 · Accepted Answer · edited Nov 28 '15 at 11:12

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);

edited Nov 28 '15 at 11:12

Alexandre Neto

14,212
2
58
85

answered Nov 25 '15 at 17:34

sokimsu

56
2

So I need to convert the mulltipoint to point for the second query to work? – Cliff Nov 25 '15 at 17:51
1

Yep, ST_Value(r.rast, 3, (ST_Dump(p.geom)).geom). – sokimsu Nov 25 '15 at 17:58
Wow! Thanks so much for your help! I will post the solution below. – Cliff Nov 25 '15 at 18:01

score 2 · Answer 2 · answered Nov 25 '15 at 18:11

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);

Point sampling raster with PostGIS

2 Answers2

Linked