import os
import numpy as np
from plio.io.io_gdal import GeoDataset
from libpyhat import spectral_data
from libpyhat.IO.cube_to_df import cube_to_df
[docs]
def read(f):
if os.path.splitext(f)[-1] == "hdr":
# GDAL wants the img, but many users aim at the .hdr
f = os.path.splitext(f)[:-1] + ".img"
geodata = GeoDataset(f)
wv = dict(
(k, v) for (k, v) in geodata.metadata.items() if "Band" in k
) # Only get the 'Band_###' keys
wavelengths = [
float(j.split(" ")[0])
for i, j in sorted(wv.items(), key=lambda x: int(x[0].split("_")[-1]))
]
# translate key names to numerical wvl values
# read each spectral band in the data cube
bands = []
for i in range(1, len(wavelengths) + 1, 1):
try:
bands.append(geodata.read_array(i))
except:
pass
cube = np.array(bands, dtype=float) # stack the bands together into a cube
cube = np.transpose(cube, (2, 1, 0))
df = cube_to_df(cube, wavelengths) # convert to a flattened data frame
data = spectral_data.SpectralData(df, name=None, geodata=geodata)
return data