encapsulating GLT

EMITL2ARFL/GLT.py

@@ -0,0 +1,30 @@
+
+from rasters import Raster, RasterGeolocation
+import numpy as np
+
+class GLT(np.ndarray):
+    """
+    Numpy ndarray subclass for the GLT array, with geolocation attribute and rows/cols properties.
+    """
+    def __new__(cls, GLT_array: np.ndarray, geolocation: RasterGeolocation = None):
+        # Cast input array to GLT subclass
+        if not isinstance(geolocation, RasterGeolocation):
+            raise TypeError(f"geolocation must be a RasterGeolocation, got {type(geolocation)}")
+        obj = np.asarray(GLT_array).view(cls)
+        obj.geolocation = geolocation
+        return obj
+
+    def __array_finalize__(self, obj):
+        if obj is None:
+            return
+        self.geolocation = getattr(obj, 'geolocation', None)
+
+    @property
+    def rows(self) -> Raster:
+        """Return the row indices (glt_y) as a Raster object with geometry from geolocation."""
+        return Raster(self[..., 0], geometry=self.geolocation)
+
+    @property
+    def cols(self) -> Raster:
+        """Return the column indices (glt_x) as a Raster object with geometry from geolocation."""
+        return Raster(self[..., 1], geometry=self.geolocation)

EMITL2ARFL/extract_GLT.py

@@ -0,0 +1,22 @@
+import xarray as xr
+
+from .read_geolocation import read_geolocation
+from .extract_GLT_array import extract_GLT_array
+from .GLT import GLT
+from .emit_xarray import emit_xarray
+from .constants import GLT_NODATA_VALUE
+
+def extract_GLT(swath_dataset: xr.Dataset, GLT_nodata_value: int = GLT_NODATA_VALUE):
+    """
+    Wrapper function to extract a GLT object from a dataset or filename.
+    """
+    # If input is a filename, load the xarray.Dataset
+    if isinstance(swath_dataset, str):
+        ds: xr.Dataset = emit_xarray(swath_dataset, ortho=False)
+    else:
+        ds: xr.Dataset = swath_dataset
+
+    GLT_array = extract_GLT_array(ds, GLT_nodata_value)
+    geolocation = read_geolocation(swath_dataset)
+
+    return GLT(GLT_array=GLT_array, geolocation=geolocation)

EMITL2ARFL/extract_GLT_array.py

@@ -3,7 +3,7 @@
 
 from .constants import *
 
-def extract_GLT_array(swath_ds: xr.Dataset, GLT_nodata_value: int = GLT_NODATA_VALUE) -> np.ndarray:
+def extract_GLT_array(swath_dataset: 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.
 
@@ -32,12 +32,12 @@ def extract_GLT_array(swath_ds: xr.Dataset, GLT_nodata_value: int = GLT_NODATA_V
         Missing values are set to GLT_nodata_value.
     """
     # Step 1: If input is a filename, load the xarray.Dataset
-    if isinstance(swath_ds, str):
+    if isinstance(swath_dataset, str):
         # Local import to avoid circular import
         from .emit_xarray import emit_xarray
-        ds: xr.Dataset = emit_xarray(swath_ds, ortho=False)
+        ds: xr.Dataset = emit_xarray(swath_dataset, ortho=False)
     else:
-        ds: xr.Dataset = swath_ds
+        ds: xr.Dataset = swath_dataset
 
     # 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

EMITL2ARFL/granule.py

@@ -4,14 +4,20 @@
 
 import numpy as np
 
-import netCDF4
+
+
 
 from rasters import Raster, RasterGeometry, RasterGeolocation
+from .read_latitude import read_latitude
+from .read_longitude import read_longitude
+from .read_geolocation import read_geolocation
 
 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
+from .extract_GLT import extract_GLT
+from .GLT import GLT
 
 class EMITL2ARFL:
     def __init__(self, directory: str):
@@ -42,29 +48,18 @@ def uncertainty_filename(self) -> str:
 
     @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
+        return read_latitude(self.reflectance_filename)
 
     @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
+        return read_longitude(self.reflectance_filename)
 
     @property
     def geolocation(self) -> RasterGeolocation:
-        return RasterGeolocation(
-            x=self.lon,
-            y=self.lat
-        )
+        return read_geolocation(self.reflectance_filename)
 
-    def GLT(self) -> Raster:
-        return Raster(extract_GLT_array(swath_ds=self.reflectance_filename), geometry=self.geolocation)
+    def GLT(self) -> GLT:
+        return extract_GLT(self.reflectance_filename)
 
     def quality_mask(self, quality_bands: List[int] = QUALITY_BANDS) -> np.ndarray:
         qmask = quality_mask(

EMITL2ARFL/ortho_xr.py

@@ -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_array(swath_ds=swath_ds)
+    GLT_array = extract_GLT_array(swath_dataset=swath_ds)
 
     # List Variables
     var_list = list(swath_ds.data_vars)

EMITL2ARFL/read_geolocation.py

@@ -0,0 +1,27 @@
+from rasters import RasterGeolocation
+from .read_latitude import read_latitude
+from .read_longitude import read_longitude
+
+def read_geolocation(filename: str) -> RasterGeolocation:
+    """
+    Reads the latitude and longitude arrays from a NetCDF reflectance file and constructs a RasterGeolocation object.
+
+    This function uses helper functions to extract the latitude and longitude arrays from the specified NetCDF file.
+    It then creates and returns a RasterGeolocation object, which is used to represent the geospatial coordinates
+    associated with the raster data. This is useful for mapping pixel locations to geographic coordinates.
+
+    Args:
+        filename (str): Path to the NetCDF reflectance file containing the 'lat' and 'lon' arrays in the 'location' group.
+
+    Returns:
+        RasterGeolocation: An object containing the longitude (x) and latitude (y) arrays for georeferencing.
+    """
+    # Read the latitude array from the NetCDF file using the helper function
+    lat = read_latitude(filename)
+    # Read the longitude array from the NetCDF file using the helper function
+    lon = read_longitude(filename)
+    # Create a RasterGeolocation object using the longitude and latitude arrays
+    geolocation = RasterGeolocation(x=lon, y=lat)
+
+    # Return the RasterGeolocation object for use in geospatial operations
+    return geolocation

EMITL2ARFL/read_latitude.py

@@ -0,0 +1,8 @@
+import netCDF4
+import numpy as np
+
+def read_latitude(reflectance_filename: str) -> np.ndarray:
+    """Read the `lat` array from the `location` group in the reflectance NetCDF file."""
+    with netCDF4.Dataset(reflectance_filename, "r") as ds:
+        lat = ds.groups["location"].variables["lat"][:]
+    return lat

EMITL2ARFL/read_longitude.py

@@ -0,0 +1,8 @@
+import netCDF4
+import numpy as np
+
+def read_longitude(reflectance_filename: str) -> np.ndarray:
+    """Read the `lon` array from the `location` group in the reflectance NetCDF file."""
+    with netCDF4.Dataset(reflectance_filename, "r") as ds:
+        lon = ds.groups["location"].variables["lon"][:]
+    return lon

pyproject.toml

@@ -20,7 +20,7 @@ dependencies = [
     "geopandas",
     "netCDF4",
     "python-dateutil",
-    "rasters>=1.5.3",
+    "rasters>=1.8.0",
     "rioxarray",
     "spectral",
     "xarray"