From ebc58764a8a2b3fa69c4d6318a49ebf73b157283 Mon Sep 17 00:00:00 2001 From: tsutterley Date: Fri, 19 Jul 2024 15:05:38 -0700 Subject: [PATCH] feat: added option to crop tides to the domain of the input data feat: added option to use JSON format definition files --- DEM/interp_ATL14_DEM_ICESat2_ATL06.py | 9 +++++ DEM/interp_ATL14_DEM_ICESat2_ATL11.py | 9 +++++ GZ/calculate_GZ_ICESat2_ATL06.py | 17 ---------- GZ/calculate_GZ_ICESat2_ATL11.py | 16 --------- GZ/calculate_GZ_ICESat_GLA12.py | 17 ---------- doc/environment.yml | 2 +- tides/adjust_tides_ICESat2_ATL11.py | 23 +++++++++++-- tides/compute_tides_ICESat2_ATL03.py | 49 +++++++++++++++++++++++---- tides/compute_tides_ICESat2_ATL06.py | 41 ++++++++++++++++++---- tides/compute_tides_ICESat2_ATL07.py | 40 ++++++++++++++++++---- tides/compute_tides_ICESat2_ATL10.py | 42 +++++++++++++++++++---- tides/compute_tides_ICESat2_ATL11.py | 49 +++++++++++++++++++++++---- tides/compute_tides_ICESat2_ATL12.py | 40 ++++++++++++++++++---- tides/compute_tides_ICESat_GLA12.py | 28 +++++++++++---- tides/compute_tides_icebridge_data.py | 32 ++++++++++++----- tides/fit_tides_ICESat2_ATL11.py | 45 +++++++++++++----------- 16 files changed, 328 insertions(+), 131 deletions(-) diff --git a/DEM/interp_ATL14_DEM_ICESat2_ATL06.py b/DEM/interp_ATL14_DEM_ICESat2_ATL06.py index 0724884..2876a89 100644 --- a/DEM/interp_ATL14_DEM_ICESat2_ATL06.py +++ b/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) diff --git a/DEM/interp_ATL14_DEM_ICESat2_ATL11.py b/DEM/interp_ATL14_DEM_ICESat2_ATL11.py index 961106b..27ddd7c 100644 --- a/DEM/interp_ATL14_DEM_ICESat2_ATL11.py +++ b/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) diff --git a/GZ/calculate_GZ_ICESat2_ATL06.py b/GZ/calculate_GZ_ICESat2_ATL06.py index 7423e8a..9058d48 100644 --- a/GZ/calculate_GZ_ICESat2_ATL06.py +++ b/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 diff --git a/GZ/calculate_GZ_ICESat2_ATL11.py b/GZ/calculate_GZ_ICESat2_ATL11.py index 1680a62..47a915e 100644 --- a/GZ/calculate_GZ_ICESat2_ATL11.py +++ b/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 diff --git a/GZ/calculate_GZ_ICESat_GLA12.py b/GZ/calculate_GZ_ICESat_GLA12.py index cd3cba0..b28e1d0 100644 --- a/GZ/calculate_GZ_ICESat_GLA12.py +++ b/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 diff --git a/doc/environment.yml b/doc/environment.yml index 5b19ff1..dc91feb 100644 --- a/doc/environment.yml +++ b/doc/environment.yml @@ -2,7 +2,7 @@ name: gz-docs channels: - conda-forge dependencies: - - docutils<0.18 + - docutils - fontconfig - freetype - graphviz diff --git a/tides/adjust_tides_ICESat2_ATL11.py b/tides/adjust_tides_ICESat2_ATL11.py index 8f1d09d..7719284 100644 --- a/tides/adjust_tides_ICESat2_ATL11.py +++ b/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) diff --git a/tides/compute_tides_ICESat2_ATL03.py b/tides/compute_tides_ICESat2_ATL03.py index 5bd1177..4a9736c 100644 --- a/tides/compute_tides_ICESat2_ATL03.py +++ b/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, diff --git a/tides/compute_tides_ICESat2_ATL06.py b/tides/compute_tides_ICESat2_ATL06.py index abcac19..cd31b4c 100644 --- a/tides/compute_tides_ICESat2_ATL06.py +++ b/tides/compute_tides_ICESat2_ATL06.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat2_ATL06.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat-2 land ice elevation data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -61,6 +61,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 @@ -135,6 +137,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, @@ -150,7 +154,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) @@ -195,26 +200,39 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE, IS2_atl06_mds,IS2_atl06_attrs,IS2_atl06_beams = \ is2tk.io.ATL06.read_granule(INPUT_FILE, ATTRIBUTES=True) + # find geospatial ranges for bounding box + BOUNDS = [np.inf, -np.inf, np.inf, -np.inf] + for gtx in IS2_atl06_beams: + lon = IS2_atl06_mds[gtx]['land_ice_segments']['longitude'] + lat = IS2_atl06_mds[gtx]['land_ice_segments']['latitude'] + BOUNDS[0] = np.minimum(BOUNDS[0], np.min(lon)) + BOUNDS[1] = np.maximum(BOUNDS[1], np.max(lon)) + BOUNDS[2] = np.minimum(BOUNDS[2], np.min(lat)) + BOUNDS[3] = np.maximum(BOUNDS[3], np.max(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 @@ -623,7 +641,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', @@ -672,6 +697,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, diff --git a/tides/compute_tides_ICESat2_ATL07.py b/tides/compute_tides_ICESat2_ATL07.py index 0afe2f3..82f5fbd 100644 --- a/tides/compute_tides_ICESat2_ATL07.py +++ b/tides/compute_tides_ICESat2_ATL07.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat2_ATL07.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat-2 sea ice height data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -59,6 +59,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 @@ -131,6 +133,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, @@ -145,7 +149,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) @@ -188,26 +193,38 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE, IS2_atl07_mds,IS2_atl07_attrs,IS2_atl07_beams = \ is2tk.io.ATL07.read_granule(INPUT_FILE, ATTRIBUTES=True) + # find geospatial ranges for bounding box + BOUNDS = [np.inf, -np.inf, np.inf, -np.inf] + for gtx in IS2_atl07_beams: + lon = IS2_atl07_mds[gtx]['sea_ice_segments']['longitude'] + lat = IS2_atl07_mds[gtx]['sea_ice_segments']['latitude'] + BOUNDS[0] = np.minimum(BOUNDS[0], np.min(lon)) + BOUNDS[1] = np.maximum(BOUNDS[1], np.max(lon)) + BOUNDS[2] = np.minimum(BOUNDS[2], np.min(lat)) + BOUNDS[3] = np.maximum(BOUNDS[3], np.max(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) + grid=model.format, crop=CROP, bounds=BOUNDS) # 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 @@ -649,7 +666,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', @@ -694,6 +718,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, diff --git a/tides/compute_tides_ICESat2_ATL10.py b/tides/compute_tides_ICESat2_ATL10.py index 2ecf9df..464b805 100644 --- a/tides/compute_tides_ICESat2_ATL10.py +++ b/tides/compute_tides_ICESat2_ATL10.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat2_ATL10.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat-2 sea ice height data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -59,6 +59,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 @@ -132,6 +134,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, @@ -146,7 +150,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) @@ -189,26 +194,40 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE, IS2_atl10_mds,IS2_atl10_attrs,IS2_atl10_beams = \ is2tk.io.ATL10.read_granule(INPUT_FILE, ATTRIBUTES=True) + # find geospatial ranges for bounding box + BOUNDS = [np.inf, -np.inf, np.inf, -np.inf] + for gtx in IS2_atl10_beams: + # for each ATL10 group + for group in ['freeboard_beam_segment','leads']: + lon = IS2_atl10_mds[gtx][group]['longitude'] + lat = IS2_atl10_mds[gtx][group]['latitude'] + BOUNDS[0] = np.minimum(BOUNDS[0], np.min(lon)) + BOUNDS[1] = np.maximum(BOUNDS[1], np.max(lon)) + BOUNDS[2] = np.minimum(BOUNDS[2], np.min(lat)) + BOUNDS[3] = np.maximum(BOUNDS[3], np.max(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) + grid=model.format, crop=CROP, bounds=BOUNDS) # 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 @@ -612,7 +631,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', @@ -657,6 +683,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, diff --git a/tides/compute_tides_ICESat2_ATL11.py b/tides/compute_tides_ICESat2_ATL11.py index 14e13a0..0066456 100644 --- a/tides/compute_tides_ICESat2_ATL11.py +++ b/tides/compute_tides_ICESat2_ATL11.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat2_ATL11.py -Written by Tyler Sutterley (06/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat-2 annual land ice height data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -62,6 +62,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 06/2024: added option to not run with crossover measurements Updated 05/2024: use wrapper to importlib for optional dependencies Updated 04/2024: use timescale for temporal operations @@ -121,6 +123,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, @@ -137,7 +141,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) @@ -184,26 +189,47 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE, ATTRIBUTES=True, CROSSOVERS=True) + # read orbit info for bounding polygons + bounding_lon = IS2_atl11_mds['orbit_info']['bounding_polygon_lon1'] + bounding_lat = IS2_atl11_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_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'] + 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 @@ -807,7 +833,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', @@ -860,6 +893,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, diff --git a/tides/compute_tides_ICESat2_ATL12.py b/tides/compute_tides_ICESat2_ATL12.py index 122faa4..06c73b7 100644 --- a/tides/compute_tides_ICESat2_ATL12.py +++ b/tides/compute_tides_ICESat2_ATL12.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat2_ATL12.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat-2 ocean surface height data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -59,6 +59,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 @@ -132,6 +134,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, @@ -146,7 +150,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) @@ -189,26 +194,38 @@ def compute_tides_ICESat2(tide_dir, INPUT_FILE, IS2_atl12_mds,IS2_atl12_attrs,IS2_atl12_beams = \ is2tk.io.ATL12.read_granule(INPUT_FILE, ATTRIBUTES=True) + # find geospatial ranges for bounding box + BOUNDS = [np.inf, -np.inf, np.inf, -np.inf] + for gtx in IS2_atl12_beams: + lon = IS2_atl12_mds[gtx]['ssh_segments']['longitude'] + lat = IS2_atl12_mds[gtx]['ssh_segments']['latitude'] + BOUNDS[0] = np.minimum(BOUNDS[0], np.min(lon)) + BOUNDS[1] = np.maximum(BOUNDS[1], np.max(lon)) + BOUNDS[2] = np.minimum(BOUNDS[2], np.min(lat)) + BOUNDS[3] = np.maximum(BOUNDS[3], np.max(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) + grid=model.format, crop=CROP, bounds=BOUNDS) # 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 @@ -612,7 +629,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', @@ -657,6 +681,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, diff --git a/tides/compute_tides_ICESat_GLA12.py b/tides/compute_tides_ICESat_GLA12.py index a023a4b..af80d57 100755 --- a/tides/compute_tides_ICESat_GLA12.py +++ b/tides/compute_tides_ICESat_GLA12.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_ICESat_GLA12.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting ICESat/GLAS L2 GLA12 Antarctic and Greenland Ice Sheet elevation data @@ -65,6 +65,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 @@ -134,6 +136,8 @@ def compute_tides_ICESat(tide_dir, INPUT_FILE, ATLAS_FORMAT=None, GZIP=True, DEFINITION_FILE=None, + DEFINITION_FORMAT='ascii', + CROP=False, METHOD='spline', EXTRAPOLATE=False, CUTOFF=None, @@ -149,7 +153,8 @@ def compute_tides_ICESat(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) @@ -234,25 +239,25 @@ def compute_tides_ICESat(tide_dir, INPUT_FILE, if model.format in ('OTIS','ATLAS','TMD3'): amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon_40HZ, lat_40HZ, model.grid_file, model.model_file, model.projection, - type=model.type, method=METHOD, extrapolate=EXTRAPOLATE, + type=model.type, crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, grid=model.format, apply_flexure=APPLY_FLEXURE) deltat = np.zeros((n_40HZ)) elif (model.format == 'netcdf'): amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon_40HZ, lat_40HZ, model.grid_file, model.model_file, type=model.type, - method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, compressed=model.compressed) deltat = np.zeros((n_40HZ)) elif (model.format == 'GOT'): amp,ph,c = pyTMD.io.GOT.extract_constants(lon_40HZ, lat_40HZ, - model.model_file, method=METHOD, extrapolate=EXTRAPOLATE, + model.model_file, crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, compressed=model.compressed) # delta time (TT - UT1) deltat = ts.tt_ut1 elif (model.format == 'FES'): amp,ph = pyTMD.io.FES.extract_constants(lon_40HZ, lat_40HZ, model.model_file, type=model.type, version=model.version, - method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, compressed=model.compressed) # available model constituents c = model.constituents @@ -519,7 +524,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', @@ -568,6 +580,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, diff --git a/tides/compute_tides_icebridge_data.py b/tides/compute_tides_icebridge_data.py index 6836542..cd3a95e 100644 --- a/tides/compute_tides_icebridge_data.py +++ b/tides/compute_tides_icebridge_data.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_icebridge_data.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Calculates tidal elevations for correcting Operation IceBridge elevation data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -69,6 +69,8 @@ read_ATM1b_QFIT_binary.py: read ATM1b QFIT binary files (NSIDC version 1) 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 @@ -140,6 +142,8 @@ def compute_tides_icebridge_data(tide_dir, arg, TIDE_MODEL, ATLAS_FORMAT=None, GZIP=True, DEFINITION_FILE=None, + DEFINITION_FORMAT='ascii', + CROP=False, METHOD='spline', EXTRAPOLATE=False, CUTOFF=None, @@ -155,7 +159,8 @@ def compute_tides_icebridge_data(tide_dir, arg, TIDE_MODEL, # 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) @@ -256,25 +261,25 @@ def compute_tides_icebridge_data(tide_dir, arg, TIDE_MODEL, if model.format in ('OTIS','ATLAS','TMD3'): amp,ph,D,c = pyTMD.io.OTIS.extract_constants(dinput['lon'], dinput['lat'], model.grid_file, model.model_file, model.projection, - type=model.type, method=METHOD, extrapolate=EXTRAPOLATE, + type=model.type, crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, grid=model.format, apply_flexure=APPLY_FLEXURE) deltat = np.zeros((file_lines)) elif model.format in ('netcdf'): amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(dinput['lon'], dinput['lat'], - model.grid_file, model.model_file, type=model.type, method=METHOD, - extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, - compressed=model.compressed) + model.grid_file, model.model_file, type=model.type, crop=CROP, + method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + scale=model.scale, compressed=model.compressed) deltat = np.zeros((file_lines)) elif (model.format == 'GOT'): amp,ph,c = pyTMD.io.GOT.extract_constants(dinput['lon'], dinput['lat'], - model.model_file, method=METHOD, extrapolate=EXTRAPOLATE, + model.model_file, crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, compressed=model.compressed) # delta time (TT - UT1) deltat = ts.tt_ut1 elif (model.format == 'FES'): amp,ph = pyTMD.io.FES.extract_constants(dinput['lon'], dinput['lat'], model.model_file, type=model.type, version=model.version, - method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + crop=CROP, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, compressed=model.compressed) # available model constituents c = model.constituents @@ -417,7 +422,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', @@ -465,6 +477,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, diff --git a/tides/fit_tides_ICESat2_ATL11.py b/tides/fit_tides_ICESat2_ATL11.py index bf23831..df49faa 100644 --- a/tides/fit_tides_ICESat2_ATL11.py +++ b/tides/fit_tides_ICESat2_ATL11.py @@ -1,30 +1,13 @@ #!/usr/bin/env python u""" fit_tides_ICESat2_ATL11.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (07/2024) Fits tidal amplitudes to ICESat-2 data in ice sheet grounding zones COMMAND LINE OPTIONS: -D X, --directory X: Working data directory -O X, --output-directory X: input/output data directory -T X, --tide X: Tide model to use in correction - CATS0201 - CATS2008 - CATS2022 - 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 @@ -52,6 +35,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 12/2023: don't have a default tide model in arguments @@ -101,6 +85,8 @@ def common_reference_points(XT, AT): def fit_tides_ICESat2(tide_dir, INPUT_FILE, OUTPUT_DIRECTORY=None, TIDE_MODEL=None, + DEFINITION_FILE=None, + DEFINITION_FORMAT='ascii', REANALYSIS=None, VERBOSE=False, MODE=0o775 @@ -111,7 +97,11 @@ def fit_tides_ICESat2(tide_dir, INPUT_FILE, logging.basicConfig(level=loglevel) # get tide model parameters - model = pyTMD.io.model(tide_dir, verify=False).elevation(TIDE_MODEL) + if DEFINITION_FILE is not None: + model = pyTMD.io.model(tide_dir, verify=False).from_file(DEFINITION_FILE, + format=DEFINITION_FORMAT) + else: + model = pyTMD.io.model(tide_dir, verify=False).elevation(TIDE_MODEL) # log input file logging.info(f'{str(INPUT_FILE)} -->') @@ -985,6 +975,8 @@ 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='+', help='ICESat-2 ATL11 file to run') @@ -1001,6 +993,19 @@ def arguments(): metavar='TIDE', type=str, choices=get_available_models(), help='Tide model to use in correction') + # tide model to use + 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') + # inverse barometer response to use parser.add_argument('--reanalysis','-R', metavar='REANALYSIS', type=str, help='Reanalysis model to use in inverse-barometer correction') @@ -1027,6 +1032,8 @@ def main(): fit_tides_ICESat2(args.directory, FILE, OUTPUT_DIRECTORY=args.output_directory, TIDE_MODEL=args.tide, + DEFINITION_FILE=args.definition_file, + DEFINITION_FORMAT=args.definition_format, REANALYSIS=args.reanalysis, VERBOSE=args.verbose, MODE=args.mode)