3

I asked the question below on StackOverflow, but haven't received any feedback there yet:

I have precipitation data from the PRISM Climate Group which are now offered in .bil format (ESRI BIL, I think) and I'd like to be able to read these datasets with Python.

I've installed the spectral package, but the open_image() method returns an error:

def ReadBilFile(bil):
    import spectral as sp
    b = sp.open_image(bil)
ReadBilFile(r'G:\truncated\ppt\1950\PRISM_ppt_stable_4kmM2_1950_bil.bil')

IOError: Unable to determine file type or type not supported.

The documentation for spectral clearly says that it supports BIL files, can anyone shed any light on what's happening here? I am also open to using GDAL, which supposedly supports the similar/equivalent ESRI EHdr format, but I can't find any good code snipets to get started.

  • 1
    Hi Jason. I see that you've got a solution based on GDAL, nice one! For ESRI BIL/BSQ files you can access without GDAL, just flat binary IO. Have a look at this post gis.stackexchange.com/questions/96050/… . I admit that binary IO in python isn't for the faint hearted but it might serve you in the future when you want to get/save individual pixels. – Michael Stimson May 24 '14 at 9:02
  • Thank you for the link, I'll keep that in mind for future reference. – Jason May 24 '14 at 12:31
  • 1
    An explanation of the reason for this error is here. – bogatron Jun 17 '14 at 16:31
8

Ok, I'm sorry to post a question and then answer it myself so quickly, but I found a nice set of course slides from Utah State University that has a lecture on opening raster image data with GDAL. For the record, here is the code I used to open the PRISM Climate Group datasets (which are in the EHdr format).

def ReadBilFile(bil):
    import gdal
    gdal.GetDriverByName('EHdr').Register()
    img = gdal.Open(bil)
    band = img.GetRasterBand(1)
    data = band.ReadAsArray()
    return data

if __name__ == '__main__':
    a = ReadBilFile(r'G:\truncated\ppt\1950\PRISM_ppt_stable_4kmM2_1950_bil.bil')
    print a[44, 565]

EDIT 5/27/2014

I've built upon my answer above and wanted to share it here since the documentation seems to be lacking. I now have a class with one main method that reads the BIL file as an array and returns some key attributes.

class BilFile(object):

    def __init__(self, bil_file):
        self.bil_file = bil_file
        self.hdr_file = bil_file.split('.')[0]+'.hdr'

    def get_array(self, mask=None):
        import gdal
        import gdalconst
        self.nodatavalue, self.data = None, None
        gdal.GetDriverByName('EHdr').Register()
        img = gdal.Open(self.bil_file, gdalconst.GA_ReadOnly)
        band = img.GetRasterBand(1)
        self.nodatavalue = band.GetNoDataValue()
        self.ncol = img.RasterXSize
        self.nrow = img.RasterYSize
        geotransform = img.GetGeoTransform()
        self.originX = geotransform[0]
        self.originY = geotransform[3]
        self.pixelWidth = geotransform[1]
        self.pixelHeight = geotransform[5]
        self.data = band.ReadAsArray()
        self.data = np.ma.masked_where(self.data==self.nodatavalue, self.data)
        if mask is not None:
            self.data = np.ma.masked_where(mask==True, self.data)
        return self.nodatavalue, self.data

I call this class using the following function where I use GDAL's vsizip function to read the BIL file directly from a zip file.

import prism
def getPrecipData(years=None):
    grid_pnts = prism.getGridPointsFromTxt()
    flrd_pnts = np.array(pd.read_csv(r'D:\truncated\PrismGridPointsFlrd.csv').grid_code)
    mask = prism.makeGridMask(grid_pnts, grid_codes=flrd_pnts)
    for year in years:
        bil = r'/vsizip/G:\truncated\PRISM_ppt_stable_4kmM2_{0}_all_bil.zip\PRISM_ppt_stable_4kmM2_{0}_bil.bil'.format(year)
        b = prism.BilFile(bil)
        nodatavalue, data = b.get_array(mask=mask)
        data *= mm_to_in
        b.write_to_csv(data, 'PrismPrecip_{}.txt'.format(year))
    return

# Get datasets
years = range(1950, 2011, 5)
getPrecipData(years=years)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.