Reading National Elevation Dataset (ArcGrid/GridFloat/IMG) with Python only tools?

I've found high precision elevation (1/3 - 1/9 arcsecond resolution) data from the National Elevation Dataset provided by the USGS. It comes int IMG, Arcgrid, and GridFloat formats. I know ArcGrid corresponds to a paid software package, but I'm trying to stick to using freely available tools.

I have GPS data that I'm trying to correlate with ground level.

Are there any python libraries which will let me transform gps data into altitude using a datafile in IMG, ArcGrid, or GridFloat formats?

Answer

Working with the IMG file directly in python is straightforward with the GDAL bindings. For example, you can read the data directly into a NumPy array:

from osgeo import gdal
geo = gdal.Open('imgn36w100_11.img')
arr = geo.ReadAsArray()
print repr(arr)
array([[ 744.31896973,  743.68762207,  743.1116333 , ...,  550.42498779,
         553.77813721,  556.18640137],
       [ 744.22955322,  743.66082764,  743.05273438, ...,  552.05706787,

         554.81365967,  557.55877686],
       [ 744.0133667 ,  743.49041748,  743.00061035, ...,  553.0123291 ,
         555.78076172,  558.01312256],
       ...,
       [ 568.70880127,  567.33666992,  566.56170654, ...,  447.68035889,
         447.68804932,  447.65426636],
       [ 568.01116943,  566.95739746,  564.23382568, ...,  447.6696167 ,
         447.71224976,  447.62734985],
       [ 565.62896729,  562.65325928,  560.78759766, ...,  447.67129517,
         447.67529297,  447.65179443]], dtype=float32)

For a more complete example of plotting IMG format data, see this script which generated the image below. For your transformation of GPS data into altitude, you'll have to sample the resulting NumPy array.