feat: added option to crop tides to the domain of the input data

DEM/interp_ATL14_DEM_ICESat2_ATL06.py

@@ -113,6 +113,15 @@ def interp_ATL14_DEM_ICESat2(INPUT_FILE,
     # convert bounding polygon coordinates to projection
     BX, BY = transformer.transform(bounding_lon, bounding_lat)
     BOUNDS = [BX.min(), BX.max(), BY.min(), BY.max()]
+    # check if bounding polygon is in multiple parts
+    if 'bounding_polygon_lon2' in IS2_atl06_mds['orbit_info']:
+        bounding_lon = IS2_atl06_mds['orbit_info']['bounding_polygon_lon2']
+        bounding_lat = IS2_atl06_mds['orbit_info']['bounding_polygon_lat2']
+        BX, BY = transformer.transform(bounding_lon, bounding_lat)
+        BOUNDS[0] = np.minimum(BOUNDS[0], BX.min())
+        BOUNDS[1] = np.maximum(BOUNDS[1], BX.max())
+        BOUNDS[2] = np.minimum(BOUNDS[2], BY.min())
+        BOUNDS[3] = np.maximum(BOUNDS[3], BY.max())
 
     # read ATL14 DEM model files within spatial bounds
     DEM = gz.io.ATL14(DEM_MODEL, BOUNDS=BOUNDS)

DEM/interp_ATL14_DEM_ICESat2_ATL11.py

@@ -113,6 +113,15 @@ def interp_ATL14_DEM_ICESat2(INPUT_FILE,
     # convert bounding polygon coordinates to projection
     BX, BY = transformer.transform(bounding_lon, bounding_lat)
     BOUNDS = [BX.min(), BX.max(), BY.min(), BY.max()]
+    # check if bounding polygon is in multiple parts
+    if 'bounding_polygon_dim2' in IS2_atl11_mds['orbit_info']:
+        bounding_lon = IS2_atl11_mds['orbit_info']['bounding_polygon_lon2']
+        bounding_lat = IS2_atl11_mds['orbit_info']['bounding_polygon_lat2']
+        BX, BY = transformer.transform(bounding_lon, bounding_lat)
+        BOUNDS[0] = np.minimum(BOUNDS[0], BX.min())
+        BOUNDS[1] = np.maximum(BOUNDS[1], BX.max())
+        BOUNDS[2] = np.minimum(BOUNDS[2], BY.min())
+        BOUNDS[3] = np.maximum(BOUNDS[3], BY.max())
 
     # read ATL14 DEM model files within spatial bounds
     DEM = gz.io.ATL14(DEM_MODEL, BOUNDS=BOUNDS)

GZ/calculate_GZ_ICESat2_ATL06.py

@@ -16,23 +16,6 @@
     -O X, --output-directory X: input/output data directory
     --mean-file X: Mean elevation file to remove from the height data
     -T X, --tide X: Tide model to use in correction
-        CATS0201
-        CATS2008
-        TPXO9-atlas
-        TPXO9-atlas-v2
-        TPXO9-atlas-v3
-        TPXO9-atlas-v4
-        TPXO9-atlas-v5
-        TPXO9.1
-        TPXO8-atlas
-        TPXO7.2
-        AODTM-5
-        AOTIM-5
-        AOTIM-5-2018
-        GOT4.7
-        GOT4.8
-        GOT4.10
-        FES2014
     -R X, --reanalysis X: Reanalysis model to run
         ERA-Interim: http://apps.ecmwf.int/datasets/data/interim-full-moda
         ERA5: http://apps.ecmwf.int/data-catalogues/era5/?class=ea

GZ/calculate_GZ_ICESat2_ATL11.py

@@ -18,22 +18,6 @@
     -O X, --output-directory X: input/output data directory
     --mean-file X: Mean elevation file to remove from the height data
     -T X, --tide X: Tide model to use in correction
-        CATS0201
-        CATS2008
-        TPXO9-atlas
-        TPXO9-atlas-v2
-        TPXO9-atlas-v3
-        TPXO9-atlas-v4
-        TPXO9.1
-        TPXO8-atlas
-        TPXO7.2
-        AODTM-5
-        AOTIM-5
-        AOTIM-5-2018
-        GOT4.7
-        GOT4.8
-        GOT4.10
-        FES2014
     -R X, --reanalysis X: Reanalysis model to run
         ERA-Interim: http://apps.ecmwf.int/datasets/data/interim-full-moda
         ERA5: http://apps.ecmwf.int/data-catalogues/era5/?class=ea

GZ/calculate_GZ_ICESat_GLA12.py

@@ -16,23 +16,6 @@
     -O X, --output-directory X: input/output data directory
     --mean-file X: Mean elevation file to remove from the height data
     -T X, --tide X: Tide model to use in correction
-        CATS0201
-        CATS2008
-        TPXO9-atlas
-        TPXO9-atlas-v2
-        TPXO9-atlas-v3
-        TPXO9-atlas-v4
-        TPXO9-atlas-v5
-        TPXO9.1
-        TPXO8-atlas
-        TPXO7.2
-        AODTM-5
-        AOTIM-5
-        AOTIM-5-2018
-        GOT4.7
-        GOT4.8
-        GOT4.10
-        FES2014
     -R X, --reanalysis X: Reanalysis model to run
         DAC: AVISO dynamic atmospheric correction (DAC) model
         ERA-Interim: http://apps.ecmwf.int/datasets/data/interim-full-moda

doc/environment.yml

@@ -2,7 +2,7 @@ name: gz-docs
 channels:
   - conda-forge
 dependencies:
-  - docutils<0.18
+  - docutils
   - fontconfig
   - freetype
   - graphviz

tides/adjust_tides_ICESat2_ATL11.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 adjust_tides_ICESat2_ATL11.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (07/2024)
 Applies interpolated tidal adjustment scale factors to
     ICESat-2 ATL11 annual land ice height data within
     ice sheet grounding zones
@@ -34,6 +34,7 @@
     io/ATL11.py: reads ICESat-2 annual land ice height data files
 
 UPDATE HISTORY:
+    Updated 07/2024: added option to use JSON format definition files
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 08/2023: create s3 filesystem when using s3 urls as input
@@ -64,6 +65,8 @@
 def adjust_tides_ICESat2_ATL11(adjustment_file, INPUT_FILE,
         OUTPUT_DIRECTORY=None,
         TIDE_MODEL=None,
+        DEFINITION_FILE=None,
+        DEFINITION_FORMAT='ascii',
         VERBOSE=False,
         MODE=0o775
     ):
@@ -73,7 +76,11 @@ def adjust_tides_ICESat2_ATL11(adjustment_file, INPUT_FILE,
     logger = pyTMD.utilities.build_logger('pytmd', level=loglevel)
 
     # get tide model parameters
-    model = pyTMD.io.model(None, verify=False).elevation(TIDE_MODEL)
+    if DEFINITION_FILE is not None:
+        model = pyTMD.io.model(None, verify=False).from_file(DEFINITION_FILE,
+            format=DEFINITION_FORMAT)
+    else:
+        model = pyTMD.io.model(None, verify=False).elevation(TIDE_MODEL)
     # source of tide model
     tide_source = TIDE_MODEL
     tide_reference = model.reference
@@ -715,6 +722,7 @@ def arguments():
     )
     parser.convert_arg_line_to_args = gz.utilities.convert_arg_line_to_args
     # command line parameters
+    group = parser.add_mutually_exclusive_group(required=True)
     # input ICESat-2 annual land ice height files
     parser.add_argument('infile',
         type=pathlib.Path, nargs='+',
@@ -728,10 +736,17 @@ def arguments():
         type=pathlib.Path,
         help='Ice flexure file to use')
     # tide model to use
-    parser.add_argument('--tide','-T',
+    group.add_argument('--tide','-T',
         metavar='TIDE', type=str,
         choices=get_available_models(),
         help='Tide model to use in correction')
+    # tide model definition file to set an undefined model
+    group.add_argument('--definition-file',
+        type=pathlib.Path,
+        help='Tide model definition file')
+    parser.add_argument('--definition-format',
+        type=str, default='ascii', choices=('ascii', 'json'),
+        help='Format for model definition file')
     # verbosity settings
     # verbose will output information about each output file
     parser.add_argument('--verbose','-V',
@@ -755,6 +770,8 @@ def main():
         adjust_tides_ICESat2_ATL11(args.flexure_file, FILE,
             OUTPUT_DIRECTORY=args.output_directory,
             TIDE_MODEL=args.tide,
+            DEFINITION_FILE=args.definition_file,
+            DEFINITION_FORMAT=args.definition_format,
             VERBOSE=args.verbose,
             MODE=args.mode)
 

tides/compute_tides_ICESat2_ATL03.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_tides_ICESat2_ATL03.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (07/2024)
 Calculates tidal elevations for correcting ICESat-2 photon height data
 Calculated at ATL03 segment level using reference photon geolocation and time
 Segment level corrections can be applied to the individual photon events (PEs)
@@ -66,6 +66,8 @@
     predict.py: predict tidal values using harmonic constants
 
 UPDATE HISTORY:
+    Updated 07/2024: added option to crop to the domain of the input data
+        added option to use JSON format definition files
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 01/2024: made the inferrence of minor constituents an option
@@ -137,6 +139,8 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE,
         ATLAS_FORMAT=None,
         GZIP=True,
         DEFINITION_FILE=None,
+        DEFINITION_FORMAT='ascii',
+        CROP=False,
         METHOD='spline',
         EXTRAPOLATE=False,
         CUTOFF=None,
@@ -152,7 +156,8 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE,
 
     # get parameters for tide model
     if DEFINITION_FILE is not None:
-        model = pyTMD.io.model(tide_dir).from_file(DEFINITION_FILE)
+        model = pyTMD.io.model(tide_dir).from_file(DEFINITION_FILE,
+            format=DEFINITION_FORMAT)
     else:
         model = pyTMD.io.model(tide_dir, format=ATLAS_FORMAT,
             compressed=GZIP).elevation(TIDE_MODEL)
@@ -197,26 +202,47 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE,
     IS2_atl03_mds,IS2_atl03_attrs,IS2_atl03_beams = \
         is2tk.io.ATL03.read_main(INPUT_FILE, ATTRIBUTES=True)
 
+    # read orbit info for bounding polygons
+    bounding_lon = IS2_atl03_mds['orbit_info']['bounding_polygon_lon1']
+    bounding_lat = IS2_atl03_mds['orbit_info']['bounding_polygon_lat1']
+    # convert bounding polygon coordinates to bounding box
+    BOUNDS = [np.inf, -np.inf, np.inf, -np.inf]
+    BOUNDS[0] = np.minimum(BOUNDS[0], np.min(bounding_lon))
+    BOUNDS[1] = np.maximum(BOUNDS[1], np.max(bounding_lon))
+    BOUNDS[2] = np.minimum(BOUNDS[2], np.min(bounding_lat))
+    BOUNDS[3] = np.maximum(BOUNDS[3], np.max(bounding_lat))
+    # check if bounding polygon is in multiple parts
+    if 'bounding_polygon_lon2' in IS2_atl03_mds['orbit_info']:
+        bounding_lon = IS2_atl03_mds['orbit_info']['bounding_polygon_lon2']
+        bounding_lat = IS2_atl03_mds['orbit_info']['bounding_polygon_lat2']
+        BOUNDS[0] = np.minimum(BOUNDS[0], np.min(bounding_lon))
+        BOUNDS[1] = np.maximum(BOUNDS[1], np.max(bounding_lon))
+        BOUNDS[2] = np.minimum(BOUNDS[2], np.min(bounding_lat))
+        BOUNDS[3] = np.maximum(BOUNDS[3], np.max(bounding_lat))
+
     # read tidal constants
     if model.format in ('OTIS','ATLAS','TMD3'):
         constituents = pyTMD.io.OTIS.read_constants(model.grid_file,
             model.model_file, model.projection, type=model.type,
-            grid=model.format, apply_flexure=APPLY_FLEXURE)
+            grid=model.format, crop=CROP, bounds=BOUNDS,
+            apply_flexure=APPLY_FLEXURE)
         # available model constituents
         c = constituents.fields
     elif (model.format == 'netcdf'):
         constituents = pyTMD.io.ATLAS.read_constants(model.grid_file,
-            model.model_file, type=model.type, compressed=model.compressed)
+            model.model_file, type=model.type, compressed=model.compressed,
+            crop=CROP, bounds=BOUNDS)
         # available model constituents
         c = constituents.fields
     elif (model.format == 'GOT'):
         constituents = pyTMD.io.GOT.read_constants(model.model_file,
-            compressed=model.compressed)
+            compressed=model.compressed, crop=CROP, bounds=BOUNDS)
         # available model constituents
         c = constituents.fields
     elif (model.format == 'FES'):
         constituents = pyTMD.io.FES.read_constants(model.model_file,
-            type=model.type, version=model.version, compressed=model.compressed)
+            type=model.type, version=model.version, compressed=model.compressed,
+            crop=CROP, bounds=BOUNDS)
         # available model constituents
         c = model.constituents
 
@@ -621,7 +647,14 @@ def arguments():
     # tide model definition file to set an undefined model
     group.add_argument('--definition-file',
         type=pathlib.Path,
-        help='Tide model definition file for use as correction')
+        help='Tide model definition file')
+    parser.add_argument('--definition-format',
+        type=str, default='ascii', choices=('ascii', 'json'),
+        help='Format for model definition file')
+    # crop tide model to (buffered) bounds of data
+    parser.add_argument('--crop', '-C',
+        default=False, action='store_true',
+        help='Crop tide model to bounds of data')
     # interpolation method
     parser.add_argument('--interpolate','-I',
         metavar='METHOD', type=str, default='spline',
@@ -670,6 +703,8 @@ def main():
             ATLAS_FORMAT=args.atlas_format,
             GZIP=args.gzip,
             DEFINITION_FILE=args.definition_file,
+            DEFINITION_FORMAT=args.definition_format,
+            CROP=args.crop,
             METHOD=args.interpolate,
             EXTRAPOLATE=args.extrapolate,
             CUTOFF=args.cutoff,