feat: add tide adjustment mosaic program

doc/source/api_reference/io/raster.rst

@@ -0,0 +1,15 @@
+=========
+io.raster
+=========
+
+Utilities for operating on raster data
+
+`Source code`__
+
+.. __: https://github.com/tsutterley/Grounding-Zones/blob/main/grounding_zones/io/raster.py
+
+General Attributes and Methods
+==============================
+
+.. autoclass:: grounding_zones.io.raster
+   :members:

doc/source/api_reference/mosaic.rst

@@ -9,8 +9,8 @@ Utilities for creating spatial mosaics
 .. __: https://github.com/tsutterley/Grounding-Zones/blob/main/grounding_zones/mosaic.py
 
 
-General Methods
-===============
+General Attributes and Methods
+==============================
 
 .. autoclass:: grounding_zones.mosaic
    :members:

doc/source/api_reference/tides/mosaic_tide_adjustment.rst

@@ -0,0 +1,19 @@
+=========================
+mosaic_tide_adjustment.py
+=========================
+
+- Creates a mosaic of interpolated tidal adjustment scale factors
+
+`Source code`__
+
+.. __: https://github.com/tsutterley/Grounding-Zones/blob/main/tides/mosaic_tide_adjustment.py
+
+Calling Sequence
+################
+
+.. argparse::
+    :filename: mosaic_tide_adjustment.py
+    :func: arguments
+    :prog: mosaic_tide_adjustment.py
+    :nodescription:
+    :nodefault:

doc/source/index.rst

@@ -18,6 +18,7 @@ Python Tools for Estimating Ice Sheet Grounding Zone Locations with data from NA
     :caption: API Reference
 
     api_reference/io/icebridge.rst
+    api_reference/io/raster.rst
     api_reference/crs.rst
     api_reference/mosaic.rst
     api_reference/utilities.rst
@@ -136,3 +137,4 @@ Python Tools for Estimating Ice Sheet Grounding Zone Locations with data from NA
     api_reference/tides/compute_tides_ICESat2_ATL12.rst
     api_reference/tides/fit_tides_ICESat2_ATL11.rst
     api_reference/tides/interpolate_tide_adjustment.rst
+    api_reference/tides/mosaic_tide_adjustment.rst

grounding_zones/io/__init__.py

@@ -1 +1,2 @@
 from .icebridge import *
+from .raster import *

grounding_zones/io/raster.py

