Merge pull request #11 from podaac/feature/cowvr

added a cowvr collection footprint generator

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- Added cowvr strategy
 ### Deprecated
 ### Removed
 ### Fixed

podaac/forge_py/forge.py

@@ -5,24 +5,18 @@
 from shapely.wkt import dumps
 
 
-def fit_footprint(lon, lat, thinning_fac=100, alpha=0.05, is360=False):
+def fit_footprint(lon, lat, thinning_fac=100, alpha=0.05):
     """
     lon, lon: list/array-like
         Latitudes and longitudes.
     thinning_fac: int
         Factor to thin out data by (makes alphashape fit faster).
     alpha: float
         The alpha parameter passed to alphashape.
-    is360: bool
-        Tell us if the logitude data is between 0-360
     """
 
-    lon_array = lon
-    if is360:
-        lon_array = ((lon + 180) % 360.0) - 180
-
     # lat, lon need to be 1D:
-    x = np.array(lon_array).flatten()
+    x = np.array(lon).flatten()
     y = np.array(lat).flatten()
 
     # Thinning out the number of data points helps alphashape fit faster
@@ -35,7 +29,7 @@ def fit_footprint(lon, lat, thinning_fac=100, alpha=0.05, is360=False):
     return alpha_shape
 
 
-def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False):
+def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035):
     """
     Fits footprint g-polygon for level 2 data set SCATSAT1_ESDR_L2_WIND_STRESS_V1.1. Uses the
     alphashape package for the fit, which returns a shapely.geometry.polygon.Polygon object.
@@ -46,16 +40,11 @@ def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False):
         Factor to thin out data by (makes alphashape fit faster).
     alpha: float
         The alpha parameter passed to alphashape.
-    is360: bool
-        Tell us if the logitude data is between 0-360
+
     """
     # lat, lon need to be 1D:
 
-    lon_array = lon
-    if is360:
-        lon_array = ((lon + 180) % 360.0) - 180
-
-    x = np.array(lon_array).flatten()
+    x = np.array(lon).flatten()
     y = np.array(lat).flatten()
 
     # Outlying data near the poles. As a quick fix, remove all data near the poles, at latitudes higher than
@@ -83,6 +72,38 @@ def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False):
     return footprint
 
 
+def cowvr_footprint(lon, lat, thinning_fac=200, alpha=0.03):
+    """
+    Uses alphashape to get the footprint from a COWVR EDR or TSDR file using the lat, lon data.
+    Returns: (1) the alpha shape object (contains a shapely object with the footprint coords),
+
+    lon, lon: list/array-like
+        Latitudes and longitudes.
+    thinning_fac: int
+        Factor to thin out data by (makes alphashape fit faster).
+    alpha: float
+        The alpha parameter passed to alphashape.
+    """
+
+    # Remove missing values:
+    lon_qc = lon.where(lon > -999999).values
+    lat_qc = lat.where(lat > -999999).values
+
+    ifinite = np.isfinite(lon_qc * lat_qc)
+    lon_qc = lon_qc[ifinite]
+    lat_qc = lat_qc[ifinite]
+
+    # Thin out arrays so that alphashape doesn't have to work as hard:
+    lon_thin = lon_qc[np.arange(0, len(lon_qc), thinning_fac)]
+    lat_thin = lat_qc[np.arange(0, len(lat_qc), thinning_fac)]
+
+    # Fit the footprint to a polygon:
+    xy = np.array(list(zip(lon_thin, lat_thin)))  # Reshape coords to use with alphashape
+    alpha_shape = alphashape.alphashape(xy, alpha=alpha)
+
+    return alpha_shape
+
+
 def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simplify=0.1, strategy=None):
     """
     Generates footprint by calling different footprint strategies
@@ -101,10 +122,18 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
         What footprint strategy to use
     """
 
-    if strategy == "scatsat":
-        alpha_shape = scatsat_footprint(lon, lat, thinning_fac=thinning_fac, alpha=alpha, is360=is360)
-    else:
-        alpha_shape = fit_footprint(lon, lat, thinning_fac=thinning_fac, alpha=alpha, is360=is360)
+    strategy_functions = {
+        "scatsat": scatsat_footprint,
+        "cowvr": cowvr_footprint
+    }
+    # Get the function corresponding to the strategy, or the default function
+    selected_function = strategy_functions.get(strategy, fit_footprint)
+    # Transform lon array if it is 360
+    lon_array = lon
+    if is360:
+        lon_array = ((lon + 180) % 360.0) - 180
+    # Call the selected function with the provided arguments
+    alpha_shape = selected_function(lon_array, lat, thinning_fac=thinning_fac, alpha=alpha)
     alpha_shape = alpha_shape.simplify(simplify)
 
     # If the polygon is not valid, attempt to fix self-intersections

podaac/lambda_handler/lambda_handler.py

@@ -188,9 +188,10 @@ def footprint_generate(self, file_, config_file, granule_id):
         alpha = read_config.get('footprint', {}).get('alpha', 0.05)
         strategy = read_config.get('footprint', {}).get('strategy', None)
         simplify = read_config.get('footprint', {}).get('simplify', 0.1)
+        group = read_config.get('footprint', {}).get('group')
 
         # Generate footprint
-        with xr.open_dataset(local_file, decode_times=False) as ds:
+        with xr.open_dataset(local_file, group=group, decode_times=False) as ds:
             lon_data = ds[longitude_var]
             lat_data = ds[latitude_var]
             wkt_representation = forge.generate_footprint(lon_data, lat_data, thinning_fac=thinning_fac, alpha=alpha, is360=is360, simplify=simplify, strategy=strategy)