Getting elevation at lat/long from raster using python?

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?

Answer

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
os.chdir(r'D:\\temp\\AHN2_060')

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


if ds is None:
    print 'Could not open image'
    sys.exit(1)

# 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