update forge-py to have one algorithm for all collections

CHANGELOG.md

@@ -8,8 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-- Added cowvr strategy
-- Update default footprint algorithm
+- Update default footprint algorithm, now there is only one algorithm for forge-py
 - Poetry update python libraries
 ### Deprecated
 ### Removed

podaac/forge_py/cli.py

@@ -7,6 +7,7 @@
 import json
 from datetime import datetime, timezone
 import xarray as xr
+import numpy as np
 
 from podaac.forge_py.args import parse_args
 from podaac.forge_py.file_util import make_absolute
@@ -71,13 +72,14 @@ def main(args=None):
     group = read_config.get('footprint', {}).get('group')
     cutoff_lat = read_config.get('footprint', {}).get('cutoff_lat', None)
     smooth_poles = read_config.get('footprint', {}).get('smooth_poles', None)
+    fill_value = read_config.get('footprint', {}).get('fill_value', np.nan)
 
     # 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,
-                                                      cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, strategy=strategy)
+                                                      cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, strategy=strategy, fill_value=fill_value)
 
     if args.output_file:
         with open(args.output_file, "w") as json_file:

podaac/forge_py/forge.py

@@ -6,83 +6,8 @@
 import shapely
 
 
-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.
-
-    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()
-
-    # Outlying data near the poles. As a quick fix, remove all data near the poles, at latitudes higher than
-    # 87 degrees. This quick fix has impact on footprint shape.
-    i_lolats = np.where(abs(y) < 86)
-    x = x[i_lolats]
-    y = y[i_lolats]
-
-    # 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)]
-
-    # Fit with alphashape
-    xy = np.array(list(zip(x_thin, y_thin)))  # Reshape coords to use with alphashape
-    alpha_shape = alphashape.alphashape(xy, alpha=alpha)
-
-    # Because of the thinning processes, the pole-edges of the footprint are jagged rather than
-    # flat, quick fix this by making all latitude points above 85 degrees a constant value:
-    fp_lon, fp_lat = alpha_shape.exterior.coords.xy
-    fp_lat = np.array(fp_lat)
-    fp_lat[np.where(fp_lat > 82)] = 88
-    fp_lat[np.where(fp_lat < -82)] = -88
-    footprint = Polygon(list(zip(fp_lon, np.asarray(fp_lat, dtype=np.float64))))
-
-    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, cutoff_lat=None, smooth_poles=None):
+def generate_footprint(lon, lat, thinning_fac=30, alpha=0.05, is360=False, simplify=0.1,
+                       strategy=None, cutoff_lat=None, smooth_poles=None, fill_value=np.nan):  # pylint: disable=unused-argument
     """
     Generates footprint by calling different footprint strategies
 
@@ -98,24 +23,27 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
         simplify polygon factor
     strategy:
         What footprint strategy to use
+    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".
+    fill_value: (optional) float
+        Fill value in the latitude, longitude arrays. Default = np.nan; the default
+        will work even if the data have no NAN's. Future functionality will accommodate
+        multiple possible fill values.
     """
 
-    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, 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
-    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)
+    thinning = {'method': 'standard', 'value': thinning_fac}
+    alpha_shape = fit_footprint(lon_array, lat, alpha=alpha, thinning=thinning, cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, fill_value=fill_value)
     alpha_shape = alpha_shape.simplify(simplify)
 
     # If the polygon is not valid, attempt to fix self-intersections
@@ -129,7 +57,7 @@ def generate_footprint(lon, lat, thinning_fac=30, alpha=0.035, is360=False, simp
 def fit_footprint(
         lon, lat, alpha=0.05,
         thinning=None, cutoff_lat=None,
-        smooth_poles=None,
+        smooth_poles=None, fill_value=np.nan,
         return_xythin=False):
     """
     Fits instrument coverage footprint for level 2 data set. Output is a polygon object for
@@ -157,14 +85,24 @@ def fit_footprint(
         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".
+    fill_value: (optional) float
+        Fill value in the latitude, longitude arrays. Default = np.nan; the default
+        will work even if the data have no NAN's. Future functionality will accommodate
+        multiple possible fill values.
     return_xythin: bool, default = False
         If True, returns the thinned out latitude, longitude arrays along with the
         footprint.
     """
 
-    x, y = np.array(lon).flatten(), np.array(lat).flatten()
-    mask = ~np.isnan(x) & ~np.isnan(y)
-    x, y = x[mask], y[mask]
+    # Prep arrays and remove missing values:
+    x = np.array(lon).flatten()
+    y = np.array(lat).flatten()
+    if fill_value is np.nan:
+        inan = np.isnan(x*y)
+    else:
+        inan = (x == fill_value) | (y == fill_value)
+    x = x[~inan]
+    y = y[~inan]
 
     # Optional thinning (typically helps alphashape fit faster):
     if thinning is not None:

podaac/lambda_handler/lambda_handler.py

@@ -8,6 +8,7 @@
 import requests
 
 import botocore
+import numpy as np
 import xarray as xr
 from cumulus_logger import CumulusLogger
 from cumulus_process import Process, s3
@@ -191,13 +192,14 @@ def footprint_generate(self, file_, config_file, granule_id):
         group = read_config.get('footprint', {}).get('group')
         cutoff_lat = read_config.get('footprint', {}).get('cutoff_lat', None)
         smooth_poles = read_config.get('footprint', {}).get('smooth_poles', None)
+        fill_value = read_config.get('footprint', {}).get('fill_value', np.nan)
 
         # 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,
-                                                          cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, strategy=strategy)
+                                                          cutoff_lat=cutoff_lat, smooth_poles=smooth_poles, strategy=strategy, fill_value=fill_value)
 
         wkt_json = {
             "FOOTPRINT": wkt_representation,