update footprint

.github/workflows/build.yml

@@ -203,8 +203,27 @@ jobs:
           poetry run flake8 podaac
           poetry run pytest --junitxml=build/reports/pytest.xml --cov=podaac/ --cov-report=html -m "not aws and not integration" tests/
 
-      ## TODO: Find out where the test report goes
-
+      - name: Run Snyk as a blocking step
+        uses: snyk/actions/python@master
+        env:
+            SNYK_TOKEN: ${{ secrets.SNYK_TOKEN }}
+        with:
+          command: test
+          args: >
+              --org=${{ secrets.SNYK_ORG_ID }}
+              --project-name=${{ github.repository }}
+              --severity-threshold=high
+              --fail-on=all
+
+      - name: Run Snyk on Python
+        uses: snyk/actions/python@master
+        env:
+          SNYK_TOKEN: ${{ secrets.SNYK_TOKEN }}
+        with:
+          command: monitor
+          args: >
+            --org=${{ secrets.SNYK_ORG_ID }}
+            --project-name=${{ github.repository }}
 
       #########################################################################
       # Build

podaac/forge_py/forge.py

@@ -1,32 +1,9 @@
 """Python footprint generator"""
 import numpy as np
 import alphashape
-from shapely.geometry import Polygon
+from shapely.geometry import Polygon, MultiPolygon
 from shapely.wkt import dumps
-
-
-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.
-    """
-
-    # lat, lon need to be 1D:
-    x = np.array(lon).flatten()
-    y = np.array(lat).flatten()
-
-    # Thinning out the number of data points helps alphashape fit faster
-    x_thin = x[np.arange(0, len(x), thinning_fac)]
-    y_thin = y[np.arange(0, len(y), thinning_fac)]
-
-    xy = np.array(list(zip(x_thin, y_thin)))  # Reshape coords to use with alphashape
-    alpha_shape = alphashape.alphashape(xy, alpha=alpha)
-
-    return alpha_shape
+import shapely
 
 
 def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035):
@@ -104,7 +81,8 @@ def cowvr_footprint(lon, lat, thinning_fac=200, alpha=0.03):
     return alpha_shape
 
 
-def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simplify=0.1, strategy=None):
+def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simplify=0.1,
+                       strategy=None, cutoff_lat=None, smooth_poles=None):
     """
     Generates footprint by calling different footprint strategies
 
@@ -127,13 +105,17 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
         "cowvr": cowvr_footprint
     }
     # Get the function corresponding to the strategy, or the default function
-    selected_function = strategy_functions.get(strategy, fit_footprint)
+    selected_function = strategy_functions.get(strategy, None)
     # 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)
+    if selected_function is not None:
+        alpha_shape = selected_function(lon_array, lat, thinning_fac=thinning_fac, alpha=alpha)
+    else:
+        thinning = {'method': 'standard', 'value': thinning_fac}
+        alpha_shape = fit_footprint(lon_array, lat, thinning=thinning, cutoff_lat=cutoff_lat, smooth_poles=smooth_poles)
     alpha_shape = alpha_shape.simplify(simplify)
 
     # If the polygon is not valid, attempt to fix self-intersections
@@ -142,3 +124,92 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
 
     wkt_alphashape = dumps(alpha_shape)
     return wkt_alphashape
