ENVI binary files to numpy arrays using GDAL

Most of my work currently revolves around raster processing of one kind or another. Being able to take advantage of the extensive libraries within Python (e.g. Numpy) is extremely helpful when carrying out such tasks. This blog really got me onto this.

Before trying to make use of the below functions, make sure you have access to the osgeo libraries - if you don’t, see here.

Below is a Python function that you can take that will get any given ENVI type binary file (so long as it has a .hdr file) into a numpy array.

 
from __future__ import division
import os
import sys

from osgeo import gdal, gdalconst
from osgeo.gdalconst import *


def ENVI_raster_binary_to_2d_array(file_name):
'''
Converts a binary file of ENVI type to a numpy array.
Lack of an ENVI .hdr file will cause this to crash.
'''


driver.Register()

inDs = gdal.Open(file_name, GA_ReadOnly)

if inDs is None:
print "Couldn't open this file: " + file_name
print '\nPerhaps you need an ENVI .hdr file?
sys.exit("Try again!")
else:
print "%s opened successfully" %file_name

print '~~~~~~~~~~~~~~'
print 'Get image size'
print '~~~~~~~~~~~~~~'
cols = inDs.RasterXSize
rows = inDs.RasterYSize
bands = inDs.RasterCount

print "columns: %i" %cols
print "rows: %i" %rows
print "bands: %i" %bands

print '~~~~~~~~~~~~~~'
print 'Get georeference information'
print '~~~~~~~~~~~~~~'
geotransform = inDs.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

print "origin x: %i" %originX
print "origin y: %i" %originY
print "width: %2.2f" %pixelWidth
print "height: %2.2f" %pixelHeight

# Set pixel offset.....
print '~~~~~~~~~~~~~~'
print 'Convert image to 2D array'
print '~~~~~~~~~~~~~~'
band = inDs.GetRasterBand(1)
image_array = band.ReadAsArray(0, 0, cols, rows)
image_array_name = file_name
print type(image_array)
print image_array.shape

return image_array, pixelWidth, (geotransform, inDs)


# The function can be called as follows:
# image_array, post, envidata = ENVI_raster_binary_to_2d_array(file_name)
#
# Notes:
# Notice the tuple (geotransform, inDs) - this contains all of your map information
# (xy tie points, postings and coordinate system information)
# pixelWidth is assumed to be the same as pixelHeight in the above example, therefore
# representing the surface posting - if this is not the case for your data then you
# must change the returns to suit

There is a lot more you can return from the above function as you can see - adjust to taste. To import files of different types, have a look here and alter the line:

 
driver = gdal.GetDriverByName('ENVI')

Once you have your numpy array and you’ve done some processing or whatever, it is then useful to be able to pull the array back out, exporting it back as an ENVI (or other) binary. A “post” input is required in case you’ve resampled the original imported binary through your processing. If you haven’t resampled just pass in the original post as returned by the first function.

 
def ENVI_raster_binary_from_2d_array(envidata, file_out, post, image_array):
util.check_output_dir(file_out)
original_geotransform, inDs = envidata

rows, cols = image_array.shape
bands = 1

# Creates a new raster data source
outDs = driver.Create(file_out, cols, rows, bands, gdal.GDT_Float32)

# Write metadata
originX = original_geotransform[0]
originY = original_geotransform[3]

outDs.SetGeoTransform([originX, post, 0.0, originY, 0.0, -post])
outDs.SetProjection(inDs.GetProjection())

#Write raster datasets
outBand = outDs.GetRasterBand(1)
outBand.WriteArray(image_array)

new_geotransform = outDs.GetGeoTransform()
new_projection = outDs.GetProjection()

print "Output binary saved: ", file_out

return new_geotransform,new_projection,file_out

Thanks to Martin for cleaning some of this up!

Written on November 17, 2014