[postgis-users] extract a set of WKT raster values from a point geometry table

[Pierre Racine]

Phil,

Glad you came with this question... No you're not the first and only one one to ask this question. Actually WKT Raster is specifically designed to answer this kind of need, but we are still working on it.

In a few weeks you will be able to do:

SELECT (gv).id, ST_X((gv).geom), ST_Y((gv).geom), (gv).val
FROM (SELECT pt.id, ST_Intersection(rt.rast, pt.geom) FROM yourrastertable rt, yourpointtable pt) gv

We are still working on the base to implement the ST_Intersection() function.

For now, if you compile the very last version with GDAL enabled, you can try to do something like:

SELECT pt.id, ST_X(ST_Intersection((gv).geom, pt.geom) x, ST_Y(ST_Intersection((gv).geom, pt.geom) y, (gv).val
FROM (SELECT ST_DumpAsPolygons(rast) gv FROM yourrastertable) rt, yourpointtable pt

Nevertheless I will add a ticket so we implement a ST_Value(rast, band, xcoord, ycoord) where xcoord and ycoord are expressed in the raster coordinate system.

We will also see if it is worth optimizing the future ST_Intersection() function when the vector part is a point table.

Thanks,

Pierre

>-----Original Message-----
>From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
>bounces at postgis.refractions.net] On Behalf Of Phil Hurvitz
>Sent: 13 avril 2010 22:39
>To: postgis-users at postgis.refractions.net
>Subject: [postgis-users] extract a set of WKT raster values from a point geometry table
>
>Hello all, I'm trying to get the values from a WKT raster for each point
>in a point geometry table and I get the feeling I'm going about this the
>hard way.
>
>I see how I can use the st_value function
>
>select st_value(rast, 1, x, y) from myraster;
>
>to return a single point value. One of the limitations of this method is
>that it uses pixel, rather than map, coordinates. The other is that it
>only gives a single output value.
>
>I've seen the usage for st_value at
>
><http://www.postgis.org/documentation/manual-svn/RT_ST_Value.html>
>
>That shows how to get a set of raster values from x and y sequences of
>integers, which would be very helpful if the x and y sequences were (a)
>image coordinates, and (b) represented the x and y values of my point data.
>
>I suppose I could kludge this into a looped stored procedure that does
>something like this, but looping is costly:
>
>get raster upper left coordinates as UL
>get raster pixel size as pxsize
>create output_table...
>for record in mypoints
>   get record(x, y)
>   transform record(x, y) to image(xt, yt) using UL and pxsize
>   select st_value(rast, 1, xt, yt) as raster_value
>   write raster_value to output_table
>
>It seems to me there should be a more efficient way (an ST_ function?)
>for this.
>
>There's also this possibility:
>
>select x, y, st_value(rast, 1, x, y) as val from myraster cross join
>(select cast((st_x(the_geom)-1220334.12)/98.43 as int) as x from
>mypoints) as x cross join (select cast((287678.357-st_y(the_geom))/98.43
>as int) as y from mypoints) as y;
>
>which seems to give me the results I'm looking for, but it still
>requires hardcoding of the UL coordinate and pixel size. Not a big deal
>for a one-off, but I'm looking at handling many point data sets and many
>rasters.
>
>I suppose I could handle those with a stored procedure that takes the
>point, raster, and output table names as arguments, reads the raster's
>UL and pixel size, then does the arithmetic, executes the select
>statment and writes out the table.
>
>I can't imagine that I'm the first person who's tried to do this;
>essentially what I'm looking for is start with a point table with a
>schema silmar to this:
>
>id | x | y
>
>and get this, extracting values from the raster:
>
>id | x | y | raster_value
>
>
>Any suggestions would be heartily welcomed!
>
>-P.
>
>**************************************************************
>     Phil Hurvitz, MFR | PhD Student, Urban Planning | CBE
>         1107 NE 45th Street, Suite 535 | Box 354802
>University of Washington, Seattle, Washington  98195-4802, USA
>phurvitz at u.washington.edu | http://gis.washington.edu/phurvitz
>"What is essential is invisible to the eye." -de Saint-Exupéry
>**************************************************************
>_______________________________________________
>postgis-users mailing list
>postgis-users at postgis.refractions.net
>http://postgis.refractions.net/mailman/listinfo/postgis-users