@@ -0,0 +1,220 @@
+#!/usr/bin/env python
+u"""
+raster.py
+Written by Tyler Sutterley (10/2023)
+
+Utilities for operating on raster data
+
+PYTHON DEPENDENCIES:
+    numpy: Scientific Computing Tools For Python
+        https://numpy.org
+        https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
+    netCDF4: Python interface to the netCDF C library
+        https://unidata.github.io/netcdf4-python/netCDF4/index.html
+    h5py: Pythonic interface to the HDF5 binary data format
+        https://www.h5py.org/
+    gdal: Pythonic interface to the Geospatial Data Abstraction Library (GDAL)
+        https://pypi.python.org/pypi/GDAL
+    PyYAML: YAML parser and emitter for Python
+        https://github.com/yaml/pyyaml
+
+PROGRAM DEPENDENCIES:
+    spatial.py: utilities for reading and writing spatial data
+
+UPDATE HISTORY:
+    Written 10/2023
+"""
+import warnings
+import numpy as np
+import scipy.interpolate
+# attempt imports
+try:
+    import pyproj
+except (AttributeError, ImportError, ModuleNotFoundError) as exc:
+    warnings.warn("pyproj not available", ImportWarning)
+try:
+    import pyTMD.spatial
+except (AttributeError, ImportError, ModuleNotFoundError) as exc:
+    warnings.warn("pyTMD not available", ImportWarning)
+
+class raster:
+    """Utilities for using raster files
+    """
+    np.seterr(invalid='ignore')
+    def __init__(self, **kwargs):
+        # set default class attributes
+        self.attributes=dict()
+
+    # PURPOSE: read a raster file
+    def from_file(self, input_file, format=None, **kwargs):
+        """
+        Read a raster file from an input format
+
+        Parameters
+        ----------
+        input_file: str
+            path or memory map for raster file
+        format: str
+            format of input file
+        **kwargs: dict
+            Keyword arguments for file reader
+        """
+        dinput = pyTMD.spatial.from_file(
+            filename=input_file,
+            format=format,
+            **kwargs
+        )
+        # separate dimensions from fields
+        self.dims, self.fields = ([], [])
+        # convert from dictionary to class attributes
+        for key,val in dinput.items():
+            if key in ('x','y','time'):
+                self.dims.append(key)
+            elif key in ('attributes','filename'):
+                pass
+            else:
+                self.fields.append(key)
+            # set attribute
+            setattr(self, key, val)
+        # return the raster object
+        return self
+
+    def warp(self, xout, yout, order=0, reducer=np.ceil):
+        """
+        Interpolate raster data to a new grid
+
+        Parameters
+        ----------
+        datain: np.ndarray
+            input data grid to be interpolated
+        xin: np.ndarray
+            input x-coordinate array (monotonically increasing)
+        yin: np.ndarray
+            input y-coordinate array (monotonically increasing)
+        xout: np.ndarray
+            output x-coordinate array
+        yout: np.ndarray
+            output y-coordinate array
+        order: int, default 0
+            interpolation order
+
+                - ``0``: nearest-neighbor interpolation
+                - ``k``: bivariate spline interpolation of degree k
+        reducer: obj, default np.ceil
+            operation for converting mask to boolean
+        """
+        temp = raster()
+        if (order == 0):
+            # interpolate with nearest-neighbors
+            xcoords = (len(self.x)-1)*(xout-self.x[0])/(self.x[-1]-self.x[0])
+            ycoords = (len(self.y)-1)*(yout-self.y[0])/(self.y[-1]-self.y[0])
+            xcoords = np.clip(xcoords,0,len(self.x)-1)
+            ycoords = np.clip(ycoords,0,len(self.y)-1)
+            XI = np.around(xcoords).astype(np.int32)
+            YI = np.around(ycoords).astype(np.int32)
+            temp.data = self.data[YI, XI]
+            if hasattr(self.data, 'mask'):
+                temp.mask = reducer(self.data.mask[YI, XI]).astype(bool)
+            elif hasattr(self, 'mask'):
+                temp.mask = reducer(self.mask[YI, XI]).astype(bool)
+        else:
+            # interpolate with bivariate spline approximations
+            s1 = scipy.interpolate.RectBivariateSpline(self.x, self.y,
+                self.data.T, kx=order, ky=order)
+            temp.data = s1.ev(xout, yout)
+            # try to interpolate the mask
+            if hasattr(self.data, 'mask'):
+                s2 = scipy.interpolate.RectBivariateSpline(self.x, self.y,
+                    self.data.mask.T, kx=order, ky=order)
+                temp.mask = reducer(s2.ev(xout, yout)).astype(bool)
+            elif hasattr(self, 'mask'):
+                s2 = scipy.interpolate.RectBivariateSpline(self.x, self.y,
+                    self.mask.T, kx=order, ky=order)
+                temp.mask = reducer(s2.ev(xout, yout)).astype(bool)
+        # return the interpolated data on the output grid
+        return temp
+
+    def get_latlon(self, srs_proj4=None, srs_wkt=None, srs_epsg=None):
+        """
+        Get the latitude and longitude of grid cells
+
+        Parameters
+        ----------
+        srs_proj4: str or NoneType, default None
+            PROJ4 projection string
+        srs_wkt: str or NoneType, default None
+            Well-Known Text (WKT) projection string
+        srs_epsg: int or NoneType, default None
+            EPSG projection code
+
+        Returns
+        -------
+        longitude: np.ndarray
+            longitude coordinates of grid cells
+        latitude: np.ndarray
+            latitude coordinates of grid cells
+        """
+        # set the spatial projection reference information
+        if srs_proj4 is not None:
+            source = pyproj.CRS.from_proj4(srs_proj4)
+        elif srs_wkt is not None:
+            source = pyproj.CRS.from_wkt(srs_wkt)
+        elif srs_epsg is not None:
+            source = pyproj.CRS.from_epsg(srs_epsg)
+        else:
+            source = pyproj.CRS.from_string(self.projection)
+        # target spatial reference (WGS84 latitude and longitude)
+        target = pyproj.CRS.from_epsg(4326)
+        # create transformation
+        transformer = pyproj.Transformer.from_crs(source, target,
+            always_xy=True)
+        # create meshgrid of points in original projection
+        x, y = np.meshgrid(self.x, self.y)
+        # convert coordinates to latitude and longitude
+        self.lon, self.lat = transformer.transform(x, y)
+        return self
+
+    def copy(self):
+        """
+        Copy a ``raster`` object to a new ``raster`` object
+        """
+        temp = raster()
+        # copy attributes or update attributes dictionary
+        if isinstance(self.attributes, list):
+            setattr(temp,'attributes', self.attributes)
+        elif isinstance(self.attributes, dict):
+            temp.attributes.update(self.attributes)
+        # get dimensions and field names
+        temp.dims = self.dims
+        temp.fields = self.fields
+        # assign variables to self
+        for key in [*self.dims, *self.fields]:
+            try:
+                setattr(temp, key, getattr(self, key))
+            except AttributeError:
+                pass
+        return temp
+
+    def flip(self, axis=0):
+        """
+        Reverse the order of data and dimensions along an axis
+
+        Parameters
+        ----------
+        axis: int, default 0
+            axis to reorder
+        """
+        # output spatial object
+        temp = self.copy()
+        # copy dimensions and reverse order
+        if (axis == 0):
+            temp.y = temp.y[::-1].copy()
+        elif (axis == 1):
+            temp.x = temp.x[::-1].copy()
+        # attempt to reverse possible data variables
+        for key in self.fields:
+            try:
+                setattr(temp, key, np.flip(getattr(self, key), axis=axis))
+            except Exception as exc:
+                pass
+        return temp

