Generating elevation 'database' from NED IMG data

Question

I'm starting to work with the NED 1/3 and 1/9 arc-second data. I've been using gdal in Python and QGIS to do some basic conversions of the raw .IMG data from USGS. My goal is to get elevation data of the San Francisco Bay Area that I can look up with coordinates. It seemed like the most straightforward way to do this was to use gdal to convert to .XYZ format - this resulted in a really big file from a small piece of the area I'm using as a test.

Eventually I'd like to do this for the state of California, the USA, etc. so scalability is a concern. Although I envision some sort of database, if gleaning the info directly from the IMG or some other format is better I'm open to the idea.

A similar poster reads data from the files directly - would this be a good route to go from a memory/speed point of view? What are the tradeoffs? If this is a better method than converting to raw data directly, can anyone point me to some resources on how to accomplish this in python?

What is the best way to get to where I want to be?

Answer 1

assuming that the goal is to determine an elevation with a user-supplied set of coordinates, there should be a way to do this with gdal - i use a function like this, which came almost verbatim from this answer (assumes raster is not rotated).

rast = 'myIMGfile.img'
mx = XCOORDINATE
my = YCOORDINATE

def readRastPix(rast,mx,my):
    src_ds=gdal.Open(rast) 
    gt=src_ds.GetGeoTransform()
    rb=src_ds.GetRasterBand(1)
    gdal.UseExceptions() #so it doesn't print to screen everytime point is outside grid
    # Convert from map to pixel coordinates.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel
    try: #in case raster isnt full extent
        structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_Float32) #Assumes 32 bit int aka 'float'
        intval = struct.unpack('f' , structval)
        #had to add 0.0000... so that it wouldn't truncate to integer and fail with constraint error my database
        val=intval[0]
        if intval[0]<-9999:
            val=-9999
    except:
       pass
    src_ds=None
    return(val)

in my experience, using gdal directly on binary formats (not that familiar with IMG) is blazing fast. I would think much faster than accessing an ASCII version of the same data...