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

Question

我有来自 PRISM Climate Group 的降水数据，现在以 .bil 格式提供（我认为是 ESRI BIL），我希望能够使用 Python 读取这些数据集。

我已经安装了 spectral 软件包，但

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.

spectrum 的文档清楚地表明它支持 BIL 文件，任何人都可以了解这里发生的情况吗？我也愿意使用 GDAL，据说它支持类似/等效的 ESRI EHdr 格式，但我找不到任何好的代码片段来开始。

Answer 1

现在已经是 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)

Answer 2

好吧，我很抱歉发布一个问题，然后自己这么快回答，但我发现了一组来自 犹他州立大学 的很好的课程幻灯片，其中有关于使用 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]

编辑 2014 年 5 月 27 日

我已经根据上面的答案进行了构建，并想在这里分享它，因为似乎缺乏文档。我现在有一个类，它有一个主要方法，该方法将 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 3

您已经找到了读取文件的良好解决方案，因此这个答案只是为了阐明您遇到的错误。

问题在于 spectral 软件包不支持 Esri 多波段栅格格式。 BIL（Band Interleaved by Line）不是一种特定的文件格式；相反，它是一种数据交错方案（如 BIP 和 BSQ），可用于多种文件格式。光谱包确实支持其识别的文件格式的 BIL（例如 ENVI、Erdas LAN），但 Esri 栅格标头不是其中之一。