tides/interpolate_tide_adjustment.py

@@ -7,9 +7,13 @@
 COMMAND LINE OPTIONS:
     --help: list the command line options
     -O X, --output-directory X: input/output data directory
+    -H X, --hemisphere X: Region of interest to run
     -W X, --width: Width of tile grid
     -s X, --subset: Width of interpolation subset
     -S X, --spacing X: Output grid spacing
+    -P X, --pad X: Tile pad for creating mosaics
+    -T X, --tide X: Tide model used in correction
+    -I X, --interpolate X: Interpolation method
     -t X, --tension X: Biharmonic spline tension
     -w X, --smooth X: Radial basis function smoothing weight
     -e X, --epsilon X: Radial basis function adjustable constant
@@ -32,7 +36,6 @@
         single implicit import of grounding zone tools
     Updated 07/2022: place some imports within try/except statements
     Updated 06/2022: use argparse descriptions within documentation
-        read mask files to not interpolate over grounded ice
     Updated 01/2022: added options for using radial basis functions
         wait if HDF5 tile file is unavailable for read or write
     Written 12/2021
@@ -96,13 +99,8 @@ def interpolate_tide_adjustment(tile_file,
     SMOOTH=0,
     EPSILON=0,
     POLYNOMIAL=0,
-    VERBOSE=False,
     MODE=0o775):
 
-    # create logger
-    loglevel = logging.INFO if VERBOSE else logging.CRITICAL
-    logging.basicConfig(level=loglevel)
-
     # input tile data file
     tile_file = pathlib.Path(tile_file).expanduser().absolute()
     tile_file_format = 'E{0:0.0f}_N{1:0.0f}.h5'
@@ -495,8 +493,8 @@ def arguments():
         help='Tile pad for creating mosaics')
     # tide model to use
     parser.add_argument('--tide','-T',
-        metavar='TIDE', type=str, default='CATS2022',
-        help='Tide model to use in correction')
+        metavar='TIDE', type=str, default='CATS2008-v2023',
+        help='Tide model used in correction')
     # interpolation method
     parser.add_argument('--interpolate','-I',
         type=str, default='radial', choices=('spline','radial'),
@@ -524,7 +522,7 @@ def arguments():
     # permissions mode of the directories and files (number in octal)
     parser.add_argument('--mode','-M',
         type=lambda x: int(x,base=8), default=0o775,
-        help='Local permissions mode of the output mosaic')
+        help='Local permissions mode of the output file')
     # return the parser
     return parser
 
@@ -534,6 +532,10 @@ def main():
     parser = arguments()
     args,_ = parser.parse_known_args()
 
+    # create logger
+    loglevel = logging.INFO if args.verbose else logging.CRITICAL
+    logging.basicConfig(level=loglevel)
+
     # run program for each file
     for FILE in args.infile:
         interpolate_tide_adjustment(FILE,
@@ -549,7 +551,6 @@ def main():
             SMOOTH=args.smooth,
             EPSILON=args.epsilon,
             POLYNOMIAL=args.polynomial,
-            VERBOSE=args.verbose,
             MODE=args.mode)
 
 # run main program