Thursday 27 December 2018

Point sampling raster with PostGIS


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?



Answer



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

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...