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

[Phil Hurvitz]

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
**************************************************************