Querying PostGIS Raster with MultiPoint

Question

For one of my projects using PostGIS 2 I have a table with GPS tracks (igc_files). It holds records including a MULTIPOINT() column that has about 1,000 points for each record and a column that has an array of the corresponding timestamps. I also have a raster table (elevations) that should at some point hold all of the SRTM elevation data. For now I have only imported a few relevant areas though.

What I want to do is get the elevation for each point of one certain record and the corresponding timestamps. I have implemented an naive approach and it does work, but I'm wondering about its performance implications. This is how it currently looks:

SELECT timestamps[location_id] AS timestamp,
       ST_VALUE(elevations.rast, subq.location) AS elevation 
FROM elevations, (
    SELECT (ST_DUMP(igc_files.locations)).path[1] AS location_id, 
           (ST_DUMP(igc_files.locations)).geom AS location 
    FROM igc_files 
    WHERE igc_files.id = 1
) AS subq, igc_files 
WHERE igc_files.id = 1 AND 
      ST_INTERSECTS(subq.location, elevations.rast) AND
      ST_VALUE(elevations.rast, subq.location) IS NOT NULL

Since this is a GPS track and the points are not that far apart, PostGIS could assume that the next point is in the same raster tile and check that tile first before checking the index again from the top. Is that actually done by PostGIS or is it just starting to search the index from the top for each of the points?

I have tried to use ST_INTERSECTION() instead of ST_VALUE() and feeding in the MULTIPOINT directly, but that was actually even slower by a few orders of magnitude. Is there a feature/function of PostGIS that I am missing to do this more efficiently?

Answer 1

Did you tile your raster when you loaded it in the db (using the -t raster2pgsql option)? Smaller tiles size = better performance. Check this Blog entry: http://duncanjg.wordpress.com/2012/10/30/tile-size-for-raster-vector-overlays-in-postgis/

Make sure to also create the index with the -I option...

Also in your query:

1) Don't call ST_DUMP() two times. Assign the result to a name and use (name).path[1] and (name).geom in the first SELECT

2) "ST_VALUE(elevations.rast, subq.location) IS NOT NULL" in the WHERE clause is useless as ST_INTERSECTS will return FALSE if the point intersects with a nodata value. This test is quite slow if there are a lot of nodata values in your raster. Is this the case?