V0.1.5

.gitignore

@@ -6,4 +6,5 @@ dist/
 .DS_STORE
 build/
 myenv
-*.so*
+*.so*
+*mlr.c*

README.md

@@ -8,6 +8,8 @@
 
 `hyperquest`: A Python package for estimating image-wide quality estimation metrics of hyperspectral imaging (imaging spectroscopy). Computations are sped up and scale with number of cpus.
 
+Important: this package assumes your hyperspectral data is in ENVI format with a .HDR file.
+
 
 ## Installation Instructions
 
@@ -20,55 +22,26 @@ pip install hyperquest
 
 ## All Methods
 
-| **Result**          | **Method**                 | **Description**                                                                                                |
-|---------------------|----------------------------|----------------------------------------------------------------------------------------------------------------|
-| **SNR**             | `hrdsdc()`                 | Homogeneous regions division and spectral de-correlation (Gao et al., 2008)                                    |
-|                     | `rlsd()`                   | Residual-scaled local standard deviation (Gao et al., 2007)                                                    |
-|                     | `ssdc()`                   | Spectral and spatial de-correlation (Roger & Arnold, 1996)                                                     |
-| **Co-Registration** | `sub_pixel_shift()`        | Computes sub pixel co-registration between the VNIR & VSWIR imagers using skimage phase_cross_correlation      |
-| **Smile**           | `smile_metric()`           | Similar to MATLAB "smileMetric". Computes derivatives of O2 and CO2 absorption features.                       |
-
-
-
-## Usage example
-
-- see [Example Using EMIT](tutorials/example_using_EMIT.ipynb) for a recent use case.
+| **Category**             | **Method**                 | **Description**                                                                                                                               |
+|--------------------------|----------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------|
+| **SNR**                  | `hrdsdc()`                 | Homogeneous regions division and spectral de-correlation (Gao et al., 2008)                                                                   |
+|                          | `rlsd()`                   | Residual-scaled local standard deviation (Gao et al., 2007)                                                                                   |
+|                          | `ssdc()`                   | Spectral and spatial de-correlation (Roger & Arnold, 1996)                                                                                    |
+| **Co-Registration**      | `sub_pixel_shift()`        | Computes sub pixel co-registration between the VNIR & VSWIR imagers using skimage phase_cross_correlation                                     |
+| **Smile**                | `smile_metric()`           | Similar to MATLAB "smileMetric". Computes derivatives of O2 and CO2 absorption features (Dadon et al., 2010).                                 |
+|                          | `nodd_o2a()`               | Similar to method in Felde et al. (2003) to solve for nm shift at O2-A. Requires radiative transfer model run.                                |
+| **Radiative Transfer**   | `run_libradtran()`         | Runs libRadtran based on user input geometry and atmosphere at 1.0 cm-1. Saves to a .csv file for use in methods requiring radiative transfer.|
 
-```python
-import hyperquest
-import matplotlib.pyplot as plt
 
 
-# Define path to envi image header file
-envi_hdr_path = '/path/my_spectral_image.hdr'
-
-# get wavelengths
-wavelengths = hyperquest.read_center_wavelengths(envi_hdr_path)
-
-# compute SNR using HRDSDC method
-snr = hyperquest.hrdsdc(envi_hdr_path, n_segments=10000, 
-                        compactness=0.1, n_pca=3, ncpus=3)
-
-plt.scatter(wavelengths, snr, color='black', s=100, alpha=0.7)
-```
-![SNR Plot](tests/plots/demo_snr.png)
 
 
+## Usage example
 
+- see [SNR example](tutorials/example_using_EMIT.ipynb) where different SNR methods are computed over Libya-4.
+- see [Smile example ](tutorials/testing_smile_methods.ipynb) where different smile methods are computed over Libya-4.
 
-## Citing HyperQuest (STILL WORKING ON THIS, TODO:)
 
-If you use HyperQuest in your research, please use the following BibTeX entry.
-
-```bibtex
-@article{wilder202x,
-  title={x},
-  author={Brenton A. Wilder},
-  journal={x},
-  url={x},
-  year={x}
-}
-```
 
 
 ## References:
@@ -79,6 +52,8 @@ If you use HyperQuest in your research, please use the following BibTeX entry.
 
 - Dadon, A., Ben-Dor, E., & Karnieli, A. (2010). Use of derivative calculations and minimum noise fraction transform for detecting and correcting the spectral curvature effect (smile) in Hyperion images. IEEE Transactions on Geoscience and Remote Sensing, 48(6), 2603-2612.
 
