I am utterly new to GIS, so please bear with me if I am providing incomplete information here.

I have a DEM file (.dt2). What I am looking to do is as follows:

  1. take a series of points specified as lat, long
  2. calculate their elevations using the DEM data
  3. Output these points in as a list of ECEF coords

I need to do this in either GRASS or QGIS (using Python scripting, not the GUI).

I don't need a complete working code. I am only looking to be pointed in the right direction.

  • 1
    look at Python Script for getting elevation difference between two points: how to extract the elevation values ​​of a DEM and treat these values in Python – gene Apr 1 '15 at 6:13
  • I was hoping for a much more high level way of doing this. Please post your response as an answer. I'll hold out for a bit longer for an answer that uses QGIS API or GRASS. If nothing of that nature materializes I will accept your answer. – balajeerc Apr 1 '15 at 11:34
  • 1
    Perhaps gdallocationinfo gdal.org/gdallocationinfo.html executed with Python could be useful for you. – user30184 Apr 1 '15 at 12:53

If you want to do that with PyQGIS, the def Val_raster(x,y,layer,bands,gt) of Python Script for getting elevation difference between two points becomes:

def Val_raster(point,raster):
  return raster.dataProvider().identify(point,QgsRaster.IdentifyFormatValue).results().values()

for point in points.getFeatures():
       geompt = point.geometry().asPoint()
       print print Val_raster(geompt,DEM)

with a multiband raster (R,G,B values) the results are, for example

[203.0, 177.0, 202.0]
[194.0, 181.0, 199.0]
[109.0, 85.0, 101.0]

Or you can use GRASS command v.drape :

In the QGIS console (processing):


In GRASS GIS with grass.script (GRASS 6.4.x, GRASS 7.x ) or Pygrass (GRASS 7.x)

  • That QGIS snippet was exactly what I was looking for. Thanks! – balajeerc Apr 2 '15 at 7:22

I would point out, in addition, that the GRASS model, v.what.rast takes a point layer and raster layer as input and uploads the raster values to a given attrib column for all points. To run this in the python console (within a GRASS session) you would do something like:

import os
import grass.script as grass

input_points = "<your point list>"
input_dem = "<your elevation raster>"
# Prepare names for the GRASS vector and raster
pts, ext = os.path.splitext(input_points)
dem, ext = os.path.splitext(input_dem)

# Set the x,y values and separator as appropriate in the input list of locations
grass.run_command('v.in.ascii', input=input_points, output=pts, x=1,y=2, separator=comma)
grass.run_command('r.in.gdal' input=input_dem, output=dem)
grass.run_command('v.db.addtable', map=pts, columns="elev double precision")
grass.run_command('v.what.rast', map=pts, raster=dem, column='elev')

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.