Importing rasters into Python

I have been prompted to post a little piece on importing rasters into Python using GDAL.

The following code is something I use personally and may not be suited to everyones needs but it should be easy to adapt it.

def rasterImport(file):

    '''Import a raster file and return grid plus meta data. Returns three objects: data, meta, metadata '''

    # Tell the user that the function is working

    print '\nImporting ',file

    # Register all of the GDAL drivers

    gdal.AllRegister()

    # Open the file

    raster = gdal.Open(file)

    if raster is None:

        print 'Could not open ',file,'\n'

        sys.exit(1)

    # Get coordinate system and affine transformation parameters

    projec = raster.GetProjection()

    srs = osr.SpatialReference(wkt=projec)

    transf = raster.GetGeoTransform()

    # Get meta data

    # Get short name and then long name of driver

    driveshrt = raster.GetDriver().ShortName

    driver = gdal.GetDriverByName(driveshrt)

    metadata = driver.GetMetadata()

    metakeys = metadata.keys()

    # Read the file band to a 'matrix' called band

    band = raster.GetRasterBand(1)

    # Access data in rastern band as array

    data = band.ReadAsArray(0,0,cols,rows)
    # gdal_grid creates "upside down" files, hence this
    Yres=yres
    if Yres/np.fabs(Yres)==-1:
        Yres=-1*Yres
        data = np.flipud(data)

    # From here you can 'return' the data to the main routine

    # and skip the rest of this.
    #

    # Get the upper left coordinates and the pixel resolution

    ul_x = transf[0]

    ul_y = transf[3]

    xres = transf[1]

    yres = transf[5]

    # Get image size

    rows = raster.RasterYSize

    cols = raster.RasterXSize

    dims = {'xres':xres,'yres':yres,'rows':rows,'cols':cols}

    # Calculate all corners

    ll_x = ul_x

    ll_y = ul_y + (rows * yres)

    ur_x = ul_x + (cols * xres)

    ur_y = ul_y

    lr_x = ur_x

    lr_y = ll_y

    # Make dictionary of corners

    corners = {'ll':(ll_x,ll_y),'ul':(ul_x,ul_y),'ur':(ur_x,ur_y),'lr':(lr_x,lr_y)}

    # Get nodata value

    nandat = band.GetNoDataValue()

    # Get minimum and maximum value

    mindat = band.GetMinimum()

    maxdat = band.GetMaximum()

    # If the stats aren't available yet calculate the min-max

    if mindat is None or maxdat is None:

            (mindat,maxdat) = band.ComputeRasterMinMax(1)

    dats = [nandat,mindat,maxdat]

    # This next bit sums the values for a volume calculation

    sumvals = data[np.where(np.logical_not(data == nandat))]

    sumval = sumvals.sum()

    numvals = len(sumvals)

    avval = sumval/numvals

    volvals = sumval*(math.fabs((transf[1]*transf[5])))

    # Organise the metadata for export

    meta = {'transform':transf,'corners':corners,'dimension':dims,'dats':dats,'sumval':sumval,'numvals':numvals,'avval':avval,'volvals':volvals}

    # Tell user about the file

    if srs.IsProjected:

        print 'projcs 1: ',srs.GetAttrValue('projcs')

        print 'projcs 2: ',srs.GetAuthorityCode('projcs')

        print 'projcs 3: ',srs.GetAuthorityName('projcs')

    print 'geogcs 1:: ',srs.GetAttrValue('geogcs')

    print 'Sum of pixel values = ',sumval,' in range ',mindat,' to ',maxdat

    print 'From ',numvals,' of ',rows,' X ',cols,' = ',rows*cols,' pixels, or ',((1.0*numvals)/(rows*cols))*100,'%'

    print 'Average then is ',avval

    print 'Pixels size, x:',transf[1],' y:',transf[5]

    print '"Volume" is then the sum * 1 pixel area: ',volvals

    # Return the data matrix plus metadata to main routine

    return data, meta, metadata

That's it. Not the greatest, most efficient code but it does what I want.