V0.1.13 patch

README.md

@@ -10,11 +10,13 @@
 
 __Important: this package assumes the following about input hyperspectral data:__ 
 - Data must be in NetCDF (.nc) or ENVI (.hdr)
-- Data not georeferenced (typically referred to as L1B before ortho)
-- Data in radiance (assumed microW/cm2/nm/sr (for now))
+- Currently data is expected in Radiance.
+  - For smile & striping methods, data must not be georeferenced (typically referred to as L1B before ortho)
 - Pushbroom imaging spectrometer, such as, but not limited to:
     - AVIRIS-NG, AVIRIS-3, DESIS, EnMAP, EMIT, GaoFen-5, HISUI, Hyperion EO-1, HySIS, PRISMA, Tanager-1
 
+NOTE: this is under active development. It is important to note that noise methods shown here do not account for spectrally correlated noise. This is a work in progress as I digest literature and translate into python.
+
 ## Installation Instructions
 The latest release can be installed via pip:
 

hyperquest/coregistration.py

@@ -7,7 +7,7 @@
 
 def sub_pixel_shift(path_to_data, band_index_vnir, band_index_vswir, no_data_value=-9999, upsample_factor=5000):
     '''
-    A wrapper function for skimage.registration's `phase_cross_correlation`
+    A wrapper function for skimage.registration's `phase_cross_correlation` 
 
     Parameters
     ----------

hyperquest/intrinsic_dimensionality.py

@@ -48,12 +48,12 @@ def random_matrix_theory(path_to_data, noise_variance, mask_waterbodies=True, al
         raise Exception('Data needs to be a 3D array.')
 
     # ensure data is hyperspectral 
-    if (np.max(w) - np.min(w)) / len(w) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w) - np.min(w)) / len(w) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # mask waterbodies
     if mask_waterbodies is True:
-        img = mask_water_using_ndwi(img, w)
+        img = mask_water_using_ndwi(img, w, no_data_value=no_data_value)
 
     # Mask no data values
     img[img <= no_data_value] = np.nan
@@ -86,22 +86,23 @@ def random_matrix_theory(path_to_data, noise_variance, mask_waterbodies=True, al
 
     # Get the noise covariance matrix of size bands x bands
     noise_variance[np.isnan(noise_variance)] = 0
-    noise_variance  = np.expand_dims(noise_variance, axis=1)
-    N = np.diag(np.diag(noise_variance.T * noise_variance / bands)) # similar to https://github.com/bearshng/LRTFL0/blob/master/estNoise.m
+    N = np.diag(noise_variance.T * noise_variance / bands) # similar to https://github.com/bearshng/LRTFL0/blob/master/estNoise.m
     N[np.isnan(N)] = 0
-    
+
     # Eigenvectors and eigenvalues of S
     # Note these are opposite of MATLAB [evec_S,eval_S]
     eval_S, evec_S = np.linalg.eig(S)
-    sortind = np.argsort(eval_S)[::-1]
-    eval_S = eval_S[sortind]
-    evec_S = evec_S[:,sortind]
+    eval_S = np.diag(eval_S)
+    sortind = np.argsort(-eval_S)
+    evec_S = evec_S[np.arange(evec_S.shape[0])[:,None], sortind]
+    eval_S = eval_S[np.arange(eval_S.shape[0])[:,None], sortind]
 
     # Eigenvectors and eigenvalues of Pi = S-N
     eval_Pi, evec_Pi = np.linalg.eig(S-N) 
-    sortind2 = np.argsort(eval_Pi)[::-1]
-    eval_Pi = eval_Pi[sortind2]
-    evec_Pi = evec_Pi[:,sortind2]
+    eval_Pi = np.diag(eval_Pi)
+    sortind2 = np.argsort(-eval_Pi)
+    evec_Pi = evec_Pi[np.arange(evec_Pi.shape[0])[:,None], sortind2]
+    eval_Pi = eval_Pi[np.arange(eval_Pi.shape[0])[:,None], sortind2]
 
     # there is a different threshold for each band to represent noise conditions
     X = []
@@ -110,9 +111,8 @@ def random_matrix_theory(path_to_data, noise_variance, mask_waterbodies=True, al
     X = np.array(X)
 
     # The ID is the number of eigenvalues exceeding the threshold X.T*sigma
-    thresholded = np.maximum(0, eval_S - X * sigma)
-    intersect_i = np.argwhere(thresholded > 0)
+    eval_S = np.diag(eval_S) # put back to 1d
+    intersect_i = np.argwhere((eval_S > X.T * sigma) & (eval_S - X.T * sigma >= 0))
     K_est = len(intersect_i)
 
-
-    return K_est
+    return K_est

hyperquest/mlr.pyx

@@ -30,6 +30,11 @@ def mlr_spectral(np.ndarray[np.float64_t, ndim=2] block):
 
     # for k in range of wavelengths (except first and last)
     for k in range(1, cols - 1):
+
+        # ensure not using bad band if reflectance product
+        if block[0,k] <= -0.01 or block[0,k-1] <= -0.01 or block[0,k+1] <= -0.01:
+            continue
+
         # create the X and y for MLR
         X = np.vstack([block[:, k - 1], block[:, k + 1]]).T
         y = block[:, k]
@@ -48,6 +53,7 @@ def mlr_spectral(np.ndarray[np.float64_t, ndim=2] block):
 def mlr_spectral_spatial(np.ndarray[np.float64_t, ndim=2] block):
     '''
     TODO
