4

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.

enter image description here

  • You have the postgis function InvDistWeight4 which I think does what you want. You can chop a raster up into tiles essentially, which will make processing faster. There is no reason to assume that raster data is bigger than storing lat/lon points plus attributes in a table, as a raster stores the meta data, ie, start coordinate, pixel size, etc and all the points are just offsets from that. From your picture, it would seem that your data are not regularly spaced though, in which case you might want to consider PDAL – John Powell Nov 17 '16 at 21:00
  • I confess I am no expert in either of these and someone else might have a more coherent answer. – John Powell Nov 17 '16 at 21:02
  • I believe the data is spaced on a regular grid (I am using DTED), so I think I should update the picture. I'm just having a hard time identifying how I should load the data into the database, i.e. should I just load the rasters in, or should I preprocess the rasters and just extract the lat,lon,alt using GDAL and upload that metadata to the database. – user985030 Nov 17 '16 at 21:05
  • Use raster2pgsql (which uses gdal under the hood). But make sure you tile it, as this will lead to much faster calculations, when you are only interested in interpolation in small areas. – John Powell Nov 17 '16 at 21:09
  • Thanks! I did indeed load it into my database that way (I wasn't really sure what tile size to use though?). But I have heard rumors that loading the raster data into the database that way can make queries very slow, hence why I was looking at possibly loading in only lat,lon,alt points using gdal2xyz.py. So currently I have a subset of the total raster data in the database, but I have no idea how to query for the elevation of interpolated points. I'll look into the function you mentioned above "InvDistWeght4". – user985030 Nov 17 '16 at 21:14
7

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, y1), 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.

  • Thanks for information. So I have successfully loaded the raster data into the database. I tried running the code in #3 but didn't have any success because I don't have pointtable. Do I have to create a table of points just to make a query? This seems odd to me... – user985030 Nov 17 '16 at 22:35
  • Why would you not want to do some kind of IDW if your point of interest falls between grid points (and you are on a slope)? – John Powell Nov 18 '16 at 11:08
  • @user985030 Right! I Removed query 2 and edited query 3 (now 2) to better fit your needs. – Pierre Racine Nov 18 '16 at 22:17
  • @JohnBarça That's a good question. I would say 1) ST_InvDistWeight4ma is misleading and was not developed to interpolate between points. It was developed to compute missing pixel values in a raster. 2) There is nothing in PostGIS to do interpolation between points. You could probably develop a ST_ValueInterpolate(rast, geom, method) that would compute an interpolated value based on the 9 or 25 pixels surrounding the point. – Pierre Racine Nov 18 '16 at 22:26
  • Pierre, thank you, that is useful to know. To be honest, for anything on a grid, I would probably use numpy (and related) these days, but I continue to monitor postgis raster with interest. – John Powell Nov 18 '16 at 22:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.