# we use NumPy to do array processing of the cube data
import numpy as np

# we use SPy to read/write ENVI image files
import spectral.io.envi as envi

# the number of spectral cubes are to be processed
imgCount = 840

# read the first cube to get the spectral dimension
img0 = envi.open("cube_000.img.hdr")
bandCount = img0.shape[2]

# create a new output ENVI image containing the spectral waterfall
hdrDict = {
  "samples": bandCount,
  "lines": imgCount,
  "bands": 1,
  "data type": 5
}
outImage = envi.create_image("waterfall.hdr",hdrDict,force=True)
outData = outImage.open_memmap(writable=True)
print(outData.shape)

# start processing cubes
for iIdx in range(0,imgCount):

  # load the spectral data cube
  imgFilename = "cube_" + f"{iIdx:03d}" + ".img"
  print(imgFilename)
  img = envi.open(imgFilename+".hdr", imgFilename)
  imgData = img.load()

  # average the spatial dimension to make a spectral curve
  imgAvg = np.average(imgData,axis=(0,1))

  # copy this curve into the waterfall image
  for bIdx in np.arange(0,bandCount):
    outData[iIdx,bIdx] = imgAvg[bIdx]
