Since it isn’t always easy to get GDAL working with Conda, I decided to try reading the PRISM BIL file using built in Python modules and NumPy.

The first example I found (Plotting USGS BIL DEMs in Python) provided 95% of the steps needed to get this working using the struct module. I only needed to change the datatype to a float32/short. I decided to hardcode the width, height, datatype, number of bands, and nodata value for now, but all of these could/should be read from the HDR file.

import struct prism_path = 'PRISM_tmean_30yr_normal_4kmM2_annual_bil.bil' prism_cols = 1405 prism_rows = 621 prism_nodata = -9999 bil_f = open(bil_path, 'rb') bil_data = bil_f.read() bil_f.close() # Unpack binary data into a flat tuple z s = '<%df' % (int(prism_cols * prism_rows),) z = struct.unpack(s, bil_data) prism_array = np.zeros((prism_rows, prism_cols), dtype=np.float32) for r in range(0, prism_rows): for c in range(0, prism_cols): prism_array[r][c] = float(z[((prism_cols) * r) + c]) # prism_array.itemset((r, c), float(z[((prism_cols) * r) + c])) prism_array[prism_array == prism_nodata] = np.nan

I then realized that NumPy should be able to easily (and quickly) read a simple single band binary array file using fromfile().

import numpy as np prism_path = 'PRISM_tmean_30yr_normal_4kmM2_annual_bil.bil' prism_cols = 1405 prism_rows = 621 prism_nodata = -9999 prism_array = np.fromfile(bil_path, dtype=np.float32) prism_array = prism_array.reshape(prism_rows, prism_cols) prism_array[prism_array == prism_nodata] = np.nan

Might as well read the width, height, and nodata values from the HDF file, since it is avaiable. You could get the datatype also.

import numpy as np def read_prism_hdr(hdr_path): """Read an ESRI BIL HDR file""" with open(hdr_path, 'r') as input_f: header_list = input_f.readlines() return dict(item.strip().split() for item in header_list) def read_prism_bil(bil_path): """Read an array from ESRI BIL raster file""" hdr_dict = read_prism_hdr(bil_path.replace('.bil', '.hdr')) # For now, only use NROWS, NCOLS, and NODATA # Eventually use NBANDS, BYTEORDER, LAYOUT, PIXELTYPE, NBITS prism_array = np.fromfile(bil_path, dtype=np.float32) prism_array = prism_array.reshape( int(hdr_dict['NROWS']), int(hdr_dict['NCOLS'])) prism_array[prism_array == float(hdr_dict['NODATA'])] = np.nan return prism_array prism_path = 'PRISM_tmean_30yr_normal_4kmM2_annual_bil.bil' prism_array = read_prism_bil(prism_path)