6

我有来自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 格式，但我找不到任何好的代码片段来开始。

4

3 回答 3

9

现在是 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)
于 2017-11-11T16:03:04.230 回答
6

好的，我很抱歉发布一个问题然后这么快就自己回答了，但是我从犹他州立大学找到了一组不错的课程幻灯片，其中有关于使用 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)
于 2014-05-23T15:20:48.230 回答
2

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

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

于 2014-06-17T13:57:48.330 回答