Vjylku

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.

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