working toward subset loading

.github/workflows/ci.yml

@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.10"]
+        python-version: ["3.10", "3.11", "3.12", "3.13"]
 
     steps:
       - name: Checkout repository

.gitignore

@@ -169,3 +169,5 @@ cython_debug/
 
 # PyPI configuration file
 .pypirc
+
+.vscode

EMITL2ARFL/EMITL2ARFL.py

@@ -10,6 +10,6 @@
 from .apply_geometry_lookup_table import *
 from .emit_ortho_raster import *
 from .emit_xarray import *
-from .extract_GLT import *
+from .extract_GLT_array import *
 from .ortho_xr import *
 from .get_pixel_center_coords import *

EMITL2ARFL/apply_geometry_lookup_table.py

@@ -3,62 +3,71 @@
 from .constants import *
 
 def apply_GLT(
-        swath_array: np.ndarray, 
-        GLT_array: np.ndarray, 
-        fill_value: int = FILL_VALUE, 
-        GLT_nodata_value: int = GLT_NODATA_VALUE) -> np.ndarray:
+            swath_array: np.ndarray, 
+            GLT_array: np.ndarray, 
+            fill_value: int = FILL_VALUE, 
+            GLT_nodata_value: int = GLT_NODATA_VALUE) -> np.ndarray:
     """
-    Applies a Geometry Lookup Table (GLT) to a numpy array representing satellite data, 
-    to orthorectify it based on the GLT. This function supports input arrays of 2 or 3 dimensions.
+    Orthorectifies satellite swath data using a Geometry Lookup Table (GLT).
+
+    This function remaps raw satellite data (swath_array) onto a georeferenced output grid using a GLT.
+    The GLT provides, for each output pixel (in geographic coordinates), the corresponding (row, column)
+    indices in the original swath data. This process corrects for geometric distortions due to sensor
+    viewing geometry, terrain, and other effects, producing an orthorectified (geospatially accurate) output.
+
+    Geospatial process overview:
+    - The GLT is a 3D array of shape (latitude, longitude, 2), where the last dimension holds the (row, column)
+        indices in the swath data for each output pixel location in the georeferenced grid.
+    - For each output pixel, if the GLT provides valid indices, the corresponding value(s) from the swath data
+        are copied to the output. If the GLT entry is invalid (nodata), the output pixel is set to a fill value.
+    - The result is an orthorectified array, spatially aligned to the geographic grid defined by the GLT.
 
     Parameters:
-    - swath_array (np.ndarray): The input satellite data array to be orthorectified. 
-                                Can be 2D (single band) or 3D (multiple bands).
-    - GLT_array (np.ndarray): The Geometry Lookup Table array, which maps the input array's 
-                              pixels to geographic locations. It is a 3-dimensional array in the 
-                              shape of (latitude, longitude, 2), with the last dimension 
-                              representing (row, column) indices.
-    - fill_value (int, optional): The value used to fill the output array wherever the GLT 
-                                  does not provide a mapping. Defaults to FILL_VALUE from constants.
-    - GLT_nodata_value (int, optional): The value in the GLT_array that indicates no data or 
-                                        invalid mapping. Pixels with this value in the GLT are 
-                                        filled with `fill_value` in the output array. Defaults to 
-                                        GLT_NODATA_VALUE from constants.
+            swath_array (np.ndarray): Raw satellite data to be orthorectified. Shape: (rows, cols) or (rows, cols, bands).
+            GLT_array (np.ndarray): Geometry Lookup Table. Shape: (latitude, longitude, 2), with (row, col) indices.
+            fill_value (int, optional): Value to use for unmapped output pixels. Defaults to FILL_VALUE.
+            GLT_nodata_value (int, optional): Value in GLT indicating invalid mapping. Defaults to GLT_NODATA_VALUE.
 
     Returns:
-    - np.ndarray: A numpy array of the same number of dimensions as `swath_array`, containing the 
-                  orthorectified data. The shape of the output array is determined by the dimensions 
-                  of the GLT_array and the number of bands in `swath_array`.
+            np.ndarray: Orthorectified data array, shape (latitude, longitude, bands).
 
     Raises:
-    - ValueError: If the dimensions of the input arrays are not compatible or if the GLT_array does 
-                  not have the last dimension of size 2.
+            ValueError: If input array dimensions are incompatible or GLT last dimension is not size 2.
     """
 
