[postgis-users] Raster Points in Polygon (was Re: PostGIS Raster for Met Data)

[Pierre Racine]

Hi,

I had to modify the coordinates of your raster and your polygon because the 900917 SRID do not exist by default in spatial_ref_sys, but this example returns the x and y coordinates (a well as the point geometry and the value) of every pixels intersecting a polygon.

I had to create the ST_PixelAsPoints(rast raster, band integer) function which is a derivate of the ST_PixelAsPolygons(rast raster, band integer) function available in this page: http://trac.osgeo.org/postgis/wiki/WKTRasterUsefulFunctions

It takes 60 seconds on my machine. I display everything in OpenJump.

CREATE TYPE geomvalxy AS (
    geom geometry,
    val double precision,
    x int,
    y int
);

-----------------------------------------------------------------------
-- ST_PixelAsPoints
-- Return all the pixels of a raster as a geomval record
-- Should be called like this:
-- SELECT (gv).geom, (gv).val FROM (SELECT ST_PixelAsPoints(rast) gv FROM mytable) foo
-----------------------------------------------------------------------
DROP FUNCTION IF EXISTS ST_PixelAsPoints(rast raster, band integer);
CREATE OR REPLACE FUNCTION ST_PixelAsPoints(rast raster, band integer) 
    RETURNS SETOF geomvalxy AS 
    $$
    DECLARE
        rast alias for $1;
        w integer;
        h integer;
        x integer;
        y integer;
        result geomvalxy;
    BEGIN
        SELECT st_width(rast), st_height(rast)
        INTO w, h;
        FOR x IN 1..w LOOP
             FOR y IN 1..h LOOP
                 SELECT ST_Centroid(ST_PixelAsPolygon(rast, x, y)), ST_Value(rast, band, x, y), x, y INTO result;
            RETURN NEXT result;
         END LOOP;
        END LOOP;
        RETURN;
    END;
    $$
    LANGUAGE 'plpgsql';

DROP FUNCTION IF EXISTS ST_PixelAsPoints(rast raster);
CREATE FUNCTION ST_PixelAsPoints(raster) RETURNS SETOF geomvalxy AS 
$$
    SELECT ST_PixelAsPoints($1, 1);
$$
LANGUAGE SQL;

-- Display raster extent in OpenJump
SELECT ST_AsBinary(ST_ConvexHull(ST_AddBand(ST_MakeEmptyRaster(850, 1010, -6.7390, -16.0390, 0.036, 0.036, 0, 0, 900913), '1BB')))

-- Display the polygon in OpenJump
SELECT ST_AsBinary(ST_GeomFromText('Polygon((-6 -10, 0 -10, 0 -15, -6 -15, -6 -10))', 900913))

-- Display raster as points in OpenJump
SELECT ST_AsBinary((gvxy).geom), (gvxy).val, (gvxy).x, (gvxy).y
FROM (SELECT ST_PixelAsPoints(ST_AddBand(ST_MakeEmptyRaster(850, 1010, -6.7390, -16.0390, 0.036, 0.036, 0, 0, 900913), '1BB')) gvxy) foo

-- Display the intersection of the raster as points with the geometry in OpenJump
SELECT ST_AsBinary((gvxy).geom), (gvxy).val, (gvxy).x, (gvxy).y
FROM ST_PixelAsPoints(ST_AddBand(ST_MakeEmptyRaster(850, 1010, -6.7390, -16.0390, 0.036, 0.036, 0, 0, 900913), '1BB')) gvxy, 
     ST_GeomFromText('Polygon((-6 -10, 0 -10, 0 -15, -6 -15, -6 -10))', 900913) geomin
WHERE ST_Intersects((gvxy).geom, geomin)

Pierre

> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
> bounces at postgis.refractions.net] On Behalf Of Michael Akinde
> Sent: Wednesday, August 24, 2011 9:46 AM
> To: PostGIS Users Discussion
> Subject: [postgis-users] Raster Points in Polygon (was Re: PostGIS Raster for Met
> Data)
> 
> Hi,
> 
> Vacation time's over, and I'm experimenting with this again.
> 
> ----- Original Message -----
> > > Not quite - I was thinking more of World2RasterCoord, but using a
> > > geometry to extract multiple raster points. Essentially,
> > > ST_INTERSECT, but returning raster coordinates rather than world
> coordinates and value.
> >
> > You can extract every point of a complex geometry using a mixture of
> > ST_PointN(),ST_NPoints(),ST_NumGeometries() and ST_GeometryN(). Have a
> > look in the PostGIS doc.
> 
> I'm probably not getting across exactly what I want to do. I'll try to give an
> example:
> 
> Given a raster rast, e.g.,
> st_makeemptyraster( 850, 1100, -6.7390, -16.0390, 0.036, 0.036, 0, 0, 900917 )
> 
> I would like to query the raster using a geom, e.g.,
> ST_GeomFromText('Polygon((10 60, 11 60, 11 61, 10 61, 10 60 ))', 4030)
> transformed to the appropriate SRID (900917 in this case) and return the result
> as a set of rows with RasterCoordX/RasterCoordY for each raster point in the
> polygon.
> 
> As an alternative, retrieving the set of point values stored in the raster would
> work if we could extract as POINT/pixel-value pairs. I suppose I may have stared
> myself blind on this problem, but I haven't found any obvious way to do this
> operation. Is it doable in any form of reasonable manner, and if so - can
> someone point me in the relevant direction?
> 
> Regards,
> 
> Michael A.
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users