5

I have about 1 mill 2D points in a table in a PostgreSQL database.

I also have a DEM (Digital Elevation Model) in raster. It’s a lot of .tiff files, bound together by a .vrt file, so they function as one raster data set, at least when seen from QGIS.

I’m looking for the best way to attach to each point, the height in the DEM, at the points location. I would prefer to have that point in a new field in the point table (i.e. as a new attribute to the point feature).

After looking at the logical (PostGIS) possibilities, ST_Value [http://postgis.net/docs/RT_ST_Value.html] seems to be the relevant function. But I can’t see any ways to point the ST_Value() to a .vrt file. And since the DEM is about 1Tb I would like to avoid pulling it inside the PostgreSQL database. I was briefly considering FDW (Foreign Data Wrapper) but have no experience in that direction.

Googling around I found GDAL’s gdallocationinfo [http://www.gdal.org/gdallocationinfo.html]. So I generated a text file with the 1 mill x y coordinates of the points, and copy-pasted that list into a batch-file with lines like this.

gdallocationinfo -geoloc -valonly my_dem.vrt 444302.48 6157429.94 > collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444322.773687121 6156609.71326993 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444514.567093169 6209236.39142913 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444516.550364779 6209265.20736336 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444523.113049805 6209178.98623131 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444528.996558214 6209137.25582758 >> collect_dtm_values.values
...

Though this seems to work, it will estimated take 3 days to run the 1 mill lines batch file (it’s running).

A part from not being an elegant solution, it also leaves me with just a simple text file, holding one number per line. I would now have to figure out a good way to connect the values back onto the points, and the only connection is the fact that their order of appearance in the output text file, is the same as the order of appearance in the input batch file.

I will have to do it again and again, so I’m willing to invest in a better solution :-)

So – the question is: What would be the good way to sample points in a .vrt raster, given a set of points in a PostgreSQL table?

I work in PostgreSQL/PostGIS because I find it fast and convenient, so a solution in this environment would be preferable, and in any case, a solution that runs considerable faster, than what I have now, would be nice.

  • I don't know much about vrt files, but you can load rasters to Postgis using raster2pgsql, which is a wrapper on GDAL functions, using the -R (out-db) switch, ie, the meta data is loaded, but the rasters stay on the file system. This should allow you to use RT_ST_Value (and other functions) without having to load 1Tb of data. I don't know how out-db raster performs in comparison to in-db -- I have heard it said much slower. It should be easy enough to test with a small subset. – John Powell Mar 29 '17 at 13:31
  • Thanks John Barca. I tried your suggestion, but it seems that the -R option have no effect?!? It still produces .sql file with a lot of data insert lines. Any idea if the -R flag is unstable in any ways? – MartinHvidberg Mar 30 '17 at 9:08
  • For reference... I do this: "C:\Program Files\PostgreSQL\9.6\bin\raster2pgsql.exe" -c -t 256x256 -R -I -T ts_raster F:\GDB\DTM.vrt > D:\temp\load_DTM2pg.sql – MartinHvidberg Mar 30 '17 at 9:11
  • Sorry - I know i forgot the schema.table. The real command is "C:\Program Files\PostgreSQL\9.6\bin\raster2pgsql.exe" -c -t 256x256 -R -I -T ts_raster F:\GDB\DTM.vrt dhm.dtm > D:\temp\load_DTM2pg.sql But -R is still doing nothing. – MartinHvidberg Mar 30 '17 at 9:19
  • The sql file would be the meta data, no? Can you post a sample.? – John Powell Mar 30 '17 at 11:11
0

Working with POSTGIS out-of-db rasters, try use full file path in both the VRT construction and when loading with raster2pgsql.

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.