Getting elevation readings at any lat,lon from shapefile with scattered irregular points

Question

I have a set of x,y coordinates, each of which I'm trying to get an elevation reading for.

The elevation data source, however, is not a regular grid, but consists of a shapefile containing scattered, irregular points, as this QGIS screenshot shows:

I'm used to using grid-based files such as SRTM or ASTER, which have a reading on a regular grid, and using a TIF reader to extract the locations. However, the above data is obviously quite different.

What is the best way to deal with this? I'm not a GIS expert - much more a coder.

Approaches that occur to me are:

• Interpolate the points into a regular grid. Would gdal_grid as noted at https://gis.stackexchange.com/a/24610/58752 be workable?

• Convert the points to triangles somehow and then interpolate the x,y from the triangle that it is in.

Crucially, I want to do this using command-line processing rather than point-and-click in GIS if possible, so that it is easily reproducible from the original data download.

Presumably the above screenshot indicates that the elevation field is called 'ALTURA'.

Any pointers would be useful.

(I'm aware that SRTM covers the area that this second data source does, but the data is a ground-based survey and avoids the problems of incorrect readings in SRTM due to building heights and trees, which are proving problematic.)

PS Is there a simple way to get the points out of the shapefile as a CSV? I can't seem to get this using either ogr2ogr or gdallocationinfo.

Answer 1

1. gdal_grid is the "best" choice. It lets you create a DEM/DSM GeoTIFF directly from your SHP. There is no need to first convert your SHP to CSV, and rightly this idea shouldn't had came up in the first place. And, there is also no need for programming.
2. gdal_grid offers many interpolation algorithms but (unfortunately) they confuse the casual users. However, for serious use cases - where and when we need to be explicit on the "right" interpolation algorithm to use - gdal_grid fits the need.
3. Many years back, I came across a helpful step-by-step tutorial in QGIS Tutorials and Tips - Interpolating Point Data. The tutorial uses QGIS' built-in raster interpolation tool. [Note: QGIS 2.18.16 has a bug that causes the resultant TIFF to be placed elsewhere on the map. This can be temporarily fixed by explicitly overriding QGIS' CRS assumption, but hopefully the QGIS team will fix this bug soon.]
4. GRASS also offers many interpolation algorithms. Do check out GRASS.
5. In gdal_grid, the TIN algorithm is called the linear algorithm. An important point to be aware when using gdal_grid is that -outsize specifies number of cells (across and down), while in other GDAL tools, similarly-named parameters specify pixel resolution.

6. Using the SHP sample (Arlington_Soundings_2007_stpl83.shp) in the aforementioned tutorial, the following command gives the expected output:-

   gdal_grid -ot float32 -outsize 2502 4096 -zfield elevation -a linear Arlington_Soundings_2007_stpl83.shp result.tif

7. Also, do be aware that the heights in your SRTM/ASTER data set may be in EGM96, while your SHP data set may in your country's vertical datum. Or, the "heights" in both data sets may be just geoid undulations.

Answer 2

I figured out that this is:

#!/bin/sh


# Usage:
#  Test, giving coarse resolution (50)
#    ./convert.sh 50
#  Standard, giving fine resolution (1024):
#    ./convert.sh


# Reproject the supplied shapefile as WGS84; will create data_wgs84.dbf.{dbf|prj|shp|shx}
ogr2ogr -f "ESRI Shapefile" data_wgs84.shp data.shp -s_srs data.prj -t_srs EPSG:4326

# Save the Shapefile as CSV
ogr2ogr -f CSV data_wgs84.csv -a_srs data_wgs84.prj -lco GEOMETRY=AS_XY data_wgs84.shp

# Create a .vrt file, which acts as the metadata definition for gdal_grid and which points to the CSV file
cat <<EOF > data_wgs84.vrt
<OGRVRTDataSource>
    <OGRVRTLayer name="data_wgs84">
        <SrcDataSource>data_wgs84.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="X" y="Y" z="ALTURA"/>
    </OGRVRTLayer>
</OGRVRTDataSource>
EOF

# Determine resolution; 50 is a good test
RESOLUTION=${1:-1024}

# Convert the points to a grid; see: http://www.gdal.org/gdal_grid.html
# Using "-outsize 1024 1024" will give better resolution than say "-outsize 256 256"
# To limit the tile location, add e.g. "-txe minLon maxLon -tye maxLat minLat"
gdal_grid -l data_wgs84 -outsize ${RESOLUTION} ${RESOLUTION} -a_srs WGS84 data_wgs84.vrt elevations.tiff

# View .tiff properties using:
gdalinfo elevations.tiff

# Clean up
rm data.*
rm data_wgs84.*

# Confirm
echo "The elevations.tiff file has been created. Testing gdallocationinfo :"

# Confirm useful elevation info using
gdallocationinfo -valonly -wgs84 elevations.tiff latValue lonValue