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: [email protected] [mailto:postgis-users- > [email protected]] 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 > [email protected] > http://postgis.refractions.net/mailman/listinfo/postgis-users _______________________________________________ postgis-users mailing list [email protected] http://postgis.refractions.net/mailman/listinfo/postgis-users