-    # Ensure GLT_array has the correct shape
+    # 1. Validate GLT shape: must be (latitude, longitude, 2) for geospatial mapping
     if GLT_array.ndim not in [2, 3] or (GLT_array.ndim == 3 and GLT_array.shape[-1] != 2):
         raise ValueError("GLT_array must be 2D or 3D with the last dimension of size 2.")
 
-    # Adjust swath_array dimensions if necessary
+    # 2. Ensure swath_array is 3D for consistent band handling
     if swath_array.ndim == 2:
+        # Convert single-band data to shape (rows, cols, 1)
         swath_array = swath_array[:, :, np.newaxis]
 
-    # Extract dimensions for the output array
+    # 3. Prepare output array shape: (latitude, longitude, bands)
     latitude_length, longitude_length = GLT_array.shape[:2]
     band_length = swath_array.shape[-1]
     ortho_array_shape = (latitude_length, longitude_length, band_length)
 
-    # Initialize the output array
+    # 4. Initialize output with fill_value (for unmapped pixels)
     ortho_array = np.full(ortho_array_shape, fill_value, dtype=np.float32)
 
-    # Identify valid GLT entries
+    # 5. Identify valid GLT entries (where both row and col indices are not nodata)
+    #    valid_GLT is a 2D boolean mask of shape (latitude, longitude)
     valid_GLT = np.all(GLT_array != GLT_nodata_value, axis=-1)
 
-    # Adjust GLT indices to zero-based
+    # 6. Convert GLT indices from 1-based to 0-based (Python convention)
     zero_based_indices = GLT_array - 1
 
-    # Apply GLT to swath_array
+    # 7. For each valid output pixel, copy the corresponding swath data using GLT indices
+    #    This remaps swath_array values to their georeferenced locations in ortho_array
+    #    - zero_based_indices[..., 0] gives swath column indices
+    #    - zero_based_indices[..., 1] gives swath row indices
+    #    - valid_GLT mask selects only valid mappings
     ortho_array[valid_GLT, :] = swath_array[zero_based_indices[valid_GLT, 1], zero_based_indices[valid_GLT, 0], :]
+
+    # 8. Replace any fill value of -9999 with np.nan for easier downstream analysis
     ortho_array = np.where(ortho_array == -9999, np.nan, ortho_array)
 
+    # 9. Return the orthorectified, geospatially aligned output array
     return ortho_array

EMITL2ARFL/constants.py

@@ -6,4 +6,6 @@
 EMIT_L2A_REFLECTANCE_DOI = "10.5067/EMIT/EMITL2ARFL.001"
 EMIT_L2A_REFLECTANCE_CONCEPT_ID = "C2408750690-LPCLOUD"
 
-DOWNLOAD_DIRECTORY = "~/data/EMIT_L2A_RFL"
+DOWNLOAD_DIRECTORY = "~/data/EMIT_L2A_RFL"
+
+QUALITY_BANDS = [0, 1, 2, 3, 4]

EMITL2ARFL/extract_GLT_array.py

