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

[Paragon Corporation]

Phil,

We updated that page to show some examples using geo coordinates.

Can you take a look and see if the additional examples are helpful to you.

http://www.postgis.org/documentation/manual-svn/RT_ST_Value.html

Thanks,
Leo and Regina
http://www.postgis.us


-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Phil
Hurvitz
Sent: Tuesday, April 13, 2010 10:39 PM
To: postgis-users at postgis.refractions.net
Subject: [postgis-users] extract a set of WKT raster values from a
pointgeometry 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