+
     '''
 
     # remove data that is NaN
@@ -69,6 +75,11 @@ def mlr_spectral_spatial(np.ndarray[np.float64_t, ndim=2] block):
 
     # for k in range of wavelengths (except first and last)
     for k in range(1, cols - 1):
+
+        # ensure not using bad band if reflectance product
+        if block[0,k] <= -0.01 or block[0,k-1] <= -0.01 or block[0,k+1] <= -0.01:
+            continue
+
         # create the X and y for MLR
         X = np.vstack([block[:, k - 1], block[:, k + 1]]).T
         neighbor_k = np.roll(block[:, k], shift=1)  # Shift 1 to find a neighbor pixel
@@ -79,7 +90,7 @@ def mlr_spectral_spatial(np.ndarray[np.float64_t, ndim=2] block):
             coef = np.linalg.lstsq(X, y, rcond=None)[0]
             y_pred = X @ coef
 
-            # 3 DOF because of MLR
+            # 4 DOF because of MLR
             sigma_block[k] = np.std(y - y_pred, ddof=4)
             mu_block[k] = np.mean(y)
 

hyperquest/smile.py

@@ -52,7 +52,7 @@ def smile_metric(path_to_data, rotate, mask_waterbodies=True):
         raise Exception('Data needs to be a 3D array.')
 
     # ensure data is hyperspectral 
-    if (np.max(w) - np.min(w)) / len(w) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w) - np.min(w)) / len(w) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # Perform rotation if needed
@@ -67,7 +67,7 @@ def smile_metric(path_to_data, rotate, mask_waterbodies=True):
 
     # mask waterbodies
     if mask_waterbodies is True:
-        array = mask_water_using_ndwi(array, w)
+        array = mask_water_using_ndwi(array, w, no_data_value=np.nan)
 
     # set up outputs
     co2_mean = np.full(array.shape[1], fill_value=np.nan)
@@ -169,7 +169,7 @@ def nodd_o2a(path_to_data, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, m
         w_sensor, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure data is hyperspectral 
-    if (np.max(w_sensor) - np.min(w_sensor)) / len(w_sensor) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w_sensor) - np.min(w_sensor)) / len(w_sensor) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # ensure data has wavelength range in O2-A band
@@ -196,7 +196,7 @@ def nodd_o2a(path_to_data, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, m
 
     # mask waterbodies
     if mask_waterbodies is True:
-        array = mask_water_using_ndwi(array, w_sensor)
+        array = mask_water_using_ndwi(array, w_sensor, no_data_value=np.nan)
 
     # Average in down-track direction (reduce to 1 row)
     array = np.nanmean(array, axis=0)

hyperquest/snr.py

@@ -12,6 +12,8 @@ def rlsd(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_w
     '''
     Residual-scaled local standard deviation (Gao et al., 2007).
 