+- Felde, G. W., Anderson, G. P., Cooley, T. W., Matthew, M. W., Berk, A., & Lee, J. (2003, July). Analysis of Hyperion data with the FLAASH atmospheric correction algorithm. In IGARSS 2003. 2003 IEEE International Geoscience and Remote Sensing Symposium. Proceedings (IEEE Cat. No. 03CH37477) (Vol. 1, pp. 90-92). IEEE.
+
 - Gao, L., Wen, J., & Ran, Q. (2007, November). Residual-scaled local standard deviations method for estimating noise in hyperspectral images. In Mippr 2007: Multispectral Image Processing (Vol. 6787, pp. 290-298). SPIE.
 
 - Gao, L. R., Zhang, B., Zhang, X., Zhang, W. J., & Tong, Q. X. (2008). A new operational method for estimating noise in hyperspectral images. IEEE Geoscience and remote sensing letters, 5(1), 83-87.
@@ -87,4 +62,5 @@ If you use HyperQuest in your research, please use the following BibTeX entry.
 
 - Scheffler, D., Hollstein, A., Diedrich, H., Segl, K., & Hostert, P. (2017). AROSICS: An automated and robust open-source image co-registration software for multi-sensor satellite data. Remote sensing, 9(7), 676.
 
-- Tian, W., Zhao, Q., Kan, Z., Long, X., Liu, H., & Cheng, J. (2022). A new method for estimating signal-to-noise ratio in UAV hyperspectral images based on pure pixel extraction. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 16, 399-408.
+- Thompson, D. R., Green, R. O., Bradley, C., Brodrick, P. G., Mahowald, N., Dor, E. B., ... & Zandbergen, S. (2024). On-orbit calibration and performance of the EMIT imaging spectrometer. Remote Sensing of Environment, 303, 113986.
+

hyperquest/libradtran.py

