Querying Elevation Data in Postgis

Question

Problem: I have a bunch of raster data (DTED) and I want to use this to make elevation queries using PostGIS (if it makes sense). E.g. I want to ask the question "what is the elevation at (lat1,lon1), (lat2, lon2)?" and the database would return the interpolated elevations based on the nearest neighbors where I have data.

• Idea 1: Should I load the raster data into the data base in its raw form and make queries on that? Immediate downside sounds like the rasters are pretty big and maybe querying would be slower?

• Idea 2: should I preprocess the raster data and load a bunch of points into the database (i.e. the database would have columns latitude, longitude, and altitude for all the measurements in the DTED data).

An example is show below.The intersection of the grid represents points that I have, the red stars represent points that I am querying for. So ideally, I would get a bilinear interpolation (or some other method) between my nearest neighbors. In the comments below, it was mentioned that the function http://postgis.net/docs/manual-dev/RT_ST_InvDistWeight4ma.html would help with this calculation. However, it is not clear to me how to use this function to get the answer I am looking for.

Answer 1

Your elevation points are definitely regularly aligned and hence there is no need to bother about nearest neighbors. The value of the pixel in which your point fall into is always the value of the nearest neighbor.

1. Load all the raster in one command like this:

   raster2pgsql -t 10x10 -I -C c:/temp/*.dem schema.table | psql -U postgresuser -d database

This will load all the rasters into one table, tile them as 10x10 pixel raster per row and index all those rows.

2. Query the raster for each point like this:

   SELECT pointid, ST_Value(rast, ST_SetSRID(ST_MakePoint(x, y), 12345)) val FROM rastertable, pointtable WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(x, y), 12345));

ST_SetSRID() is to build your points in the same coordinate system as the rasters. Just do:

SELECT ST_SRID(rast) FROM rastertable LIMIT 1;

to know the right number of the coordinate system and replace '12345' with that number.

ST_Intersects(rast, geom) takes avantage of the indexes and make sure the query will return quite fast.

You might be very unlucky and some of your points will fall exactly on the edge of two tiles. In that case the query will return more rows than the number of points. You can deal with that by GROUPING BY pointid and compute an average of the two values or by adding DISTINCT ON (pointid) to the query to select one over the other.