使用Python读取压缩的ESRI BIL文件

Question

我有来自PRISM气候组的降水数据,现在以.bil格式提供(ESRI BIL,我认为),我希望能够用Python读取这些数据集.

我已安装了光谱包,但该open_image()方法返回错误:

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.

光谱文档明确表示它支持BIL文件,任何人都可以了解这里发生的事情吗？我也愿意使用GDAL,据说它支持类似/等效的ESRI EHdr格式,但我找不到任何好的代码snipets来开始.

Answer 1

好的,我很抱歉发布一个问题然后自己快速回答,但是我找到了一套很好的犹他州立大学课程幻灯片,其中有一个关于用GDAL打开光栅图像数据的讲座.为了记录,这里是我用来打开PRISM气候组数据集的代码(它们是EHdr格式).

import gdal

def ReadBilFile(bil):

    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]

编辑5/27/2014

我已经在上面的答案的基础上,并希望在这里分享它,因为文档似乎缺乏.我现在有一个带有一个main方法的类,它将BIL文件作为数组读取并返回一些关键属性.

import gdal
import gdalconst

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):
        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

我使用以下函数调用此类,其中我使用GDAL的vsizip函数直接从zip文件读取BIL文件.

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)

Answer 2

它现在是2017年,有一个稍微好一点的选择.包rasterio支持bil文件.

>>>import rasterio
>>>tmean = rasterio.open('PRISM_tmean_stable_4kmD1_20060101_bil.bil')
>>>tmean.affine
Affine(0.041666666667, 0.0, -125.0208333333335,
       0.0, -0.041666666667, 49.9375000000025)
>>> tmean.crs
CRS({'init': 'epsg:4269'})
>>> tmean.width
1405
>>> tmean.height
621
>>> tmean.read().shape
(1, 621, 1405)