V0.1.7

README.md

@@ -36,7 +36,6 @@ pip install hyperquest
 
 
 
-
 ## Usage example
 
 - see [SNR example](tutorials/example_using_EMIT.ipynb) where different SNR methods are computed over Libya-4.

hyperquest/smile.py

@@ -34,20 +34,15 @@ def smile_metric(hdr_path, mask_waterbodies=True, no_data_value=-9999):
 
     # Mask no data values
     array[array <= no_data_value] = np.nan
-
-    # create an interior-mask for keeping only columns where it is not any NaN
-    # NOTE this may cause data loss if masking waterbodies
-    array_for_mask = np.nanmean(array[:,:,1], axis=0)
-    interior_ix = ~np.isnan(array_for_mask).any(axis=0)
 
     # get wavelengths
     w, fwhm, obs_time = read_hdr_metadata(hdr_path)
 
     # set up outputs
-    co2_mean_output = np.full(array.shape[1], fill_value=np.nan)
-    co2_std_ouput = np.full(array.shape[1], fill_value=np.nan)
-    o2_mean_output = np.full(array.shape[1], fill_value=np.nan)
-    o2_std_ouput = np.full(array.shape[1], fill_value=np.nan)
+    co2_mean = np.full(array.shape[1], fill_value=np.nan)
+    co2_std = np.full(array.shape[1], fill_value=np.nan)
+    o2_mean = np.full(array.shape[1], fill_value=np.nan)
+    o2_std = np.full(array.shape[1], fill_value=np.nan)
 
     #  first, ensure the wavelengths covered the span of o2 and co2 features
     if np.max(w) < 800:
@@ -66,26 +61,20 @@ def smile_metric(hdr_path, mask_waterbodies=True, no_data_value=-9999):
     fwhm_bar_o2 = np.nanmean([fwhm[o2_index], fwhm[o2_index+1]])
     o2_dband = (o2_b1 + o2_b2) / fwhm_bar_o2
 
-    # Compute cross-track (columnwise) means and standard deviation
-    o2_mean = np.nanmean(o2_dband, axis=0)
-    o2_std = np.nanstd(o2_dband, axis=0)
-
-    # only include interior columns (not allow all nan)
-    o2_mean_output[interior_ix] = o2_mean
-    o2_std_ouput[interior_ix] = o2_std
-
+    # Compute cross-track (columnwise) means and standard deviation (w/respect to camera)
+    o2_mean, o2_std = cross_track_stats(o2_dband)
+    o2_mean = o2_mean.flatten()
+    o2_std = o2_std.flatten()
 
     # likely has enough data to find CO2
     if np.max(w)>2100: 
         co2_b1 = array[:, :, co2_index] 
         co2_b2 = array[:, :, co2_index+1]
         fwhm_bar_co2 = np.nanmean([fwhm[co2_index], fwhm[co2_index+1]])
         co2_dband = (co2_b1 + co2_b2) / fwhm_bar_co2
-        co2_mean = np.nanmean(co2_dband, axis=0)
-        co2_std = np.nanstd(co2_dband, axis=0)
-        co2_mean_output[interior_ix] = co2_mean
-        co2_std_ouput[interior_ix] = co2_std
-
+        co2_mean, co2_std = cross_track_stats(co2_dband)
+        co2_mean = co2_mean.flatten()
+        co2_std = co2_std.flatten()
 
     return o2_mean, co2_mean, o2_std, co2_std
 
@@ -108,7 +97,7 @@ def nodd_o2a(hdr_path, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_waterbod
     array[array <= no_data_value] = np.nan
 
     # Average in down-track direction (reduce to 1 row)
-    array = np.nanmean(array, axis=0)
+    array,_ = cross_track_stats(array)
 
     # Get data from hdr
     w_sensor, fwhm, obs_time = read_hdr_metadata(hdr_path)

hyperquest/utils.py

@@ -3,7 +3,9 @@
 from os.path import abspath, exists
 from dateutil import parser
 from datetime import timezone
-
+import numpy as np
+from spectral import *
+import cv2
 
 def binning(local_mu, local_sigma, nbins):
     '''
@@ -228,4 +230,48 @@ def mask_atmos_windows(value, wavelengths):
 
     value[mask] = np.nan
 
-    return value
+    return value
+
+
+def cross_track_stats(image):
+    '''
+    TODO
+    '''
+    # get rows, cols
+    r, c = image.shape[0], image.shape[1]
+
+    # make a 3d array if not
+    if len(image.shape) == 2:
+        image = image[:, :, np.newaxis]
+
+    # get top left adn bottom left in relation to camera
+    top_left = next((x, y) for y in range(r) for x in range(c) if image[y, x, :].sum() > 0)
+    bottom_left = next((x, y) for x in range(c) for y in range(r-1, -1, -1) if image[y, x, :].sum() > 0)
+
+    # slicing direction going left to right --> 
+    dx, dy = bottom_left[0] - top_left[0], bottom_left[1] - top_left[1]
+    perp_dx, perp_dy = dy / (dx**2 + dy**2)**0.5, -dx / (dx**2 + dy**2)**0.5 
+
+    mean_along_line = []
+    std_along_line = []
+
+    for i in range(0, c):
+        # ends of linear line to sample along
+        x0, y0 = int(top_left[0] + i * perp_dx), int(top_left[1] + i * perp_dy)
+        x1, y1 = int(bottom_left[0] + i * perp_dx), int(bottom_left[1] + i * perp_dy)
+
+        # use CV2 to sample
+        mask = np.zeros_like(image[:, :, 0], dtype=np.uint8)
+        cv2.line(mask, (x0, y0), (x1, y1), color=255, thickness=1)
+        data = image[mask > 0, :]
+
+        # take mean, similar to taking it along columns , but now with respect to this line
+        mean_data = np.nanmean(data, axis=0) 
+        std_data = np.nanstd(data, axis=0)
+        mean_along_line.append(mean_data)
+        std_along_line.append(std_data)
+
+    mean_along_line = np.array(mean_along_line)
+    std_along_line = np.array(std_along_line)
+
+    return mean_along_line, std_along_line

setup.py

@@ -16,7 +16,7 @@
 
 setup(
     name="hyperquest", 
-    version="0.1.6", 
+    version="0.1.7", 
     author="Brent Wilder", 
     author_email="[email protected]",
     description=" A Python package for Hyperspectral quality estimation in hyperspectral imaging (imaging spectroscopy)",
@@ -33,7 +33,8 @@
         "joblib>=1.0.0", 
         "cython>=3.0.11",
         "spectral>=0.23.0",
-        "dem_stitcher>=2.5"
+        "dem_stitcher>=2.5",
+        "opencv-python>=4.11"
         ],
     classifiers=[
         "Programming Language :: Python :: 3",