@@ -0,0 +1,222 @@
+# Import libraries
+import os
+import numpy as np
+import pandas as pd
+from spectral import *
+
+
+
+def get_libradtran_install_path(user_input_path):
+    '''
+    TODO
+    '''
+    user_input_path = os.path.abspath(user_input_path)
+
+    # Check if they already gave the bin directory
+    if os.path.basename(user_input_path) == 'bin' and os.path.isdir(user_input_path):
+        return user_input_path
+
+    # Check if 'bin' exists inside the given directory
+    bin_path = os.path.join(user_input_path, 'bin')
+    if os.path.isdir(bin_path):
+        return bin_path
+    else:
+        raise FileNotFoundError('Could not determine location of libRadtran installation.')
+
+
+    return
+
+
+
+def get_libradtran_output_dir(hdr_path):
+    '''
+    TODO
+    '''
+
+    filename = os.path.splitext(os.path.basename(hdr_path))[0]
+
+    # directory where hdr_path is located
+    parent_dir = os.path.dirname(os.path.abspath(hdr_path))
+
+    # Define the rtm output directory
+    lrt_out_dir = os.path.join(parent_dir, f'rtm-{filename}')
+
+    # Create the directory if it does not exist
+    os.makedirs(lrt_out_dir, exist_ok=True)
+
+    return lrt_out_dir
+
+
+
+
+def write_lrt_inp(o3, h , aod, a, out_str, umu, phi0, phi, sza, lat_inp, lon_inp, doy, altitude_km,
+              atmos, path_to_libradtran_bin, lrt_dir, path_to_libradtran_base):
+    '''
+
+    adapted from: https://github.com/MarcYin/libradtran
+
+        Cm-1 resolution of libRadtran ... which is output from lrt as 0.1 nm spectral resolution.
+
+
+    '''
+    foutstr = out_str[0] + out_str[1]
+    fname = f'{lrt_dir}/lrt_h20_{h}_aot_{aod}_alb_{a}_alt_{round(altitude_km*1000)}_{foutstr}'
+    with open(f'{fname}.INP', 'w') as f:
+        f.write(f'source solar {path_to_libradtran_base}/data/solar_flux/kurudz_0.1nm.dat\n') 
+        f.write('wavelength 549 860\n')  
+        f.write(f'atmosphere_file {path_to_libradtran_base}/data/atmmod/afgl{atmos}.dat\n')
+        f.write(f'albedo {a}\n') 
+        f.write(f'umu {umu}\n') # Cosine of the view zenith angle
+        f.write(f'phi0 {phi0}\n') # SAA
+        f.write(f'phi {phi}\n') # VAA
+        f.write(f'sza {sza}\n')  # solar zenith angle
+        f.write('rte_solver disort\n')
+        f.write(f'number_of_streams 8\n')
+        f.write('pseudospherical\n')
+        f.write(f'latitude {lat_inp}\n')
+        f.write(f'longitude {lon_inp}\n')
+        f.write(f'day_of_year {doy}\n') 
+        f.write(f'mol_modify O3 {o3} DU\n')   
+        f.write(f'mol_abs_param reptran fine\n')   # Fine cm-1
+        f.write(f'mol_modify H2O {h} MM\n')   
+        f.write(f'crs_model rayleigh bodhaine \n')  
+        f.write(f'zout {out_str[0]}\n')  
+        f.write(f'altitude {altitude_km}\n')    
+        f.write(f'aerosol_default\n')  
+        f.write(f'aerosol_species_file continental_average\n') 
+        f.write(f'aerosol_set_tau_at_wvl 550 {aod}\n')  
+        f.write(f'output_quantity transmittance\n')  
+        f.write(f'output_user lambda {out_str[1]}\n')  
+        f.write('quiet')
+    cmd = f'{path_to_libradtran_bin}/uvspec < {fname}.INP > {fname}.out'
+    return cmd
+
+
+
+
+def write_lrt_inp_irrad(o3, h , aod, a, out_str, umu, phi0, phi, sza, lat_inp, lon_inp, doy, altitude_km,
+                        atmos, path_to_libradtran_bin, lrt_dir, path_to_libradtran_base):
+    # Run here manually for irrad
+    fname = f'{lrt_dir}/lrt_h20_{h}_aot_{aod}_alt_{round(altitude_km*1000)}_IRRAD'
+    with open(f'{fname}.INP', 'w') as f:
+        f.write(f'source solar {path_to_libradtran_base}/data/solar_flux/kurudz_0.1nm.dat\n') 
+        f.write('wavelength 549 860\n')  
+        f.write(f'atmosphere_file {path_to_libradtran_base}/data/atmmod/afgl{atmos}.dat\n')
+        f.write(f'albedo {a}\n')  
+        f.write(f'sza {sza}\n')  
+        f.write('rte_solver disort\n')  
+        f.write(f'number_of_streams 8\n')
+        f.write('pseudospherical\n')
+        f.write(f'latitude {lat_inp}\n')
+        f.write(f'longitude {lon_inp}\n')
+        f.write(f'day_of_year {doy}\n')  
+        f.write(f'zout {altitude_km}\n')  
+        f.write(f'aerosol_default\n')  
+        f.write(f'aerosol_species_file continental_average\n')  
+        f.write(f'aerosol_set_tau_at_wvl 550 {aod}\n')  
+        f.write(f'mol_modify O3 {o3} DU\n')  
+        f.write(f'mol_abs_param reptran fine\n')   # Fine cm-1
+        f.write(f'mol_modify H2O {h} MM\n')    
+        f.write(f'crs_model rayleigh bodhaine \n')  
+        f.write(f'output_user lambda edir edn \n') 
+        f.write('quiet')
+    cmd = f'{path_to_libradtran_bin}/uvspec < {fname}.INP > {fname}.out'
+    return cmd
+
+
+
+def lrt_create_args_for_pool(h,
+                             aod,
+                             altitude_km,
+                             umu, phi0, 
+                             phi,vza,
+                             sza, lat_inp,
+                             lon_inp, doy, atmos, 
+                             o3, albedo,
+                             lrt_dir, path_to_libradtran_bin):
+    '''
+    TODO
+    '''
+    # Run the LRT LUT pipeline
+    path_to_libradtran_base = os.path.dirname(path_to_libradtran_bin)
+
+    lrt_inp = []
+    lrt_inp_irrad = []
+
+    # path radiance run
+    cmd = write_lrt_inp(o3, h,aod,0, ['toa','uu'], umu, phi0, phi, sza, 
+                    lat_inp, lon_inp, doy, altitude_km, atmos, path_to_libradtran_bin, 
+                    lrt_dir, path_to_libradtran_base)
+    lrt_inp.append([cmd,path_to_libradtran_bin])
+
+    # upward transmittance run
+    cmd = write_lrt_inp(o3, h,aod,0, ['sur','eglo'], umu, phi0, phi, vza, 
+                    lat_inp, lon_inp, doy, altitude_km, atmos, path_to_libradtran_bin, 
+                    lrt_dir, path_to_libradtran_base)
+    lrt_inp.append([cmd,path_to_libradtran_bin])
+
+    # spherical albedo run 1
+    cmd = write_lrt_inp(o3, h,aod,0.15, ['sur','eglo'], umu, phi0, phi, sza, 
+                    lat_inp, lon_inp, doy, altitude_km, atmos, path_to_libradtran_bin, 
+                    lrt_dir, path_to_libradtran_base)
+    lrt_inp.append([cmd,path_to_libradtran_bin])
+
+    # spherical albedo run 2
+    cmd = write_lrt_inp(o3, h,aod,0.5, ['sur','eglo'], umu, phi0, phi, sza, 
+                    lat_inp, lon_inp, doy, altitude_km, atmos, path_to_libradtran_bin, 
+                    lrt_dir, path_to_libradtran_base)   
+    lrt_inp.append([cmd,path_to_libradtran_bin])
+
+    # incoming solar irradiance run
+    cmd = write_lrt_inp_irrad(o3, h,aod, albedo, ['toa','uu'], umu, phi0, phi, sza, 
+                    lat_inp, lon_inp, doy, altitude_km, atmos, path_to_libradtran_bin, 
+                    lrt_dir, path_to_libradtran_base)
+    lrt_inp_irrad.append([cmd,path_to_libradtran_bin])
+
+    return lrt_inp_irrad, lrt_inp
+
+
+
+
+
+def lrt_to_pandas_dataframe(h,aod, altitude_km, sza, lrt_out_dir):
+    '''
+
+    Save the data to panadas
+
+    '''
+
+    # Now load in each of them into pandas to perform math.
+    df_r = pd.read_csv(f'{lrt_out_dir}/lrt_h20_{h}_aot_{aod}_alb_0_alt_{round(altitude_km*1000)}_toauu.out', sep='\s+', header=None)
+    df_r.columns = ['Wavelength','uu']
+
+    df_t = pd.read_csv(f'{lrt_out_dir}/lrt_h20_{h}_aot_{aod}_alb_0_alt_{round(altitude_km*1000)}_sureglo.out', sep='\s+', header=None)
+    df_t.columns = ['Wavelength', 'eglo']
+
+    df_s1 = pd.read_csv(f'{lrt_out_dir}/lrt_h20_{h}_aot_{aod}_alb_0.15_alt_{round(altitude_km*1000)}_sureglo.out', sep='\s+', header=None)
+    df_s1.columns = ['Wavelength', 'eglo']
+
+    df_s2 = pd.read_csv(f'{lrt_out_dir}/lrt_h20_{h}_aot_{aod}_alb_0.5_alt_{round(altitude_km*1000)}_sureglo.out', sep='\s+', header=None)
+    df_s2.columns = ['Wavelength', 'eglo']
+
+    df_irr = pd.read_csv(f'{lrt_out_dir}/lrt_h20_{h}_aot_{aod}_alt_{round(altitude_km*1000)}_IRRAD.out', sep='\s+', header=None)
+    df_irr.columns = ['Wavelength', 'edir', 'edn']
+
+    # Compute S (atmos sphere albedo)
+    df_s2['sph_alb'] = (df_s2['eglo'] - df_s1['eglo']) / (0.5 * df_s2['eglo'] -  0.15 * df_s1['eglo'])
+
+    # to one pandas dataframe
+    df = pd.DataFrame(data=df_irr['Wavelength'], columns=['Wavelength'])
+    df['l0'] = df_r['uu']
+    df['t_up'] = df_t['eglo'] / np.cos(np.radians(sza)) 
+    df['s'] = df_s2['sph_alb']
+    df['e_dir'] = df_irr['edir']
+    df['e_diff'] = df_irr['edn']
+
+    # Set units to be microW/cm2/nm/sr (common for EMIT & PRISMA)
+    df['e_dir'] = df['e_dir'] / 10
+    df['e_diff'] = df['e_diff'] / 10
+    df['l0'] = df['l0'] / 10
+
+
+    return  df