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

【问题标题】：Reading zipped ESRI BIL files with Python使用 Python 读取压缩的 ESRI BIL 文件
【发布时间】：2014-07-12 01:33:07
【问题描述】：

我有来自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.

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

【问题讨论】：

标签： python-2.7 image-processing zipfile gdal spectral

【解决方案1】：

好的，很抱歉发布一个问题然后自己这么快回答，但我从Utah State University 找到了一组不错的课程幻灯片，其中有关于使用 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 函数来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)

【讨论】：

【解决方案2】：

您已经找到了读取文件的好方法，所以这个答案只是为了阐明您遇到的错误。

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

【讨论】：

感谢您的解释，这很有道理。

【解决方案3】：

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

【讨论】：