diff --git a/podaac/forge_py/cli.py b/podaac/forge_py/cli.py index de33591..8957bc1 100644 --- a/podaac/forge_py/cli.py +++ b/podaac/forge_py/cli.py @@ -8,7 +8,6 @@ import os 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 @@ -59,46 +58,14 @@ def main(args=None): config_file = args.config local_file = args.granule - with open(config_file) as config_f: - read_config = json.load(config_f) - - longitude_var = read_config.get('lonVar') - latitude_var = read_config.get('latVar') - is360 = read_config.get('is360', False) - - # Get footprint configuration - footprint_config = read_config.get('footprint', {}) - - thinning_fac = footprint_config.get('thinning_fac', 100) - thinning_method = footprint_config.get('thinning_method', 'standard') - alpha = footprint_config.get('alpha', 0.05) - strategy = footprint_config.get('strategy', None) - simplify = footprint_config.get('simplify', 0.1) - group = footprint_config.get('group') - cutoff_lat = footprint_config.get('cutoff_lat', None) - smooth_poles = footprint_config.get('smooth_poles', None) - fill_value = footprint_config.get('fill_value', np.nan) - width = footprint_config.get('width', 3600) - height = footprint_config.get('height', 1800) - - # 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] + strategy, footprint_params = forge.load_footprint_config(config_file) + footprint_params["path"] = os.getcwd() + with xr.open_dataset(local_file, group=footprint_params.get('group'), decode_times=False) as ds: + lon_data = ds[footprint_params['longitude_var']] + lat_data = ds[footprint_params['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, - fill_value=fill_value, - thinning_method=thinning_method, - width=width, - height=height, - path=os.getcwd() + lon_data, lat_data, strategy=strategy, **footprint_params ) if args.output_file: diff --git a/podaac/forge_py/forge.py b/podaac/forge_py/forge.py index 4104a33..2b25a1d 100644 --- a/podaac/forge_py/forge.py +++ b/podaac/forge_py/forge.py @@ -1,4 +1,5 @@ """Python footprint generator""" +import json import numpy as np import alphashape from shapely.geometry import Polygon, MultiPolygon @@ -8,6 +9,40 @@ from podaac.forge_py import open_cv_footprint +def load_footprint_config(config_file): + """ + Load and process the footprint configuration from a JSON file. + + Parameters: + ---------- + config_file : str + Path to the JSON configuration file. + + Returns: + ------- + tuple + (strategy, config), where `strategy` is the footprint strategy to use, + and `config` is a dictionary of parameters specific to that strategy. + """ + with open(config_file) as config_f: + read_config = json.load(config_f) + + # Select the specified strategy and its parameters + footprint_config = read_config.get('footprint', {}) + strategy = footprint_config.get('strategy', 'alpha_shape') # Default to 'alpha_shape' + strategy_params = footprint_config.get(strategy, {}) # Get params for chosen strategy + + # Include general options like is360 if needed outside strategies + common_params = { + 'is360': read_config.get('is360', False), + 'longitude_var': read_config.get('lonVar'), + 'latitude_var': read_config.get('latVar') + } + + # Merge strategy parameters with any additional common parameters + return strategy, {**common_params, **strategy_params} + + def thinning_bin_avg(x, y, rx, ry): """ One of the options for fit_footprint(). @@ -30,64 +65,54 @@ def thinning_bin_avg(x, y, rx, ry): return xy_thinned["x"].values, xy_thinned["y"].values -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, - thinning_method="standard", width=3600, height=1800, path=None): +def generate_footprint(lon, lat, strategy=None, is360=False, path=None, **kwargs): """ - Generates footprint by calling different footprint strategies - - 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 - simplify: - 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. + Generates a geographic footprint using a specified strategy. + + Parameters: + ---------- + lon, lat : list or array-like + Latitude and longitude values. + strategy : str, optional + The footprint strategy to use ('open_cv' or 'alpha_shape'). + is360 : bool, default=False + If True, adjusts longitude values from 0-360 to -180-180 range. + path : str, optional + File path for saving output if the strategy requires it. + **kwargs : dict, optional + Additional parameters to be passed to the chosen strategy function. + + Returns: + ------- + str + The footprint as a WKT (Well-Known Text) string. """ - - # Transform lon array if it is 360 - lon_array = lon + # Adjust longitude for 0-360 to -180-180 range if needed if is360: - lon_array = ((lon + 180) % 360.0) - 180 - thinning = {'method': thinning_method, 'value': thinning_fac} + lon = ((lon + 180) % 360.0) - 180 + + # Dispatch to the correct footprint strategy based on `strategy` if strategy == "open_cv": - alpha_shape_wkt = open_cv_footprint.footprint_open_cv(lon, lat, width=width, height=height, path=path) - return alpha_shape_wkt + return open_cv_footprint.footprint_open_cv(lon, lat, path=path, **kwargs) + + # Default to alpha_shape strategy if no strategy is specified + alpha_shape = fit_footprint(lon, lat, **kwargs) - 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) + # Simplify and validate the alpha shape if requested + if 'simplify' in kwargs: + alpha_shape = alpha_shape.simplify(kwargs['simplify']) - # If the polygon is not valid, attempt to fix self-intersections if not alpha_shape.is_valid: alpha_shape = alpha_shape.buffer(0) - wkt_alphashape = dumps(alpha_shape) - return wkt_alphashape + return dumps(alpha_shape) def fit_footprint( lon, lat, alpha=0.05, thinning=None, cutoff_lat=None, smooth_poles=None, fill_value=np.nan, - return_xythin=False): + return_xythin=False, **kwargs): """ 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, @@ -174,6 +199,7 @@ def pole_smoother(fp_lon, fp_lat, lat_thresh, lat_target): # Return the smoothed polygon return Polygon(zip(fp_lon, fp_lat)) + smooth_poles = None if smooth_poles is not None: if isinstance(alpha_shape, shapely.geometry.polygon.Polygon): footprint = pole_smoother(*alpha_shape.exterior.coords.xy, *smooth_poles) diff --git a/podaac/forge_py/open_cv_footprint.py b/podaac/forge_py/open_cv_footprint.py index e0c89dd..bfc436c 100644 --- a/podaac/forge_py/open_cv_footprint.py +++ b/podaac/forge_py/open_cv_footprint.py @@ -396,7 +396,7 @@ def round_coords(coords, precision): raise ValueError("Unsupported geometry type") -def footprint_open_cv(lon, lat, width=3600, height=1800, path=None, threshold_value=185): +def footprint_open_cv(lon, lat, width=3600, height=1800, path=None, threshold_value=185, **kwargs): """ Main pipeline for processing geographic coordinates to create a footprint polygon using image processing techniques. @@ -447,7 +447,6 @@ def footprint_open_cv(lon, lat, width=3600, height=1800, path=None, threshold_va simplified_polygon = simplify_polygon(polygon_structure) reduced_precision = reduce_precision(simplified_polygon) counter_clockwise = ensure_counter_clockwise(reduced_precision) - print(counter_clockwise.wkt) return counter_clockwise.wkt raise Exception("No valid polygons found.") diff --git a/podaac/lambda_handler/lambda_handler.py b/podaac/lambda_handler/lambda_handler.py index 3e02397..6fc64ab 100644 --- a/podaac/lambda_handler/lambda_handler.py +++ b/podaac/lambda_handler/lambda_handler.py @@ -8,7 +8,6 @@ import requests import botocore -import numpy as np import xarray as xr from cumulus_logger import CumulusLogger from cumulus_process import Process, s3 @@ -178,46 +177,15 @@ def footprint_generate(self, file_, config_file, granule_id): self.logger.error("Error downloading granule from s3: {}".format(ex), exc_info=True) raise ex - with open(config_file) as config_f: - read_config = json.load(config_f) - - longitude_var = read_config.get('lonVar') - latitude_var = read_config.get('latVar') - is360 = read_config.get('is360', False) - - # Get footprint configuration - footprint_config = read_config.get('footprint', {}) - - thinning_fac = footprint_config.get('thinning_fac', 100) - thinning_method = footprint_config.get('thinning_method', 'standard') - alpha = footprint_config.get('alpha', 0.05) - strategy = footprint_config.get('strategy', None) - simplify = footprint_config.get('simplify', 0.1) - group = footprint_config.get('group') - cutoff_lat = footprint_config.get('cutoff_lat', None) - smooth_poles = footprint_config.get('smooth_poles', None) - fill_value = footprint_config.get('fill_value', np.nan) - width = footprint_config.get('width', 3600) - height = footprint_config.get('height', 1800) - - # 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] + strategy, footprint_params = forge.load_footprint_config(config_file) + footprint_params["path"] = self.path + print(footprint_params) + with xr.open_dataset(local_file, group=footprint_params.get('group'), decode_times=False) as ds: + lon_data = ds[footprint_params['longitude_var']] + lat_data = ds[footprint_params['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, - fill_value=fill_value, - thinning_method=thinning_method, - width=width, - height=height, - path=self.path + lon_data, lat_data, strategy=strategy, **footprint_params ) wkt_json = { diff --git a/tests/input/SCATSAT1_ESDR_L2_WIND_STRESS_V1.1.cfg b/tests/input/SCATSAT1_ESDR_L2_WIND_STRESS_V1.1.cfg index 4bdcda0..2bec2a0 100644 --- a/tests/input/SCATSAT1_ESDR_L2_WIND_STRESS_V1.1.cfg +++ b/tests/input/SCATSAT1_ESDR_L2_WIND_STRESS_V1.1.cfg @@ -4,61 +4,22 @@ "lonVar": "lon", "timeVar": "time", "is360": true, - "tiles": { - "steps": [ - 30, - 14 - ] + "footprint":{ + "strategy":"alpha_shape", + "alpha_shape":{ + "thinning":{ + "method":"standard", + "value":30 + }, + "alpha":0.035, + "simplify":0.1, + "cutoff_lat":88, + "smooth_poles": [82,90] + }, + "open_cv":{ + "width":3600, + "height":1800 + } }, - "footprint": { - "thinning_fac": 30, - "alpha": 0.035, - "simplify": 0.1, - "cutoff_lat": 88, - "smooth_poles": [82,90] - }, - "footprinter": "forge-py", - "imgVariables": [], - "image": { - "ppd": 16, - "res": 8 - }, - "variables": [ - "time", - "lat", - "lon", - "flags", - "quality_indicator", - "nudge_wind_speed", - "nudge_wind_direction", - "cross_track_wind_speed_bias", - "rain_speed_bias", - "gmf_sst", - "distance_from_coast", - "en_wind_speed", - "en_wind_direction", - "en_wind_u", - "en_wind_v", - "en_wind_speed_uncorrected", - "en_wind_direction_uncorrected", - "en_wind_speed_error", - "en_wind_direction_error", - "en_wind_u_error", - "en_wind_v_error", - "wind_stress_magnitude", - "wind_stress_direction", - "wind_stress_u", - "wind_stress_v", - "wind_stress_magnitude_error", - "wind_stress_u_error", - "wind_stress_v_error", - "real_wind_speed", - "real_wind_direction", - "real_wind_u", - "real_wind_v", - "real_wind_speed_error", - "real_wind_direction_error", - "real_wind_u_error", - "real_wind_v_error" - ] + "footprinter": "forge-py" } \ No newline at end of file diff --git a/tests/test_footprint_generator.py b/tests/test_footprint_generator.py index cf3e15c..f317cb1 100644 --- a/tests/test_footprint_generator.py +++ b/tests/test_footprint_generator.py @@ -175,7 +175,7 @@ def test_lambda_handler_cumulus(mocked_get): input_dir = f'{test_dir}/input' config_dir = f'{test_dir}/configs' nc_file = f'{input_dir}/measures_esdr_scatsat_l2_wind_stress_23433_v1.1_s20210228-054653-e20210228-072612.nc' - cfg_file = f'{config_dir}/PODAAC-CYGNS-C2H10.cfg' + config_file = f'{input_dir}/SCATSAT1_ESDR_L2_WIND_STRESS_V1.1.cfg' with open(nc_file, 'rb') as data: aws_s3.Bucket(bucket).put_object(Key='test_folder/test_granule.nc', Body=data) @@ -191,7 +191,7 @@ def test_lambda_handler_cumulus(mocked_get): aws_s3.create_bucket(Bucket='internal-bucket') - with open(cfg_file, 'rb') as data: + with open(config_file, 'rb') as data: s3_client.put_object(Bucket='internal-bucket', Key='dataset-config/JASON-1_L2_OST_GPN_E.cfg', Body=data) @@ -263,26 +263,19 @@ def test_forge_py(): with open(result_file, "r") as file: polygon_shape = file.read() - with open(config_file) as config_f: - read_config = json.load(config_f) - - longitude_var = read_config.get('lonVar') - latitude_var = read_config.get('latVar') - is360 = read_config.get('is360', False) - thinning_fac = read_config.get('footprint', {}).get('thinning_fac', 100) - 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) + strategy, footprint_params = forge.load_footprint_config(config_file) # Generate footprint with xr.open_dataset(nc_file, decode_times=False) as ds: - lon_data = ds[longitude_var] - lat_data = ds[latitude_var] - wkt_alphashape = forge.generate_footprint(lon_data, lat_data, thinning_fac=thinning_fac, alpha=alpha, is360=is360, simplify=simplify, strategy=strategy) + lon_data = ds[footprint_params['longitude_var']] + lat_data = ds[footprint_params['latitude_var']] + wkt_alphashape = forge.generate_footprint( + lon_data, lat_data, strategy=strategy, **footprint_params + ) + print(wkt_alphashape) assert compare_shapes_similarity(wkt_alphashape, polygon_shape) - def test_forge_py_open_cv(): test_dir = os.path.dirname(os.path.realpath(__file__)) @@ -294,27 +287,14 @@ def test_forge_py_open_cv(): with open(result_file, "r") as file: polygon_shape = file.read() - with open(config_file) as config_f: - read_config = json.load(config_f) - - longitude_var = read_config.get('lonVar') - latitude_var = read_config.get('latVar') - is360 = read_config.get('is360', False) - thinning_fac = read_config.get('footprint', {}).get('thinning_fac', 100) - alpha = read_config.get('footprint', {}).get('alpha', 0.05) + strategy, footprint_params = forge.load_footprint_config(config_file) + footprint_params["path"] = test_dir strategy = "open_cv" - simplify = read_config.get('footprint', {}).get('simplify', 0.1) - width = read_config.get('footprint', {}).get('width', 3600) - height = read_config.get('footprint', {}).get('height', 1800) - path = test_dir - - # Generate footprint with xr.open_dataset(nc_file, decode_times=False) as ds: - lon_data = ds[longitude_var] - lat_data = ds[latitude_var] - wkt_alphashape = forge.generate_footprint(lon_data, lat_data, thinning_fac=thinning_fac, alpha=alpha, is360=is360, simplify=simplify, strategy=strategy, path=path) - assert compare_shapes_similarity(wkt_alphashape, polygon_shape) - - - + lon_data = ds[footprint_params['longitude_var']] + lat_data = ds[footprint_params['latitude_var']] + wkt_alphashape = forge.generate_footprint( + lon_data, lat_data, strategy=strategy, **footprint_params + ) + assert compare_shapes_similarity(wkt_alphashape, polygon_shape)