+    Neighbor pixel not included in MLR, ( p_k-1 , p_k+1).
+
     Parameters
     ----------
     path_to_data : str
@@ -48,7 +50,7 @@ def rlsd(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_w
         w, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure data is hyperspectral 
-    if (np.max(w) - np.min(w)) / len(w) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w) - np.min(w)) / len(w) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # mask waterbodies
@@ -101,7 +103,9 @@ def rlsd(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_w
 
 def ssdc(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_waterbodies=True, no_data_value=-9999):
     '''
-    Spectral and spatial de-correlation (Roger & Arnold, 1996)
+    Spectral and spatial de-correlation (Roger & Arnold, 1996).
+
+    Neighbor pixel with band K is included in MLR, ( p_k-1 , p_k+1, p+1_k ).
 
     Parameters
     ----------
@@ -139,12 +143,12 @@ def ssdc(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_w
         w, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure data is hyperspectral 
-    if (np.max(w) - np.min(w)) / len(w) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w) - np.min(w)) / len(w) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # mask waterbodies
     if mask_waterbodies is True:
-        array = mask_water_using_ndwi(array, w)
+        array = mask_water_using_ndwi(array, w, no_data_value=no_data_value)
 
     # Mask no data values
     array[array <= no_data_value] = np.nan
@@ -190,10 +194,13 @@ def ssdc(path_to_data, block_size, nbins=150, ncpus=1, snr_in_db = False, mask_w
     return snr, noise_variance
 
 
-def hrdsdc(path_to_data, n_segments=200, compactness=0.1, n_pca=3, ncpus=1, include_neighbor_pixel_in_mlr=True,
+def hrdsdc(path_to_data, n_segments=200, compactness=0.1, n_pca=3, ncpus=1, 
            snr_in_db=False, mask_waterbodies=True, no_data_value=-9999):
     '''
-    Homogeneous regions division and spectral de-correlation (Gao et al., 2008)
+    Homogeneous regions division and spectral de-correlation (Gao et al., 2008). 
+
+    Neighbor pixel with band K is included in MLR, ( p_k-1 , p_k+1, p+1_k ).
+
 
     Parameters
     ----------
@@ -233,12 +240,12 @@ def hrdsdc(path_to_data, n_segments=200, compactness=0.1, n_pca=3, ncpus=1, incl
         w, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure data is hyperspectral 
-    if (np.max(w) - np.min(w)) / len(w) < 50: # assume hyperspectral data not coarser than 50 nm spec res
+    if (np.max(w) - np.min(w)) / len(w) > 50: # assume hyperspectral data not coarser than 50 nm spec res
         raise Exception('Data needs to be a hyperspectral image.')
 
     # mask waterbodies
     if mask_waterbodies is True:
-        array = mask_water_using_ndwi(array, w)
+        array = mask_water_using_ndwi(array, w, no_data_value=no_data_value)
 
     # Mask no data values (to negative 9999 for PCA and SLIC to work)
     array[array <= no_data_value] = -9999
@@ -261,24 +268,19 @@ def hrdsdc(path_to_data, n_segments=200, compactness=0.1, n_pca=3, ncpus=1, incl
     def process_segment(u):
         test_mask = (segments == u)
         test_segment = array[test_mask]
-        test_segment = test_segment[test_segment[:, 0] > -99]
+        test_segment = test_segment[test_segment[:, 0] > -99] #spatial nan
         if test_segment.shape[0] != 0:
             return test_segment
         else:
             return None
     segment_data = Parallel(n_jobs=ncpus)(delayed(process_segment)(u) for u in unique_segments)
     segment_data = [seg for seg in segment_data if seg is not None]
 
-    # Parallel processing of all segments depending on method selected
-    if include_neighbor_pixel_in_mlr == False:
-        # Perform just spectral MLR
-        results = Parallel(n_jobs=ncpus, 
-                           timeout=None)(delayed(mlr_spectral)(segment) for segment in segment_data)
-
-    else: # perform spectral-spatial MLR using k` nearby neighbor.
-        results = Parallel(n_jobs=ncpus, 
-                           timeout=None)(delayed(mlr_spectral_spatial)(segment) for segment in segment_data) 
-
+    # Parallel processing of all segments 
+    # perform spectral-spatial MLR using k` nearby neighbor.
+    results = Parallel(n_jobs=ncpus, 
+                        timeout=None)(delayed(mlr_spectral_spatial)(segment) for segment in segment_data) 
+
     # Aggregate results
     local_mu = np.array([res[0] for res in results])
     local_sigma = np.array([res[1] for res in results])

hyperquest/striping.py

@@ -6,7 +6,7 @@
 
 def sigma_threshold(path_to_data, rotate, sigma_multiplier = 3):
     '''
-    Uses a sigma threshold counting neighboring pixels to determine if striping is present, as presented in,
+    Uses a sigma threshold counting neighboring pixels to determine if striping is present, as presented in, 
 
     Yokoya 2010, Preprocessing of hyperspectral imagery with consideration of smile and keystone properties,
 

hyperquest/utils.py

@@ -195,7 +195,11 @@ def retrieve_data_from_nc(path_to_data):
     ds = h5netcdf.File(path_to_data, mode='r')
 
     # get radiance and fhwm and wavelength
-    array = np.array(ds['radiance'], dtype=np.float64)
+    if 'rad' in path_to_data.lower():
+        array = np.array(ds['radiance'], dtype=np.float64)
+    else:
+        array = np.array(ds['reflectance'], dtype=np.float64)
+
     fwhm = np.array(ds['sensor_band_parameters']['fwhm'][:].data.tolist(), dtype=np.float64)
     wave = np.array(ds['sensor_band_parameters']['wavelengths'][:].data.tolist(), dtype=np.float64)
 
@@ -208,7 +212,7 @@ def linear_to_db(snr_linear):
     return 10 * np.log10(snr_linear)
 
 
-def mask_water_using_ndwi(array, wavelengths, ndwi_threshold=0.25):
+def mask_water_using_ndwi(array, wavelengths, no_data_value = -9999, ndwi_threshold=0.25):
     '''
     Returns array where NDWI greater than a threshold are set to NaN.
 
@@ -228,7 +232,7 @@ def mask_water_using_ndwi(array, wavelengths, ndwi_threshold=0.25):
     nir = array[:, :, nir_index] 
     ndwi = (green - nir) / (green + nir)
 
-    array[(ndwi > ndwi_threshold)] = np.nan
+    array[(ndwi > ndwi_threshold)] = no_data_value
 
     return array
 
@@ -245,7 +249,7 @@ def mask_atmos_windows(spectra, wavelengths):
         spectra: TOA radiance data but with noisy atmospheric bands masked out
     '''
 
-    mask = ((wavelengths >= 1250) & (wavelengths <= 1450)) | ((wavelengths >= 1780) & (wavelengths <= 1950))
+    mask = ((wavelengths >= 1250) & (wavelengths <= 1450)) | ((wavelengths >= 1750) & (wavelengths <= 1970))
 
     spectra[mask] = np.nan
 

setup.py

@@ -5,7 +5,6 @@
 with open("README.md", "r") as fh:
     long_description = fh.read()
 
-
 ext_modules = [
     Extension(
         'hyperquest.mlr',
@@ -16,7 +15,7 @@
 
 setup(
     name="hyperquest", 
-    version="0.1.12", 
+    version="0.1.13", 
     author="Brent Wilder", 
     author_email="[email protected]",
     description=" A Python package for Hyperspectral quality estimation in hyperspectral imaging (imaging spectroscopy)",