@@ -0,0 +1,55 @@
+import numpy as np
+import xarray as xr
+
+from .constants import *
+
+def extract_GLT_array(swath_ds: xr.Dataset, GLT_nodata_value: int = GLT_NODATA_VALUE) -> np.ndarray:
+    """
+    Extracts the EMIT Geometry Lookup Table (GLT) index pairs from an xarray.Dataset or NetCDF file.
+
+    The GLT is a geospatial mapping array that, for each pixel in the output geographic grid, provides
+    the corresponding (row, column) indices in the original satellite swath data. This enables orthorectification:
+    the process of transforming raw sensor data into a georeferenced product by mapping each output pixel
+    to its correct location in the input data.
+
+    This function:
+    - Loads the swath dataset (if a filename is provided).
+    - Extracts the 'glt_x' and 'glt_y' arrays, which contain the swath column and row indices for each output pixel.
+    - Stacks these into a single array of shape (latitude, longitude, 2), where the last dimension holds (row, column) pairs: (row, column) = (glt_y, glt_x).
+    - Replaces any missing values (NaN) with the specified GLT_nodata_value, ensuring all output indices are valid integers.
+
+    Parameters
+    ----------
+    swath_ds : xr.Dataset | str
+        EMIT swath xarray Dataset containing 'glt_x' and 'glt_y' arrays, or a filename to load.
+    GLT_nodata_value : int, optional
+        Value to use for missing GLT indices (default: GLT_NODATA_VALUE).
+
+    Returns
+    -------
+    np.ndarray
+        Array of shape (latitude, longitude, 2) with GLT index pairs (row, column), dtype=int.
+        Missing values are set to GLT_nodata_value.
+    """
+    # Step 1: If input is a filename, load the xarray.Dataset
+    if isinstance(swath_ds, str):
+        # Local import to avoid circular import
+        from .emit_xarray import emit_xarray
+        ds: xr.Dataset = emit_xarray(swath_ds, ortho=False)
+    else:
+        ds: xr.Dataset = swath_ds
+
+    # Step 2: Extract GLT x (column) and y (row) indices from the dataset
+    # These arrays map each output pixel to its location in the original swath
+    GLT_x: np.ndarray = ds["glt_x"].data  # swath column indices
+    GLT_y: np.ndarray = ds["glt_y"].data  # swath row indices
+
+    # Step 3: Stack row and column indices into a single array of shape (latitude, longitude, 2)
+    # The last dimension holds (row, column) pairs for each output pixel: (row, column) = (glt_y, glt_x)
+    GLT_array: np.ndarray = np.nan_to_num(
+        np.stack([GLT_y, GLT_x], axis=-1),
+        nan=GLT_nodata_value
+    ).astype(int)
+
+    # Step 4: Return the GLT array, ready for use in geospatial orthorectification
+    return GLT_array

EMITL2ARFL/find_EMIT_L2A_RFL_granule.py

@@ -0,0 +1,30 @@
+from .search_EMIT_L2A_RFL_granules import search_EMIT_L2A_RFL_granules
+import earthaccess
+
+def find_EMIT_L2A_RFL_granule(granule: earthaccess.search.DataGranule = None, orbit: int = None, scene: int = None):
+    """
+    Find an EMIT L2A Reflectance granule by granule object or by orbit and scene.
+
+    Args:
+        granule (earthaccess.search.DataGranule, optional): The granule to retrieve. Defaults to None.
+        orbit (int, optional): The orbit number to search for the granule. Defaults to None.
+        scene (int, optional): The scene number to search for the granule. Defaults to None.
+
+    Returns:
+        earthaccess.search.DataGranule: The found granule.
+
+    Raises:
+        ValueError: If no granule is found for the provided orbit and scene, or if neither granule nor orbit/scene are provided.
+    """
+    if granule is None and orbit is not None and scene is not None:
+        remote_granules = search_EMIT_L2A_RFL_granules(orbit=orbit, scene=scene)
+
+        if len(remote_granules) == 0:
+            raise ValueError(f"no EMIT L2A RFL granule found for orbit {orbit} and scene {scene}")
+
+        granule = remote_granules[0]
+
+    if granule is None:
+        raise ValueError("either granule or orbit and scene must be provided")
+
+    return granule

EMITL2ARFL/granule.py

@@ -1,10 +1,17 @@
 from glob import glob
 from os.path import join, abspath, dirname, expanduser
+from typing import List
 
-from rasters import Raster, RasterGeometry
+import numpy as np
 
+import netCDF4
+
+from rasters import Raster, RasterGeometry, RasterGeolocation
+
+from .constants import *
 from .emit_ortho_raster import emit_ortho_raster
 from .quality_mask import quality_mask
+from .extract_GLT_array import extract_GLT_array
 
 class EMITL2ARFL:
     def __init__(self, directory: str):
@@ -33,21 +40,47 @@ def mask_filename(self) -> str:
     def uncertainty_filename(self) -> str:
         return glob(join(self.directory_absolute, "*_RFLUNCERT_*.nc"))[0]
 
