Sample raster cell value at coordinates?

I have many rasters in ESRI ASCII format (http://resources.arcgis.com/en/help/main/10.1/index.html#/Esri_ASCII_raster_format/009t0000000z000000/).

I need to extract cell values at given locations / coordinates. Can anyone suggest a python package to achieve this? I suspect there may be something in gdal tools but I have been unable to find anything so far.

I am looking for similar functionality to GMT grdtrack, with which you can pass a table of coordinates and retrieve the cell value. http://gmt.soest.hawaii.edu/doc/5.1.0/grdtrack.html

However I was hoping / wondering if this is possible in python as my previous and later stages of my analysis are all in python.

Solution

Of course, you can do this with GDAL. Provided that you have the coordinate where you want to sample your raster in the same projection you can do something like this:

from osgeo import gdal

def world2Pixel(gt, x, y):
  ulX = gt[0]
  ulY = gt[3]
  xDist = gt[1]
  yDist = gt[5]
  rtnX = gt[2]
  rtnY = gt[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / yDist)
  return (pixel, line)

dataset = gdal.Open(filename)
gt = dataset.GetGeoTransform()

pixel, line = world2Pixel(gt, x, y)

band = dataset.GetRasterBand(1)
value = band.ReadAsArray(pixel, line, 1, 1)[0, 0]