update footprint code

.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,35 +1,13 @@
 """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
+import shapely
 
 
-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
-
-
-def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035):
+def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035, poleoutliers_thresh_lat=None, smooth_poles=None):
+    # pylint: disable=unused-argument
     """
     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.
@@ -72,7 +50,8 @@ def scatsat_footprint(lon, lat, thinning_fac=30, alpha=0.035):
     return footprint
 
 
-def cowvr_footprint(lon, lat, thinning_fac=200, alpha=0.03):
+def cowvr_footprint(lon, lat, thinning_fac=200, alpha=0.03, poleoutliers_thresh_lat=None, smooth_poles=None):
+    # pylint: disable=unused-argument
     """
     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),
@@ -104,7 +83,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, poleoutliers_thresh_lat=None, smooth_poles=None):
     """
     Generates footprint by calling different footprint strategies
 
@@ -133,7 +113,8 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
     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 = selected_function(lon_array, lat, thinning_fac=thinning_fac, alpha=alpha,
+                                    poleoutliers_thresh_lat=poleoutliers_thresh_lat, smooth_poles=smooth_poles)
     alpha_shape = alpha_shape.simplify(simplify)
 
     # If the polygon is not valid, attempt to fix self-intersections
@@ -142,3 +123,83 @@ 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, thinning_fac=100, alpha=0.05,
+                  poleoutliers_thresh_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.
+    thinning_fac: int
+        Factor to thin out data by, e.g. a value of 100 will take every 100th point in the
+        lon, lat arrays. Thinning the data out typically makes alphashape fit faster.
+    alpha: float
+        The alpha parameter passed to alphashape. Typical values that work for
+        L2 footprinting are in the range 0.02 - 0.06.
+    poleoutliers_thresh_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]
+    # 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)]
+
+    # 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 poleoutliers_thresh_lat is not None:
+        i_lolats = np.where(abs(y_thin) < poleoutliers_thresh_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)
+
+    # Because of the thinning processes, the pole-edges of the footprint may be jagged
+    # rather than flat. This is optionally addressed by making all latitudes higher
+    # than some thresholda 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
+            ])
+
+    wkt_alphashape = dumps(footprint)
+
+    if return_xythin:
+        return wkt_alphashape, x_thin, y_thin
+    return wkt_alphashape

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')
+        poleoutliers_thresh_lat = read_config.get('poleoutliers_thresh_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,
+                                                          poleoutliers_thresh_lat=poleoutliers_thresh_lat, smooth_poles=smooth_poles, strategy=strategy)
 
         wkt_json = {
             "FOOTPRINT": wkt_representation,