1

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:

enter image description here

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.

  • What flavour "coder" are you? Python/Java/R/Javascript? There's ways to access shapefile data from those... – Spacedman Dec 9 '15 at 11:12
  • PHP but can easily make sense of Python and call that via shell. But not really clear where to start as am not very familiar with actual GIS formats or the relevant technique to use here. Any pointers appreciated! – fooquency Dec 9 '15 at 11:16
  • ogr2ogr -f CSV pts.csv -lco GEOMETRY=AS_XY pts.shp gets you point features out of a shapefile as a CSV. As to the main question, some kind of spatial interpolation (nearest neighbour, IDW, or Kriging) is needed to interpolate altitude values at non-measured locations... – Spacedman Dec 9 '15 at 11:28
  • I put together this example for interpolation based on triangles, it might be start for you: github.com/mdsumner/xyztrisurf/blob/master/xyztrisurf.md Basically it samples from a grid to simulate points with elevation, then reconstructs the grid with different sample sizes - the method is the same as Matlab's default for "griddata", but otherwise not packaged in R in an easy form afaik. – mdsumner Dec 9 '15 at 13:51
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.

0

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
  • You don't need to convert your point shapefile to CSV. The gdal_grid utility can read any OGR supported datasource, i.e your WGS84 shapefile. From the documentation: gdal_grid [options] <src_datasource> <dst_filename> where "src_datasource: Any OGR supported readable datasource." – user2856 Jul 1 '17 at 10:19

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.