2

I try for too ling to read the envi file in python. Normally this is very easy. First, I extract necessary values from the .hdr:

for line in open("../../local_storage_only/aviris/ang20180626t145910_rfl_v2q2/ang20180626t145910_corr_v2q2_img.hdr"):

  if "data type" in line:
    _data_type = int(line.split("= ")[-1])
  if "interleave" in line:
    _interleave = line.split("= ")[-1]
  if "byte order " in line:
    _byte_order = int(line.split("= ")[-1])
  if "wavelength units" in line:
    _wavelength_units = line.split("= ")[-1]
  if "correction factors" in line:
    _correction_factors = line.split(" { ")[-1][:-2]
    _correction_factors = [float(element) for element in _correction_factors.split(" , ")]
  if "wavelength = " in line:
    _wavelength = line.split(" { ")[-1][:-2]
    _wavelength = [float(element) for element in _wavelength.split(" , ")]
    _wavelength_dict = dict(zip( _wavelength, range(len( _wavelength))))
  if "bands = " in line:
    _bands = int(line.split(" = ")[-1])
  if "lines =" in line:
    _rows = int(line.split(" = ")[-1])
  if "samples =" in line:
    _columns = int(line.split(" = ")[-1])
  if "data ignore value " in line:
    _ignore_value = int(line.split(" = ")[-1])
print(line)

Secondly, I open the file and convert it in to the right format:

tmp_file = open("../../local_storage_only/aviris/ang20180626t145910_rfl_v2q2/ang20180626t145910_corr_v2q2_img",
               encoding="utf-8")
tmp_img = np.fromfile(tmp_file,
                      np.dtype(np.float32).char)

tmp_img[tmp_img == _ignore_value] = np.nan

However, if I plot the results it looks like this:

enter image description here

I also tried to change the order of height, width and bands. I altered the sorting mechanism of np.reshape("A", "C", "F"). Does someone know what I do wrong?

3
+100

You can use numpy's reshape and transpose functions to reconstruct the desired result. And the dimensions of the "desired result" is used in one of two forms which is often up to the user to decide:

  • brc[band, row, col], e.g. to index band 1 use brc[0]
  • rcb[row, col, band], e.g. to index band 1 use rcb[:, :, 0]

Esri has a good example of BIL, BIP and BSQ to follow:

bil_src = """\
0 0 0 0 0 0 0 0 0 0 64 64 128 128 255 255 255 255 255 255 255 255 255 255
0 0 0 0 0 0 0 0 0 0 64 64 128 128 255 255 255 255 255 255 255 255 255 255
64 64 64 64 64 64 64 64 0 0 64 64 128 128 255 255 128 128 128 128 128 128 128 128
64 64 64 64 64 64 64 64 0 0 64 64 128 128 255 255 128 128 128 128 128 128 128 128
128 128 128 128 128 128 128 128 0 0 64 64 128 128 255 255 64 64 64 64 64 64 64 64
128 128 128 128 128 128 128 128 0 0 64 64 128 128 255 255 64 64 64 64 64 64 64 64
255 255 255 255 255 255 255 255 0 0 64 64 128 128 255 255 0 0 0 0 0 0 0 0
255 255 255 255 255 255 255 255 0 0 64 64 128 128 255 255 0 0 0 0 0 0 0 0
"""
bip_src = """\
0 0 255 0 0 255 0 64 255 0 64 255 0 128 255 0 128 255 0 255 255 0 255 255
0 0 255 0 0 255 0 64 255 0 64 255 0 128 255 0 128 255 0 255 255 0 255 255
64 0 128 64 0 128 64 64 128 64 64 128 64 128 128 64 128 128 64 255 128 64 288 128
64 0 128 64 0 128 64 64 128 64 64 128 64 128 128 64 128 128 64 255 128 64 288 128
128 0 64 128 0 64 128 64 64 128 64 64 128 128 64 128 128 64 128 255 64 128 255 64
128 0 64 128 0 64 128 64 64 128 64 64 128 128 64 128 128 64 128 255 64 128 255 64
255 0 0 255 0 0 255 64 0 255 64 0 255 128 0 255 128 0 255 255 0 255 255 0
255 0 0 255 0 0 255 64 0 255 64 0 255 128 0 255 128 0 255 255 0 255 255 0
""".replace('288', '255')  # Esri made a boo-boo
bsq_src = """\
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
64 64 64 64 64 64 64 64
64 64 64 64 64 64 64 64
128 128 128 128 128 128 128 128
128 128 128 128 128 128 128 128
255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
0 0 64 64 128 128 255 255
255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255
128 128 128 128 128 128 128 128
128 128 128 128 128 128 128 128
64 64 64 64 64 64 64 64
64 64 64 64 64 64 64 64
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
"""

These three examples describe the same data in three ways:

Esri example

Read the data as flat arrays:

import numpy as np

# Read text to flat arrays
bil = np.fromstring(bil_src, np.uint8, sep=' ')
bip = np.fromstring(bip_src, np.uint8, sep=' ')
bsq = np.fromstring(bsq_src , np.uint8, sep=' ')

Now reshape into the brc[band, row, col] form:

nrow = 8
ncol = 8
nband = 3  # i.e., RGB

bil_shape = (nrow, nband, ncol)
bip_shape = (nrow, ncol, nband)
bsq_shape = (nband, nrow, ncol)

bil_brc = bil.reshape(bil_shape).transpose((1, 0, 2))
bip_brc = bip.reshape(bip_shape).transpose((2, 0, 1))
bsq_brc = bsq.reshape(bsq_shape)

# these three arrays are identical
np.testing.assert_equal(bil_brc, bip_brc)
np.testing.assert_equal(bil_brc, bsq_brc)

Or into rcb[row, col, band] form, the transpose is slightly different:

bil_rcb = bil.reshape(bil_shape).transpose((0, 2, 1))
bip_rcb = bip.reshape(bip_shape)
bsq_rcb = bsq.reshape(bsq_shape).transpose((1, 2, 0))

# these three arrays are identical
np.testing.assert_equal(bil_rcb, bsq_rcb)
np.testing.assert_equal(bil_rcb, bip_rcb)

And here's what it looks like using matplotlib's imshow, which expects the last dimension to have RGB:

enter image description here


So to directly answer the question, with a few improvements to use a masked array instead of NaNs:

raw_fname = "../../local_storage_only/aviris/ang20180626t145910_rfl_v2q2/ang20180626t145910_corr_v2q2_img"
bil = np.fromfile(raw_fname, np.float32)

mask = bil == _ignore_value
if mask.any():
    bil = np.ma.array(bil, mask=mask)

bil_shape = (_rows, _bands, _columns)

# take your pick
brc = bil.reshape(bil_shape).transpose((1, 0, 2))
rcb = bil.reshape(bil_shape).transpose((0, 2, 1))
2

Ok, after some more research I found the answer: Unfortunately, numpy has no reshape order that correspond to BIL (Band-interleaved-by-line). However, you can resort it by yourself with numpy. This is the first draft from my jupyter notebook, but you can do it more optimized:

band_count = 0

for row_count in range(_rows):
    if bands == 0:
        arr = np.array(a[0: _columns])
        print(arr.shape)
    else:
        print(a[band_count * _bands *_columns : band_count * _bands * _columns + _columns].shape)
        arr = np.dstack((arr,
                       a[band_count * _bands *_columns : bands * _bands * _columns + _columns]))
        print(arr.shape)
    band_count += 1

if you want to use spectralpy, you can also do it and it will read produce the same output.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.