-    def quality_mask(self, geometry: RasterGeometry) -> Raster:
-        raster = quality_mask(
+    @property
+    def lat(self) -> np.ndarray:
+        # read the `lat` array from the `location` group in the reflectance NetCDF file    
+        with netCDF4.Dataset(self.reflectance_filename, "r") as ds:
+            lat = ds.groups["location"].variables["lat"][:]
+
+        return lat
+
+    @property
+    def lon(self) -> np.ndarray:
+        # read the `lon` array from the `location` group in the reflectance NetCDF file    
+        with netCDF4.Dataset(self.reflectance_filename, "r") as ds:
+            lon = ds.groups["location"].variables["lon"][:]
+
+        return lon
+
+    @property
+    def geolocation(self) -> RasterGeolocation:
+        return RasterGeolocation(
+            x=self.lon,
+            y=self.lat
+        )
+
+    def GLT(self) -> Raster:
+        return Raster(extract_GLT_array(swath_ds=self.reflectance_filename), geometry=self.geolocation)
+
+    def quality_mask(self, quality_bands: List[int] = QUALITY_BANDS) -> np.ndarray:
+        qmask = quality_mask(
             filepath=self.mask_filename,
-            quality_bands=[0, 1, 2, 3, 4]
+            quality_bands=quality_bands
         )
 
-        if geometry is not None:
-            raster = raster.to_geometry(geometry)
-        
-        return raster
+        return qmask
+
+    def reflectance(self, geometry: RasterGeometry = None) -> Raster:
+        qmask = self.quality_mask()
 
-    def reflectance(self, geometry: RasterGeometry) -> Raster:
         raster = emit_ortho_raster(
             filepath=self.reflectance_filename,
-            layer_name="reflectance"
+            layer_name="reflectance",
+            qmask=qmask
         )
 
         if geometry is not None:

EMITL2ARFL/ortho_xr.py

@@ -2,7 +2,7 @@
 import xarray as xr
 
 from .constants import *
-from .extract_GLT import extract_GLT
+from .extract_GLT_array import extract_GLT_array
 from .apply_geometry_lookup_table import apply_GLT
 from .get_pixel_center_coords import get_pixel_center_coords
 
@@ -19,7 +19,7 @@ def ortho_xr(swath_ds: xr.Dataset, GLT_nodata_value: int = GLT_NODATA_VALUE, fil
     ortho_ds: an orthocorrected xarray dataset.
     """
     # extract GLT
-    GLT_array = extract_GLT(swath_ds)
+    GLT_array = extract_GLT_array(swath_ds=swath_ds)
 
     # List Variables
     var_list = list(swath_ds.data_vars)

EMITL2ARFL/quality_mask.py

@@ -1,7 +1,13 @@
+from typing import List
 import numpy as np
 import xarray as xr
 
-def quality_mask(filepath, quality_bands):
+from .constants import *
+
+def quality_mask(
+        filepath: str, 
+        quality_bands: List[str] = QUALITY_BANDS, 
+        engine: str = ENGINE) -> np.ndarray:
     """
     This function builds a single layer mask to apply based on the bands selected from an EMIT L2A Mask file.
 
@@ -13,19 +19,24 @@ def quality_mask(filepath, quality_bands):
     qmask: a numpy array that can be used with the emit_xarray function to apply a quality mask.
     """
     # Open Dataset
-    mask_ds = xr.open_dataset(filepath, engine="h5netcdf")
+    mask_ds = xr.open_dataset(filepath, engine=ENGINE)
+
     # Open Sensor band Group
     mask_parameters_ds = xr.open_dataset(
-        filepath, engine="h5netcdf", group="sensor_band_parameters"
+        filepath, 
+        engine=engine, 
+        group="sensor_band_parameters"
     )
+
     # Print Flags used
     flags_used = mask_parameters_ds["mask_bands"].data[quality_bands]
-    print(f"Flags used: {flags_used}")
+
     # Check for data bands and build mask
     if any(x in quality_bands for x in [5, 6]):
         err_str = f"Selected flags include a data band (5 or 6) not just flag bands"
         raise AttributeError(err_str)
     else:
         qmask = np.sum(mask_ds["mask"][:, :, quality_bands].values, axis=-1)
         qmask[qmask > 1] = 1
-    return qmask
+
+    return qmask