+
+
+def fit_footprint(
+        lon, lat, alpha=0.05,
+        thinning=None, cutoff_lat=None,
+        smooth_poles=None,
+        return_xythin=False):
+    """
+    Fits instrument coverage footprint for level 2 data set. Output is a polygon object for
+    the indices of the footprint outline. Uses the alphashape package for the fit,
+    which returns a shapely.geometry.polygon.Polygon or
+    shapely.geometry.multipolygon.MultiPolygon object.
+
+    lon, lat: list/array-like's
+        Latitudes and longitudes of instrument coverage. Should be the same shape and size.
+    alpha: float
+        The alpha parameter passed to alphashape. Typical values that work for
+        L2 footprinting are in the range 0.02 - 0.06.
+    thinning: (optional) dictionary
+        Optional method for removing some of the data points in the lon, lat arrays. It is
+        often handy because thinning out the data makes the fit faster. Dict keys are
+        "method" and "value". If "method" is set to "standard", then e.g. a "value" of
+        100 will thin out the arrays to every 100th point; in this case "value" should be
+        an int.
+    cutoff_lat: (optional) float, default = None
+        If specified, latitudes higher than this threshold value will be
+        removed before the fit is performed. This works in both the north and
+        south direction, e.g. for a value of x, both points north of x and south
+        of -x will be removed.
+    smooth_poles: (optional) 2-tuple of floats, default = None
+        If specified, the first element gives the threshold latitude above which
+        any footprint indicies will have their latitudes set to the value of the
+        second element in "smooth_poles".
+    return_xythin: bool, default = False
+        If True, returns the thinned out latitude, longitude arrays along with the
+        footprint.
+    """
+
+    # Prep arrays:
+    x = np.array(lon).flatten()
+    y = np.array(lat).flatten()
+    inan = np.isnan(x*y)
+    x = x[~inan]
+    y = y[~inan]
+
+    # Optional thinning (typically helps alphashape fit faster):
+    if thinning is not None:
+        if thinning["method"] == "standard":
+            x_thin = x[np.arange(0, len(x), thinning["value"])]
+            y_thin = y[np.arange(0, len(y), thinning["value"])]
+    else:
+        x_thin = x
+        y_thin = y
+
+    # Optional removal of "outlying" data near the poles. Removes data at latitudes
+    # higher than the specified value. This will have an impact on footprint shape.
+    if cutoff_lat is not None:
+        i_lolats = np.where(abs(y_thin) < cutoff_lat)
+        x_thin = x_thin[i_lolats]
+        y_thin = y_thin[i_lolats]
+
+    # Fit with alphashape
+    xy = np.array(list(zip(x_thin, y_thin)))  # Reshape coords to use with alphashape
+    footprint = alpha_shape = alphashape.alphashape(xy, alpha=alpha)
+
+    # Optional pole smoothing: if the data was thinned, the fitted footprint may
+    # have jagged pole-edges. This can be optionally smoothed by making all
+    # latitudes higher than some threshold a constant value:
+    def pole_smoother(fp_lon, fp_lat, lat_thresh, lat_target):
+        """
+        Takes longitude, latitude array-likes from a single Polygon representing a footprint.
+        """
+        fp_lat = np.array(fp_lat)
+        fp_lat[np.where(fp_lat > lat_thresh)] = lat_target
+        fp_lat[np.where(fp_lat < -lat_thresh)] = -lat_target
+        return Polygon(list(zip(fp_lon, np.asarray(fp_lat, dtype=np.float64))))
+
+    if smooth_poles is not None:
+        if isinstance(alpha_shape, shapely.geometry.polygon.Polygon):
+            footprint = pole_smoother(*alpha_shape.exterior.coords.xy, *smooth_poles)
+        elif isinstance(alpha_shape, shapely.geometry.multipolygon.MultiPolygon):
+            footprint = MultiPolygon([
+                pole_smoother(*p.exterior.coords.xy, *smooth_poles)
+                for p in alpha_shape.geoms
+            ])
+
+    if return_xythin:
+        return footprint, x_thin, y_thin
+    return footprint

podaac/lambda_handler/lambda_handler.py

@@ -189,12 +189,15 @@ def footprint_generate(self, file_, config_file, granule_id):
         strategy = read_config.get('footprint', {}).get('strategy', None)
         simplify = read_config.get('footprint', {}).get('simplify', 0.1)
         group = read_config.get('footprint', {}).get('group')
+        cutoff_lat = read_config.get('cutoff_lat', None)
+        smooth_poles = read_config.get('smooth_poles', None)
 
         # Generate footprint
         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)
+            wkt_representation = forge.generate_footprint(lon_data, lat_data, thinning_fac=thinning_fac, alpha=alpha, is360=is360, simplify=simplify,
+                                                          cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, strategy=strategy)
 
         wkt_json = {
             "FOOTPRINT": wkt_representation,