Merge pull request #9 from brentwilder/v0.1.10

V0.1.10
brentwilder · Feb 15, 2025 · 9196347 · 9196347
2 parents 62f8525 + 02b4a26
commit 9196347
Show file tree

Hide file tree

Showing 13 changed files with 217 additions and 479 deletions.
diff --git a/README.md b/README.md
@@ -10,7 +10,7 @@
 
 __Important: this package assumes the following about input hyperspectral data:__ 
 - Radiance (if using radiative transfer model must be in microW/cm2/nm/sr (for now))
-- ENVI format with a .HDR file
+- Data must be in NetCDF (.nc) or ENVI (.hdr)
 - Pushbroom imaging spectrometer, such as:
     - AVIRIS-Classic, AVIRIS-NG, AVIRIS-3, DESIS, EnMAP, EMIT, GaoFen-5, HISUI, Hyperion EO-1, HySIS, PRISMA, SEBASS, Tanager
 - For smile estimation you must use data that has not been georeferenced.
@@ -39,7 +39,7 @@ pip install hyperquest
 
 ## Usage example
 
-- see [SNR example](tutorials/example_using_EMIT.ipynb) where different SNR methods are computed over Libya-4.
+- see [SNR example](tutorials/example_using_EMIT.ipynb) and  [Smile example](tutorials/testing_smile_methods.ipynb) which has different methods computed over Libya-4.
 
 ## libRadtran install instructions
 - Can be installed on Unix type system using the following link:

diff --git a/hyperquest/coregistration.py b/hyperquest/coregistration.py
@@ -6,14 +6,14 @@
 
 
 
