[postgis-users] Retrieving Data from PostGIS Raster using ST_Value

[elliott]

Hello Everyone,

I am very new to using PostGIS Raster. I need to retrieve data fields
from a tile loaded into the database. Using ST_Height and ST_Width, I
get the number of rows and columns and iterate over these using the
ST_Value to retrieve the data. This is extremely slow taking 1/2 hour
to get all of the data required. Is there a faster way to accomplish this?

Thanks!
_______________________________________________
postgis-users mailing list
postgi...@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Nordgren, Bryce]

I'd be interested to know the size of the raster tile (width x height) which took 1/2 hour to retrieve all the info.

On the other hand, the usage paradigm for raster is to use the utility functions, because this is what is fast. Plpgsql is not compiled, so it'll be slow to iterate over all the pixels. What are you trying to do with the data once you get it? Can you accomplish the same thing by making an "expression" and using ST_MapAlgebra?

Bryce

[elliott]

It is just one SRTM3 tile which is 1201 x 1201. The data is used by
scientific codes to determine the elevation at a given latitude and
longitude.

[Nordgren, Bryce]

1.44 million pixel accesses does not sound like something that should take 1/2 hour.

For better or worse, raster depends on Postgresql's underlying caching mechanism to make data I/O efficient; and the unit of caching is one "raster tile". If I read the SRTM3 data description correctly, each pixel is a 16 bit int, making for a 2.88 Mb "raster tile" (if you loaded it all into one row.) I'm a little out of my depth here, but your tile may be too big for postgresql to cache. Does performance improve if you load the SRTM data into many raster tiles?

Bryce

-----Original Message-----
From: postgis-us...@postgis.refractions.net [mailto:postgis-us...@postgis.refractions.net] On Behalf Of elliott
Sent: Wednesday, October 12, 2011 8:28 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Retrieving Data from PostGIS Raster using ST_Value

It is just one SRTM3 tile which is 1201 x 1201. The data is used by scientific codes to determine the elevation at a given latitude and longitude.

_______________________________________________

[Paragon Corporation]

Guys -- I'd like to point out that the issue with ST_value performance I
think is the same issue we have with ST_GeometryN and ST_PointN.

That being the memory copy. For each ST_Value call the whole raster tile
needs to be copied into memory. Paul can correct me if I am mistaken, but
tht is how he explained the issue to me. This I consider a very fundamental
limitation / problem with the plpgsql / PostGIS C whatever interface which
causes a lot of wasted time since some things are just faster programming
wise to implement in plgsql than in C and its esier to take advantage of
PostgreSQL SQL power in plpgsql /sql than it is in C.

This is BTW why plpgsql for this kind of task is really slow. C doesn't
have the same issue if you are doing the steps in one call since its dealing
with the same pointer.

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

[Paragon Corporation]

Elliot,

Why are you computing this for every pixel. You should be using
ST_Intersects geometry/raster for that. That is the beauty of postgis
raster is that you can do a spatial query between the 2.

Like I mentioned in other post in response to Bryce, the performance hit in
plpgsql is the memory copy of the raster tile for each call.
If you want all, you should be using ST_Intersection or ST_DumpPolygons.

[Pierre Racine]

If you want to get all the values for every pixel do:

SELECT (ST_DumpAsPolyogns(rast)).val
FROM your table

Can I ask why you want to do that? Normally people want to extract values intersecting with a polygon or a point.

You can also get the histogram of the raster using ST_Histogram().

Pierre
________________________________________
From: postgis-us...@postgis.refractions.net [postgis-us...@postgis.refractions.net] On Behalf Of elliott [ell...@cpi.com]
Sent: Wednesday, October 12, 2011 9:04 AM

To: PostGIS Users Discussion
Subject: [postgis-users] Retrieving Data from PostGIS Raster using ST_Value

Hello Everyone,

[elliott]

The table contains SRTM3 data. What I would really like to do is
retrieve the altitude for any given lat/lon. We currently have
scientific codes that read the SRTM3 flat file into an array and
interpolate the value for the lat/lon pair. We would like to retrieve
the data from a database instead of flat file thinking the performance
would be better. However, the times retrieving the array from the
database are many orders of magnitude worse than just reading the array
in from a flat file.

[Puneet Kishor]

On Oct 13, 2011, at 8:38 AM, elliott wrote:

> We would like to retrieve the data from a database instead of flat file thinking the performance would be better.

General rule, without any metrics backing it, but backed by experience -- a database will almost never be as fast as a flat file. Esp. a flat file that holds binary data whose structure you know, and that you can access by specific offsets, will be tremendously quicker than accessing the same data held in a db. A db has way too much overhead.

You use a db not for the performance gains but for the flexibility, standardization of access, and other things you can do with the data once they are in the db. Then you throw extra hardware, memory and tune the db to bring up its performance to an acceptable point, if required.

--
Puneet Kishor

[Pierre Racine]

If you want to get the pixel value under a point just do:

SELECT ST_Value(rast, ST_Transform(ST_SetSRID(ST_Point(yourLongitude, yourLatitude), 4326), ST_SRID(rast))
FROM yourTable

That should be pretty fast.

In any case you should load your raster when loading it (-k100x100 or -k 20x20) and then index it:

CREATE INDEX srtm_tiled_rast_gist_idx ON srtm_tiled USING GIST (ST_ConvexHull(rast));

Interpolating the elevation of a point is another story. You could build a grid of 8, 15 or 24 points around your lat/long, query the elevation at those points and use your algorithm to compute the interpolation.

Really getting the elevation for 1 or 25 points should be a matter of microseconds.

Pierre

[Pierre Racine]

Please read "In any case you should TILE your raster when loading it " in the previous mail...

[elliott]

Thanks for the response.  A follow up question:  what if the yourLongitude, yourLatitude values provided in the query are not grid points in that they fall within a bounding pixel grid.  Only values for the intersection points are known so how does a non-intersecting point get assigned a value from the query?

[Pierre Racine]

> Thanks for the response. A follow up question: what if the yourLongitude,
> yourLatitude values provided in the query are not grid points in that they fall
> within a bounding pixel grid. Only values for the intersection points are known
> so how does a non-intersecting point get assigned a value from the query?

Point outside the area covered by the grid are assigned NULL. There is not inter- or extra- polation done.

I forgot to mention that if you tile your input raster then you should add a WHERE clause:

SELECT ST_Value(rast, ST_Transform(ST_SetSRID(ST_Point(yourLongitude, yourLatitude), 4326), ST_SRID(rast))
FROM yourTable

WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_Point(yourLongitude, yourLatitude), 4326), ST_SRID(rast)))

In this case all points outside the area covered by the raster will not be returned.

Pierre