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. Commented Mar 29, 2017 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?
    – Martin
    Commented Mar 30, 2017 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
    – Martin
    Commented Mar 30, 2017 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.
    – Martin
    Commented Mar 30, 2017 at 9:19
  • The sql file would be the meta data, no? Can you post a sample.? Commented Mar 30, 2017 at 11:11

2 Answers 2


One solution would be to work with PostGIS Raster, using the claimed out-db way. So you don't need to store that huge dataset into Postgres, only the metadata of the raster and the paths of the files will be stored.

You should take a deeper look at the raster2pgsql documentation which support almost every file type that the GDAL driver does (including the VRT extension). An example:

$ raster2pgsql -c -s 4326 -t 256x256 -R -I -Y -C path/to/file.vrt schema.raster_table > insert_script.sql

I strongly recommend you to use the -Y param as it will make the execution WAY faster (I get to execute the script in just a few minutes for a 1TB DEM dataset).

Then, you will have to execute the generated SQL script as follows (note that since PostGIS 3.0 was released, the raster support has been moved into an extension. Ensure that your PostGIS database already supports rasters):

$ psql -U your_user -d your_db -f insert_script.sql

Now you can create SQL queries for accessing your data. E.g.:

    ST_Value(rast, 1, ST_SetSRID(ST_MakePoint(-100, 35), 4326))
    ST_Intersects(rast, ST_SetSRID(ST_MakePoint(-100, 35), 4326));

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

  • since the DEM is about 1Tb I would like to avoid pulling it inside the PostgreSQL database
    – Martin
    Commented Oct 24, 2019 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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