-def sub_pixel_shift(hdr_path, band_index_vnir, band_index_vswir, no_data_value=-9999, upsample_factor=5000):
+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`
 
 
 
     Parameters:
-        hdr_path (str): Path to the .hdr file..
+        path_to_data (str): Path to the .hdr or .nc
         band_index_vnir (int): Band index for VNIR camera , assuming the first band is 0.
         band_index_vswir (int): Band index for VSWIR camera , assuming the first band is 0.
         no_data_value (int): Assumed to be -9999.
@@ -23,9 +23,13 @@ def sub_pixel_shift(hdr_path, band_index_vnir, band_index_vswir, no_data_value=-
         Tuple containing shift in the X direction and shift in the Y direction (in pixels)
     '''
 
-    # Load image data
-    img_path = get_img_path_from_hdr(hdr_path)
-    array = np.array(envi.open(hdr_path, img_path).load(), dtype=np.float64)
+    # Identify data type
+    if path_to_data.lower().endswith('.nc'):
+        array, _, _, _ = retrieve_data_from_nc(path_to_data)
+    else:
+        # Load raster
+        img_path = get_img_path_from_hdr(path_to_data)
+        array = np.array(envi.open(path_to_data, img_path).load(), dtype=np.float64)
 
     # Select the desired bands (VNIR and VSWIR)
     vnir_band = array[:, :, band_index_vnir]

diff --git a/hyperquest/libradtran.py b/hyperquest/libradtran.py
@@ -28,15 +28,15 @@ def get_libradtran_install_path(user_input_path):
 
 
 
-def get_libradtran_output_dir(hdr_path):
+def get_libradtran_output_dir(data_path):
     '''
     TODO
     '''
 
-    filename = os.path.splitext(os.path.basename(hdr_path))[0]
+    filename = os.path.splitext(os.path.basename(data_path))[0]
 
     # directory where hdr_path is located
-    parent_dir = os.path.dirname(os.path.abspath(hdr_path))
+    parent_dir = os.path.dirname(os.path.abspath(data_path))
 
     # Define the rtm output directory
     lrt_out_dir = os.path.join(parent_dir, f'rtm-{filename}')

diff --git a/hyperquest/radiative_transfer.py b/hyperquest/radiative_transfer.py
@@ -1,7 +1,4 @@
-from dem_stitcher import stitch_dem
-import rasterio as rio
 from pysolar import solar
-from rasterio.warp import transform_bounds
 from joblib import Parallel, delayed
 from os.path import abspath
 import subprocess
@@ -10,45 +7,34 @@
 from .utils import *
 
 def run_libradtran(h2o_mm, aod_at_550nm, sensor_zenith_angle, sensor_azimith_angle,
-                   hdr_path, libradtran_path, ncpus=1, o3_DU=300, albedo=0.15):
+                   path_to_data, average_elevation_meters, lat, lon, libradtran_path, ncpus=1, o3_DU=300, albedo=0.15):
     '''
     TODO
     '''
 
     # Get absolute path
-    hdr_path = abspath(hdr_path)
+    path_to_data = abspath(path_to_data)
     libradtran_path = abspath(libradtran_path)
 
     # path_to_libradtran_install
     # get abs, get bin directory... throw error if not found
     path_to_libradtran_bin = get_libradtran_install_path(libradtran_path)
 
-    # Get data from hdr
-    wavelength, fwhm, obs_time = read_hdr_metadata(hdr_path)
+    # Identify data type
+    if path_to_data.lower().endswith('.nc'):
+        _, _, _, obs_time = retrieve_data_from_nc(path_to_data)
+    else:
+        # get wavelengths
+        _, _, obs_time = read_hdr_metadata(path_to_data)
+
+    # convert to doy
     doy = obs_time.timetuple().tm_yday
 
     # path to where runs are saved
-    lrt_out_dir = get_libradtran_output_dir(hdr_path)
-
-    # Get bounds
-    img_path = get_img_path_from_hdr(hdr_path)
-    with rio.open(img_path) as dataset:
-        bounds_utm = dataset.bounds
-        # to EPSG:4326 , as xmin, ymin, xmax, ymax in epsg:4326
-        bounding_box = transform_bounds(dataset.crs, 'EPSG:4326', *bounds_utm)
+    lrt_out_dir = get_libradtran_output_dir(path_to_data)
 
-    # Get Copernicus DEM data based on bounding box in hdr 
-    # (assume mus = mu0 ; flat assumption) ...  X is an mxn numpy array
-    X, _ = stitch_dem(bounding_box, dem_name='glo_90')
-    X = X[:, :].flatten()
-    X = X[X < 8848] #mt everest
-    X = X[X > -430] #deadsea
-    X = X[~np.isnan(X)]
-    altitude_km = np.nanmean(X) / 1000
-
-    # Get average lat and lon
-    lon = np.mean([bounding_box[0], bounding_box[2]])
-    lat = np.mean([bounding_box[1], bounding_box[3]])
+    # average altitude in km
+    altitude_km = average_elevation_meters / 1000
 
     # use pysolar compute saa and sza
     phi0 = solar.get_azimuth(lat,lon, obs_time)

diff --git a/hyperquest/smile.py b/hyperquest/smile.py
@@ -7,7 +7,7 @@
 from .optimization import *
 
 
-def smile_metric(hdr_path, rotate, mask_waterbodies=True):
+def smile_metric(path_to_data, rotate, mask_waterbodies=True):
     '''
     This is an exact remake of the MATLAB method smileMetric(), https://www.mathworks.com/help/images/ref/smilemetric.html.
 
@@ -18,7 +18,7 @@ def smile_metric(hdr_path, rotate, mask_waterbodies=True):
     IEEE Transactions on Geoscience and Remote Sensing, 48(6), 2603-2612.
 
     Parameters: 
-        hdr_path (str): Path to the .hdr file.
+        path_to_data (str): Path to the .hdr or .nc
         rotate (int): rotate counter clockwise, either 0, 90, 180, or 270.
         mask_waterbodies (bool, optional): Whether to mask water bodies based on NDWI threshold of 0. Default is True.
 
@@ -27,9 +27,16 @@ def smile_metric(hdr_path, rotate, mask_waterbodies=True):
 
     '''
 
-    # Load raster
-    img_path = get_img_path_from_hdr(hdr_path)
-    array = np.array(envi.open(hdr_path, img_path).load(), dtype=np.float64)
+    # Identify data type
+    if path_to_data.lower().endswith('.nc'):
+        array, fwhm, w, obs_time = retrieve_data_from_nc(path_to_data)
+    else:
+        # Load raster
+        img_path = get_img_path_from_hdr(path_to_data)
+        array = np.array(envi.open(path_to_data, img_path).load(), dtype=np.float64)
+
+        # get wavelengths
+        w, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure this is the raw data and has not been georeferenced.
     if (array[:,:,0]<0).sum() > 0:
@@ -43,19 +50,16 @@ def smile_metric(hdr_path, rotate, mask_waterbodies=True):
     if rotate != 0 and rotate != 90 and rotate != 180 and rotate != 270:
         raise ValueError('rotate must be 90, 180, or 270.')
     if rotate>=90:
-        array = np.rot90(array, axes=(0,1))
-        if rotate>=180:
-            array = np.rot90(array, axes=(0,1))
-            if rotate>=270:
-                array = np.rot90(array, axes=(0,1))
+        array = np.rot90(array, axes=(0,1), k=1)
+    elif rotate>=180:
+        array = np.rot90(array, axes=(0,1), k=2)
+    elif rotate>=270:
+        array = np.rot90(array, axes=(0,1), k=3)
 
     # mask waterbodies
     if mask_waterbodies is True:
-        array = mask_water_using_ndwi(array, hdr_path)
+        array = mask_water_using_ndwi(array, w)
 
-    # get wavelengths
-    w, fwhm, obs_time = read_hdr_metadata(hdr_path)
-
     # set up outputs
     co2_mean = np.full(array.shape[1], fill_value=np.nan)
     co2_std = np.full(array.shape[1], fill_value=np.nan)
@@ -99,7 +103,7 @@ def smile_metric(hdr_path, rotate, mask_waterbodies=True):
     return o2_mean, co2_mean, o2_std, co2_std
 
 
-def nodd_o2a(hdr_path, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_waterbodies=True):
+def nodd_o2a(path_to_data, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_waterbodies=True):
     '''
     Similar to method in Felde et al. (2003) to solve for nm shift at O2-A across-track. Requires radiative transfer model run.
 
@@ -124,7 +128,7 @@ def nodd_o2a(hdr_path, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_
 
     
     Parameters: 
-        hdr_path (str): Path to the .hdr file.
+        path_to_data (str): Path to the .hdr or .nc
         rotate (int): rotate counter clockwise, either 0, 90, 180, or 270.
         path_to_rtm_output_csv (str): Path to output from radiative transfer.
         ncpus (int, optional): Number of CPUs for parallel processing. Default is 1.
@@ -136,9 +140,16 @@ def nodd_o2a(hdr_path, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_
 
     '''
 
-    # Load raster
-    img_path = get_img_path_from_hdr(hdr_path)
-    array = np.array(envi.open(hdr_path, img_path).load(), dtype=np.float64)
+    # Identify data type
+    if path_to_data.lower().endswith('.nc'):
+        array, fwhm, w_sensor, obs_time = retrieve_data_from_nc(path_to_data)
+    else:
+        # Load raster
+        img_path = get_img_path_from_hdr(path_to_data)
+        array = np.array(envi.open(path_to_data, img_path).load(), dtype=np.float64)
+
+        # get wavelengths
+        w_sensor, fwhm, obs_time = read_hdr_metadata(path_to_data)
 
     # ensure this is the raw data and has not been georeferenced.
     if (array[:,:,0]<0).sum() > 0:
@@ -152,18 +163,19 @@ def nodd_o2a(hdr_path, rotate, path_to_rtm_output_csv, ncpus=1,rho_s=0.15, mask_
     if rotate != 0 and rotate != 90 and rotate != 180 and rotate != 270:
         raise ValueError('rotate must be 90, 180, or 270.')
     if rotate>=90:
-        array = np.rot90(array, axes=(0,1))
-        if rotate>=180:
-            array = np.rot90(array, axes=(0,1))
-            if rotate>=270:
-                array = np.rot90(array, axes=(0,1))
+        array = np.rot90(array, axes=(0,1), k=1)
+    elif rotate>=180:
+        array = np.rot90(array, axes=(0,1), k=2)
+    elif rotate>=270:
+        array = np.rot90(array, axes=(0,1), k=3)
+
+    # mask waterbodies
+    if mask_waterbodies is True:
+        array = mask_water_using_ndwi(array, w_sensor)
 
     # Average in down-track direction (reduce to 1 row)
     array = np.nanmean(array, axis=0)
 
-    # Get data from hdr
-    w_sensor, fwhm, obs_time = read_hdr_metadata(hdr_path)
-
     # Only include window for o2-a
     window = (w_sensor >= 730) & (w_sensor <= 790)
     w_sensor = w_sensor[window]