I was wondering if anyone has some experience in getting elevation data from a raster without using ArcGIS, but rather get the information as a python list or dict?

I get my XY data as a list of tuples:

xy34 =[perp_obj[j].CalcPnts(float(i.dist), orientation) for j in range (len(perp_obj))]

I'd like to loop through the list or pass it to a function or class-method to get the corresponding elevation for the xy-pairs.

I did some research on the topic and the gdal API sounds promising. Can anyone advice me how to go about things, pitfalls, sample code?

GDAL is not an option as I can't edit the system path variable on the machine I'm working on!

Does anyone know about a different approach?

  • 2
    unfortunately, you really do need to get GDAL running on your system to do anything with a raster in Python. With "can't edit the system path variable on the machine", are you referring to these instructions? I find this installation method very poor, and I don't use it nor recommend it. If you are using Windows, install GDAL/Python the simple way.
    – Mike T
    Jul 18, 2012 at 22:12
  • Yes, I was indeed. I'm not at work right now but I'll check out the link you posted. Looks promising! Thanks for coming back to my question!
    – LarsVegas
    Jul 23, 2012 at 10:19
  • I've used the installer by Christoph Gohlke (linked above) on many work computers, and it is really simple. You only need to ensure that you match the version of Python and either 32- or 64-bit Windows. While you are at it, you should also get NumPy from the same place, since that is needed by GDAL, as shown in the answers below.
    – Mike T
    Jul 23, 2012 at 20:11

4 Answers 4


Here's a more programmatic way of using GDAL than @Aragon's answer. I've not tested it, but it is mostly boiler-plate code that has worked for me in the past. It relies on Numpy and GDAL bindings, but that's about it.

import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
from numpy import ma

def maFromGDAL(filename):
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    if dataset is None:
        raise Exception()

    # Get the georeferencing metadata.
    # We don't need to know the CRS unless we want to specify coordinates
    # in a different CRS.
    #projection = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    # We need to know the geographic bounds and resolution of our dataset.
    if geotransform is None:
        dataset = None
        raise Exception()

    # Get the first band.
    band = dataset.GetRasterBand(1)
    # We need to nodata value for our MaskedArray later.
    nodata = band.GetNoDataValue()
    # Load the entire dataset into one numpy array.
    image = band.ReadAsArray(0, 0, band.XSize, band.YSize)
    # Close the dataset.
    dataset = None

    # Create a numpy MaskedArray from our regular numpy array.
    # If we want to be really clever, we could subclass MaskedArray to hold
    # our georeference metadata as well.
    # see here: http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
    # For details.
    masked_image = ma.masked_values(image, nodata, copy=False)
    masked_image.fill_value = nodata

    return masked_image, geotransform

def pixelToMap(gt, pos):
    return (gt[0] + pos[0] * gt[1] + pos[1] * gt[2],
            gt[3] + pos[0] * gt[4] + pos[1] * gt[5])

# Reverses the operation of pixelToMap(), according to:
# https://en.wikipedia.org/wiki/World_file because GDAL's Affine GeoTransform
# uses the same values in the same order as an ESRI world file.
# See: http://www.gdal.org/gdal_datamodel.html
def mapToPixel(gt, pos):
    s = gt[0] * gt[4] - gt[3] * gt[1]
    x = (gt[4] * pos[0] - gt[1] * pos[1] + gt[1] * gt[5] - gt[4] * gt[2]) / s
    y = (-gt[3] * pos[0] + gt[0] * pos[1] + gt[3] * gt[2] - gt[0] * gt[5]) / s
    return (x, y)

def valueAtMapPos(image, gt, pos):
    pp = mapToPixel(gt, pos)
    x = int(pp[0])
    y = int(pp[1])

    if x < 0 or y < 0 or x >= image.shape[1] or y >= image.shape[0]:
        raise Exception()

    # Note how we reference the y column first. This is the way numpy arrays
    # work by default. But GDAL assumes x first.
    return image[y, x]

    image, geotransform = maFromGDAL('myimage.tif')
    val = valueAtMapPos(image, geotransform, (434323.0, 2984745.0))
    print val
    print('Something went wrong.')
  • 1
    see the edit to my question...thanks for posting anyway! I upvoted it.
    – LarsVegas
    Jul 16, 2012 at 13:19
  • 1
    Ah damn! Well at least it's here for posterity. TBH, the maths in mapToPixel() and pixelToMap() are the important bit, as long as you can create a numpy array (or a regular Python one, but they're generally not as efficient for this sort of thing), and get the array's geographical bounding box. Jul 16, 2012 at 13:28
  • 1
    +1 for the comment (and code) about reversing the parameters to the numpy array. I was searching everywhere for a bug in my code, and this swap fixed it!
    – aldo
    Feb 19, 2015 at 23:15
  • 1
    Then I suggest your matrix (gt in the example) is wrong. An affine matrix as used in CGAL (see: gdal.org/gdal_datamodel.html) is generally invertible (otherwise you have some funky scaling values going on). So where we have g = p.A we can also do p = g.A^-1 Numpy.linalg is a bit heavyweight for our purposes - we can reduce everything to two simple equations. Jan 10, 2017 at 17:44
  • 1
    I have re-edited the code to use simple algebra rather than numpy linalg. If the maths is wrong, fix the Wikipedia page. Jan 11, 2017 at 16:53

Check out my answer here... and read here for some information. The following info was taken from Geotips:

With gdallocationinfo, we can query the elevation at a point:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

The output of the above command has the form:

   Location: (87360P,19679L)
Band 1:
   Value: 1418

This means, that the elevation value at the provided geolocation is 1418.

  • Just found out I can't use GDAL as I'm not able to edit my system variable on the machine I'm working on. Thanks for the input though.
    – LarsVegas
    Jul 16, 2012 at 13:03

See e.g. this code which is based on GDAL (and Python, no numpy needed): https://github.com/geometalab/retrieve-height-service

  • It's unfortunate that the code doesn't seem to be open source licensed.
    – user10112
    Jun 19, 2019 at 22:48
  • Now it has :-).
    – Stefan
    Jun 21, 2019 at 5:31

The provided python code extracts the value data of a raster cell based on given x,y coords. It is a slightly alterd version of an example from this excellent source. It is based on GDAL and numpy which are not part of the standard python distribution. Thanks to @Mike Toews for pointing out the Unofficial Windows Binaries for Python Extension Packages to make installation and use quick and easy.

import os, sys, time, gdal
from gdalconst import *

# coordinates to get pixel values for
xValues = [122588.008]
yValues = [484475.146]

# set directory

# register all of the drivers
# open the image
ds = gdal.Open('i25gn1_131.img', GA_ReadOnly)

if ds is None:
    print 'Could not open image'

# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# loop through the coordinates
for xValue,yValue in zip(xValues,yValues):
    # get x,y
    x = xValue
    y = yValue

    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    # create a string to print out
    s = "%s %s %s %s " % (x, y, xOffset, yOffset)

    # loop through the bands
    for i in xrange(1,bands):
        band = ds.GetRasterBand(i) # 1-based index
        # read data and add the value to the string
        data = band.ReadAsArray(xOffset, yOffset, 1, 1)
        value = data[0,0]
        s = "%s%s " % (s, value) 
    # print out the data string
    print s
    # figure out how long the script took to run
  • It seems like this is just a less generic, less flexible version of what MerseyViking offered above?
    – WileyB
    Feb 17, 2016 at 14:08

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.