diff --git a/doc/source/api_reference/arguments.rst b/doc/source/api_reference/arguments.rst index 111b3c18..6e6de179 100644 --- a/doc/source/api_reference/arguments.rst +++ b/doc/source/api_reference/arguments.rst @@ -3,7 +3,7 @@ arguments ========= - Calculates the nodal corrections for tidal constituents -- Based on Richard Ray's ``ARGUMENTS`` fortran subroutine +- Originally based on Richard Ray's ``ARGUMENTS`` fortran subroutine Calling Sequence ---------------- @@ -11,11 +11,23 @@ Calling Sequence .. code-block:: python import pyTMD.arguments - pu,pf,G = pyTMD.arguments(MJD, constituents, + pu,pf,G = pyTMD.arguments.arguments(MJD, constituents, deltat=DELTAT, corrections=CORRECTIONS) `Source code`__ .. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/arguments.py -.. autofunction:: pyTMD.arguments +.. autofunction:: pyTMD.arguments.arguments + +.. autofunction:: pyTMD.arguments.minor_arguments + +.. autofunction:: pyTMD.arguments.doodson_number + +.. autofunction:: pyTMD.arguments._arguments_table + +.. autofunction:: pyTMD.arguments._minor_table + +.. autofunction:: pyTMD.arguments._to_doodson_number + +.. autofunction:: pyTMD.arguments._from_doodson_number diff --git a/doc/source/api_reference/astro.rst b/doc/source/api_reference/astro.rst index 230e7f24..ad9001bc 100644 --- a/doc/source/api_reference/astro.rst +++ b/doc/source/api_reference/astro.rst @@ -2,8 +2,7 @@ astro ===== - Astronomical and nutation routines -- Computes the basic astronomical mean longitudes: `S`, `H`, `P`, `N` and `PP` -- Computes astronomical phase angles for the sun and moon: `S`, `H`, `P`, `TAU`, `ZNS` and `PS` +- Computes the basic astronomical mean longitudes and other fundamental orbital parameters - Computes the solar and lunar positions in Earth-Centered Earth-Fixed (ECEF) coordinates Calling Sequence @@ -24,7 +23,7 @@ Calling Sequence .. autofunction:: pyTMD.astro.mean_longitudes -.. autofunction:: pyTMD.astro.phase_angles +.. autofunction:: pyTMD.astro.doodson_arguments .. autofunction:: pyTMD.astro.delaunay_arguments diff --git a/doc/source/api_reference/convert_crs.rst b/doc/source/api_reference/convert_crs.rst deleted file mode 100644 index d2110a60..00000000 --- a/doc/source/api_reference/convert_crs.rst +++ /dev/null @@ -1,36 +0,0 @@ -=========== -convert_crs -=========== - -- Uses `pyproj `_ to convert points to and from the Coordinates Reference Systems (CRS) of tide models - -Calling Sequence ----------------- - -.. code-block:: python - - import pyTMD.convert_crs - x, y = pyTMD.convert_crs(lon, lat, PROJ, 'F') - lon, lat = pyTMD.convert_crs(x, y, PROJ, 'B') - -`Source code`__ - -.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/convert_crs.py - -.. autofunction:: pyTMD.convert_crs - -.. autofunction:: pyTMD.convert_crs.crs_from_input - -.. autofunction:: pyTMD.convert_crs._EPSG3031 - -.. autofunction:: pyTMD.convert_crs._EPSG3413 - -.. autofunction:: pyTMD.convert_crs._CATS2008 - -.. autofunction:: pyTMD.convert_crs._EPSG3976 - -.. autofunction:: pyTMD.convert_crs._PSNorth - -.. autofunction:: pyTMD.convert_crs._EPSG4326 - -.. autofunction:: pyTMD.convert_crs._custom diff --git a/doc/source/api_reference/crs.rst b/doc/source/api_reference/crs.rst new file mode 100644 index 00000000..2b692d01 --- /dev/null +++ b/doc/source/api_reference/crs.rst @@ -0,0 +1,16 @@ +=== +crs +=== + +Coordinates Reference System (CRS) routines + +`Source code`__ + +.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/crs.py + +General Attributes and Methods +============================== + +.. autoclass:: pyTMD.crs + :members: + :private-members: diff --git a/doc/source/api_reference/solve.rst b/doc/source/api_reference/solve.rst new file mode 100644 index 00000000..3d8c5cd1 --- /dev/null +++ b/doc/source/api_reference/solve.rst @@ -0,0 +1,19 @@ +===== +solve +===== + +- Routines for estimating the harmonic constants for ocean tides + +Calling Sequence +---------------- + +.. code-block:: python + + import pyTMD.solve + amp, phase = pyTMD.solve.constants(time, h, con) + +`Source code`__ + +.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/solve.py + +.. autofunction:: pyTMD.solve.constants diff --git a/doc/source/getting_started/Background.rst b/doc/source/getting_started/Background.rst index 7d7ee3a6..0434ab29 100644 --- a/doc/source/getting_started/Background.rst +++ b/doc/source/getting_started/Background.rst @@ -8,7 +8,7 @@ Ocean and Load Tides The rise and fall of the oceanic tides are a major source of the vertical variability of the ocean surface. Ocean tides are typically observed using float gauges, GPS stations, gravimeters, tiltmeters, pressure recorders, and satellite altimeters. For each of these measurements, it is important to note the `vertical datum of the measurement technique `_. -Ocean tides are driven by gravitational undulations due to the relative positions of the Earth, moon and sun, and the centripetal acceleration due to the Earth's rotation [Meeus1998]_. +Ocean tides are driven by gravitational undulations due to the relative positions of the Earth, moon and sun, and the centripetal acceleration due to the Earth's rotation [Doodson1921]_ [Meeus1998]_. A secondary tidal effect, known as load tides, is due to the elastic response of the Earth's crust to ocean tidal loading, which produces deformation of both the sea floor and adjacent land areas. Tidal oscillations for both ocean and load tides can be decomposed into a series of tidal constituents (or partial tides) of particular frequencies. @@ -21,7 +21,7 @@ or 3) an unconstrained hydrodynamic model [Stammer2014]_. Solid Earth Tides ################# -Similar to ocean tides, solid Earth tides are tidal deformations due to gravitational undulations based on the relative positions of the Earth, moon and sun [Agnew2015]_ [Meeus1998]_ [Montenbruck1989]_. +Similar to ocean tides, solid Earth tides are tidal deformations due to gravitational undulations based on the relative positions of the Earth, moon and sun [Agnew2015]_ [Doodson1921]_ [Meeus1998]_ [Montenbruck1989]_. However, while ocean tides are apparent to observers on the coast, solid Earth tides are typically more difficult to observe due to the reference frame of the observer moving. The total gravitational potential at a position on the Earth's surface due to a celestial object is directly related to the distance between the Earth and the object, and the mass of that object [Agnew2015]_ [Wahr1981]_. Analytical approximate positions for the sun and moon can be calculated within ``pyTMD``, and high-resolution numerical ephemerides for the sun and moon can be downloaded from the `Jet Propulsion Laboratory `_. @@ -81,6 +81,30 @@ As opposed to simple vertical offsets, changing the terrestial reference system This involves converting from a geographic coordinate system into a Cartesian coordinate system. Within ``pyTMD``, solid Earth tides are calculated using ECEF coordinates, and pole tides are calculated using geocentric coordinates. +Nutation is the periodic oscillation of the Earth's rotation axis around its mean position. +Nutation is often split into two components, the nutation in longitude and the nutation in obliquity. +The angle between the equator and the orbital plane of Earth around the Sun (the ecliptic) defines the inclination of the Earth's rotation axis (obliquity of the ecliptic). + +Time +#### + +The Julian Day (JD) is the continuous count of days starting at noon on January 1, 4713 B.C (-4712-01-01T12:00:00). +The Modified Julian Day (MJD) differs from the Julian Day by reducing the number of digits for modern periods, and by beginning at midnight. +The MJD is calculated from the Julian Day by + +.. math:: + :label: 4 + + MJD = JD - 2400000.5 + +The start of the Modified Julian Day calendar is 1858-11-17T00:00:00. +Time in Julian centuries (36525 days) are calculated relative to noon on January 1, 2000 (2000-01-01T12:00:00). + +.. math:: + :label: 5 + + T = \frac{JD - 2451545.0}{36525} + References ########## @@ -90,6 +114,8 @@ References .. [Desai2015] S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from satellite altimetry", *Journal of Geodesy*, 89(12), p1233-1243, (2015). `doi: 10.1007/s00190-015-0848-7 `_ +.. [Doodson1921] A. T. Doodson and H. Lamb, "The harmonic development of the tide-generating potential", *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 100(704), 305--329, (1921). `doi: 10.1098/rspa.1921.0088 `_ + .. [Egbert2002] G. D. Egbert and S. Y. Erofeeva, "Efficient Inverse Modeling of Barotropic Ocean Tides", *Journal of Atmospheric and Oceanic Technology*, 19(2), 183--204, (2002). `doi: 10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2`__ .. [Mathews1997] P. M. Mathews, V. Dehant and J. M. Gipson, "Tidal station displacements", *Journal of Geophysical Research: Solid Earth*, 102(B9), 20469--20477, (1997). `doi: 10.1029/97JB01515 `_ diff --git a/doc/source/getting_started/Getting-Started.rst b/doc/source/getting_started/Getting-Started.rst index 12f1d68f..3324ae31 100644 --- a/doc/source/getting_started/Getting-Started.rst +++ b/doc/source/getting_started/Getting-Started.rst @@ -3,7 +3,7 @@ Getting Started =============== This documentation is intended to explain how to compute ocean, solid Earth, load and pole tide variations using the set of ``pyTMD`` programs. -See the `background material <./Background.html>`_ for more information on the theory and methods used in ``pyTMD``. +See the `background material <./Background.html>`_ and `glossary <./Glossary.html>`_ for more information on the theory and methods used in ``pyTMD``. Tide Model Formats ################## @@ -43,15 +43,15 @@ Presently, the following models and their directories parameterized within ``pyT - TOPEX/POSEIDON global tide models [Egbert2002]_ + * TPXO7.2: ``/TPXO7.2_tmd/`` + * TPXO7.2_load: ``/TPXO7.2_load/`` + * `TPXO8-atlas `_: ``/tpxo8_atlas/`` + * `TPXO9.1 `_: ``/TPXO9.1/DATA/`` * `TPXO9-atlas `_: ``/TPXO9_atlas/`` * `TPXO9-atlas-v2 `_: ``/TPXO9_atlas_v2/`` * `TPXO9-atlas-v3 `_: ``/TPXO9_atlas_v3/`` * `TPXO9-atlas-v4 `_: ``/TPXO9_atlas_v4/`` * `TPXO9-atlas-v5 `_: ``/TPXO9_atlas_v5/`` - * `TPXO9.1 `_: ``/TPXO9.1/DATA/`` - * `TPXO8-atlas `_: ``/tpxo8_atlas/`` - * TPXO7.2: ``/TPXO7.2_tmd/`` - * TPXO7.2_load: ``/TPXO7.2_load/`` - Global Ocean Tide models [Ray1999]_ diff --git a/doc/source/getting_started/Glossary.rst b/doc/source/getting_started/Glossary.rst new file mode 100644 index 00000000..adfa4bc1 --- /dev/null +++ b/doc/source/getting_started/Glossary.rst @@ -0,0 +1,108 @@ +======== +Glossary +======== + +.. glossary:: + + Anomaly + angular distance between the perihelion and the position of a celestial body + + Aphelion + point of an orbit where a celestial body is furthest from the sun + + Ascending Node + point of an orbit where a celestial body intersects the ecliptic, and the latitude coordinate is increasing + + Barycenter + center of mass of a system of bodies + + Body Tide + see :term:`Solid Earth Tide` + + Chandler Wobble + small, semi-periodic deviations in the motion of the pole of rotation + + Constituent + tidal oscillation corresponding with a specific periodic motion of the Earth, Sun and Moon + + Descending Node + point of an orbit where a celestial body intersects the ecliptic, and the latitude coordinate is decreasing + + Diurnal Tide + tidal oscillations with a near-daily period + + Ecliptic + plane of the orbit of the Earth around the sun + + Ephemeris + table of positions and velocities of a celestial body at given instances in time + + Ephemerides + plural of :term:`Ephemeris` + + Epoch + fixed point in time used as a reference value + + Load Tide + elastic deformation of the solid Earth due to ocean and atmospheric tides + + Long Period Tide + tidal oscillations with periods much greater than one day + + Love and Shida Numbers + dimensionless parameters relating the vertical (`h`), horizontal (`l`) and gravitational (`k`) elastic responses to tidal loading + + Mean Tide + model with both direct and indirect permanent tidal effects retained + + Nutation + short-period oscillations in the motion of the pole of rotation + + Obliquity + angle between the equatorial and orbital planes + + Ocean Tide + periodic movement in the level of sea surface due to gravitational and rotational forces + + Perigee + point of an orbit where a celestial body is closest to Earth + + Perihelion + point of an orbit where a celestial body is closest to the sun + + Period + time it takes to make one complete revolution + + Permanent Tide + permanent deformation of the Earth caused by the presence of the Sun and the Moon + + see :term:`Mean Tide`, :term:`Tide-Free`, and :term:`Zero Tide` + + Pole Tide + apparent tide due to variations in the Earth's figure axis about its mean + + Semi-diurnal Tide + tidal oscillations with an approximate half-day period + + Solid Earth Tide + deformation of the solid Earth due to gravitational forces + + Tidal Current + horizontal movement of water due to periodic forces + + Tidal Species + classification of tidal constituents based on period + + see :term:`Semi-diurnal Tide`, :term:`Diurnal Tide`, and :term:`Long Period Tide` + + Tidal Stream + see :term:`Tidal Current` + + Tide-Free + model with direct and indirect permanent tidal effects removed + + Vertical Datum + reference coordinate surface used for vertical positions + + Zero Tide + model with permanent direct tidal effects removed, but indirect loading effects retained diff --git a/doc/source/getting_started/Install.rst b/doc/source/getting_started/Install.rst index d3598e5c..65dd19b0 100644 --- a/doc/source/getting_started/Install.rst +++ b/doc/source/getting_started/Install.rst @@ -34,7 +34,7 @@ the `Python Package Index (pypi) `_, and from `conda-forge `_. -The simplest installation for most users will likely be using ``conda``: +The simplest installation for most users will likely be using ``conda`` or ``mamba``: .. code-block:: bash diff --git a/doc/source/getting_started/Resources.rst b/doc/source/getting_started/Resources.rst index 9f00a67f..0662159f 100644 --- a/doc/source/getting_started/Resources.rst +++ b/doc/source/getting_started/Resources.rst @@ -29,6 +29,7 @@ Software - `ESR Tide Model Driver (TMD) Matlab Toolbox `_ - `TMD Matlab Toolbox Repository `_ +- `TMD3 Matlab Toolbox Repository `_ - `Antarctic Tide Gauge (AntTG) Database Tools `_ - `OSU Tidal Prediction Software (OTPS) `_ - `FES (Finite Element Solution) tide model software `_ diff --git a/doc/source/index.rst b/doc/source/index.rst index 48a1f9ac..d39ae9f1 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -12,6 +12,7 @@ ocean, load, solid Earth and pole tides getting_started/Install.rst getting_started/Getting-Started.rst getting_started/Background.rst + getting_started/Glossary.rst getting_started/Contributing.rst getting_started/Code-of-Conduct.rst getting_started/Resources.rst @@ -34,7 +35,7 @@ ocean, load, solid Earth and pole tides api_reference/check_points.rst api_reference/compute_tide_corrections.rst api_reference/constants.rst - api_reference/convert_crs.rst + api_reference/crs.rst api_reference/ellipse.rst api_reference/eop.rst api_reference/interpolate.rst @@ -46,6 +47,7 @@ ocean, load, solid Earth and pole tides api_reference/io/model.rst api_reference/io/ocean_pole_tide.rst api_reference/predict.rst + api_reference/solve.rst api_reference/spatial.rst api_reference/time.rst api_reference/utilities.rst diff --git a/notebooks/Check Tide Map.ipynb b/notebooks/Check Tide Map.ipynb index 1dd2f7fd..58484fec 100644 --- a/notebooks/Check Tide Map.ipynb +++ b/notebooks/Check Tide Map.ipynb @@ -30,7 +30,7 @@ "\n", "#### Program Dependencies\n", "\n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines\n", "- `io.model.py`: retrieves tide model parameters for named tide models\n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from netcdf models \n", @@ -55,17 +55,15 @@ "source": [ "from __future__ import print_function\n", "\n", - "import os\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "import ipywidgets as widgets\n", "import IPython.display\n", "\n", "# import tide programs\n", + "import pyTMD.crs\n", "import pyTMD.io\n", "import pyTMD.time\n", "import pyTMD.tools\n", - "import pyTMD.convert_crs\n", "# autoreload\n", "%load_ext autoreload\n", "%autoreload 2" @@ -83,7 +81,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'GOT4.10'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.atlas,\n", @@ -160,7 +158,7 @@ " # if reading a single OTIS solution\n", " xi,yi,hz,mz,iob,dt = pyTMD.io.OTIS.read_otis_grid(model.grid_file)\n", " # convert coordinate systems of input latitude and longitude\n", - " x,y = pyTMD.convert_crs(np.atleast_1d(LON), np.atleast_1d(LAT),\n", + " x,y = pyTMD.crs().convert(np.atleast_1d(LON), np.atleast_1d(LAT),\n", " model.projection, 'F')\n", " # adjust longitudinal convention of input latitude and longitude\n", " # to fit tide model convention (if global)\n", diff --git a/notebooks/Plot ATLAS Compact.ipynb b/notebooks/Plot ATLAS Compact.ipynb index a1162d66..c202786b 100644 --- a/notebooks/Plot ATLAS Compact.ipynb +++ b/notebooks/Plot ATLAS Compact.ipynb @@ -22,7 +22,7 @@ "\n", "#### Program Dependencies\n", "\n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models\n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `predict.py`: predict tidal values using harmonic constants \n", @@ -49,7 +49,6 @@ "import os\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "import ipywidgets as widgets\n", "\n", "import pyTMD.time\n", "import pyTMD.tools\n", @@ -73,7 +72,7 @@ "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'TPXO8-atlas'\n", "TMDwidgets.compress.value = False\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.compress,\n", diff --git a/notebooks/Plot Antarctic Tidal Currents.ipynb b/notebooks/Plot Antarctic Tidal Currents.ipynb index e47c67a8..bf63f7d9 100644 --- a/notebooks/Plot Antarctic Tidal Currents.ipynb +++ b/notebooks/Plot Antarctic Tidal Currents.ipynb @@ -30,7 +30,7 @@ "\n", "- `arguments.py`: load the nodal corrections for tidal constituents \n", "- `astro.py`: computes the basic astronomical mean longitudes \n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models\n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netcdf models \n", @@ -65,7 +65,6 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "import cartopy.crs as ccrs\n", - "import ipywidgets as widgets\n", "from IPython.display import HTML\n", "\n", "# import tide programs\n", @@ -103,7 +102,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'CATS2008'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.atlas,\n", diff --git a/notebooks/Plot Antarctic Tide Range.ipynb b/notebooks/Plot Antarctic Tide Range.ipynb index ee991a58..7ddd3713 100644 --- a/notebooks/Plot Antarctic Tide Range.ipynb +++ b/notebooks/Plot Antarctic Tide Range.ipynb @@ -30,7 +30,7 @@ "\n", "#### Program Dependencies\n", "\n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models \n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netcdf models \n", @@ -59,7 +59,6 @@ "matplotlib.rcParams['axes.linewidth'] = 2.0\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", - "import ipywidgets as widgets\n", "\n", "# import tide programs\n", "import pyTMD.io\n", @@ -92,7 +91,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'GOT4.10'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.atlas,\n", diff --git a/notebooks/Plot Arctic Ocean Map.ipynb b/notebooks/Plot Arctic Ocean Map.ipynb index b14d2aa3..a709e5b1 100644 --- a/notebooks/Plot Arctic Ocean Map.ipynb +++ b/notebooks/Plot Arctic Ocean Map.ipynb @@ -32,7 +32,7 @@ "\n", "- `arguments.py`: load the nodal corrections for tidal constituents \n", "- `astro.py`: computes the basic astronomical mean longitudes \n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models \n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netcdf models \n", @@ -67,7 +67,6 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "import cartopy.crs as ccrs\n", - "import ipywidgets as widgets\n", "from IPython.display import HTML\n", "\n", "# import tide programs\n", @@ -105,7 +104,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'TPXO8-atlas'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.compress,\n", diff --git a/notebooks/Plot Ross Ice Shelf Map.ipynb b/notebooks/Plot Ross Ice Shelf Map.ipynb index 914b9a82..58d8a7a9 100644 --- a/notebooks/Plot Ross Ice Shelf Map.ipynb +++ b/notebooks/Plot Ross Ice Shelf Map.ipynb @@ -32,7 +32,7 @@ "\n", "- `arguments.py`: load the nodal corrections for tidal constituents \n", "- `astro.py`: computes the basic astronomical mean longitudes \n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models\n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netcdf models \n", @@ -67,7 +67,6 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "import cartopy.crs as ccrs\n", - "import ipywidgets as widgets\n", "from IPython.display import HTML\n", "\n", "# import tide programs\n", @@ -105,7 +104,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'GOT4.10'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.atlas,\n", diff --git a/notebooks/Plot Tide Forecasts.ipynb b/notebooks/Plot Tide Forecasts.ipynb index e2f927ad..4b0e1c74 100644 --- a/notebooks/Plot Tide Forecasts.ipynb +++ b/notebooks/Plot Tide Forecasts.ipynb @@ -32,7 +32,7 @@ "\n", "- `arguments.py`: load the nodal corrections for tidal constituents \n", "- `astro.py`: computes the basic astronomical mean longitudes \n", - "- `convert_crs.py`: convert points to and from Coordinates Reference Systems \n", + "- `crs.py`: Coordinate Reference System (CRS) routines \n", "- `io.model.py`: retrieves tide model parameters for named tide models\n", "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models \n", "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netCDF4 tide models \n", @@ -65,7 +65,6 @@ "import datetime\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "import ipywidgets as widgets\n", "import IPython.display\n", "\n", "# import tide programs\n", @@ -92,7 +91,7 @@ "TMDwidgets = pyTMD.tools.widgets()\n", "TMDwidgets.model.options = model_list\n", "TMDwidgets.model.value = 'GOT4.10'\n", - "widgets.VBox([\n", + "TMDwidgets.VBox([\n", " TMDwidgets.directory,\n", " TMDwidgets.model,\n", " TMDwidgets.atlas,\n", diff --git a/pyTMD/__init__.py b/pyTMD/__init__.py index c6185e02..55141db0 100644 --- a/pyTMD/__init__.py +++ b/pyTMD/__init__.py @@ -11,17 +11,18 @@ Documentation is available at https://pytmd.readthedocs.io """ +import pyTMD.arguments import pyTMD.astro import pyTMD.eop import pyTMD.interpolate import pyTMD.predict +import pyTMD.solve import pyTMD.spatial import pyTMD.time import pyTMD.tools import pyTMD.utilities import pyTMD.version from pyTMD import io -from pyTMD.arguments import arguments from pyTMD.check_points import check_points from pyTMD.compute_tide_corrections import ( compute_corrections, @@ -36,6 +37,7 @@ _ellipsoids ) from pyTMD.convert_crs import convert_crs +from pyTMD.crs import crs from pyTMD.ellipse import ellipse # Deprecated functions diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index d1313bf7..e4b7f2d2 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -1,19 +1,19 @@ #!/usr/bin/env python u""" arguments.py -Written by Tyler Sutterley (08/2023) +Written by Tyler Sutterley (01/2024) Calculates the nodal corrections for tidal constituents Modification of ARGUMENTS fortran subroutine by Richard Ray 03/1999 CALLING SEQUENCE: - pu,pf,G = arguments(MJD, constituents) + pu, pf, G = arguments(MJD, constituents) INPUTS: MJD: Modified Julian Day of input date constituents: tidal constituent IDs OUTPUTS: - pu,pf: nodal corrections for the constituents + pu, pf: nodal corrections for the constituents G: phase correction in degrees OPTIONS: @@ -38,6 +38,11 @@ Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002). UPDATE HISTORY: + Updated 01/2024: add function to create arguments coefficients table + add function to calculate the arguments for minor constituents + include multiples of 90 degrees as variable following Ray 2017 + add function to calculate the Doodson numbers for constituents + Updated 12/2023: made keyword argument for selecting M1 coefficients Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 04/2023: using renamed astro mean_longitudes function function renamed from original load_nodal_corrections @@ -68,18 +73,24 @@ def arguments( **kwargs ): """ - Calculates the nodal corrections for tidal constituents [1]_ [2]_ [3]_ + Calculates the nodal corrections for tidal constituents + [1]_ [2]_ [3]_ [4]_ Parameters ---------- MJD: np.ndarray - modified julian day of input date + modified Julian day of input date constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT models + M1: str, default 'Ray' + coefficients to use for M1 tides + + - ``'Doodson'`` + - ``'Ray'`` Returns ------- @@ -110,110 +121,51 @@ def arguments( # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') + kwargs.setdefault('M1', 'Ray') # constituents array (not all are included in tidal program) - cindex = ['sa','ssa','mm','msf','mf','mt','alpha1','2q1','sigma1','q1', - 'rho1','o1','tau1','m1','chi1','pi1','p1','s1','k1','psi1','phi1', - 'theta1','j1','oo1','2n2','mu2','n2','nu2','m2a','m2','m2b','lambda2', - 'l2','t2','s2','r2','k2','eta2','mns2','2sm2','m3','mk3','s3','mn4', - 'm4','ms4','mk4','s4','s5','m6','s6','s7','s8','m8','mks2','msqm','mtm', - 'n4','eps2','z0'] + cindex = ['sa', 'ssa', 'mm', 'msf', 'mf', 'mt', 'alpha1', '2q1', 'sigma1', + 'q1', 'rho1', 'o1', 'tau1', 'm1', 'chi1', 'pi1', 'p1', 's1', 'k1', + 'psi1', 'phi1', 'theta1', 'j1', 'oo1', '2n2', 'mu2', 'n2', 'nu2', 'm2a', + 'm2', 'm2b', 'lambda2', 'l2', 't2', 's2', 'r2', 'k2', 'eta2', 'mns2', + '2sm2', 'm3', 'mk3', 's3', 'mn4', 'm4', 'ms4', 'mk4', 's4', 's5', 'm6', + 's6', 's7', 's8', 'm8', 'mks2', 'msqm', 'mtm', 'n4', 'eps2', 'z0'] # set function for astronomical longitudes - ASTRO5 = True if kwargs['corrections'] in ('GOT','FES') else False + ASTRO5 = True if kwargs['corrections'] in ('GOT', 'FES') else False # convert from Modified Julian Dates into Ephemeris Time - s, h, p, omega, pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'], + s, h, p, n, pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'], ASTRO5=ASTRO5) + # number of temporal values nt = len(np.atleast_1d(MJD)) # initial time conversions hour = 24.0*np.mod(MJD, 1) - t1 = 15.0*hour - t2 = 30.0*hour - + # convert from hours solar time into mean lunar time in degrees + tau = 15.0*hour - s + h + # variable for multiples of 90 degrees (Ray technical note 2017) + k = 90.0 + np.zeros((nt)) # degrees to radians dtr = np.pi/180.0 - # Determine equilibrium arguments - arg = np.zeros((nt, 60)) - arg[:,0] = h - pp # Sa - arg[:,1] = 2.0*h # Ssa - arg[:,2] = s - p # Mm - arg[:,3] = 2.0*s - 2.0*h # MSf - arg[:,4] = 2.0*s # Mf - arg[:,5] = 3.0*s - p # Mt - arg[:,6] = t1 - 5.0*s + 3.0*h + p - 90.0 # alpha1 - arg[:,7] = t1 - 4.0*s + h + 2.0*p - 90.0 # 2Q1 - arg[:,8] = t1 - 4.0*s + 3.0*h - 90.0 # sigma1 - arg[:,9] = t1 - 3.0*s + h + p - 90.0 # q1 - arg[:,10] = t1 - 3.0*s + 3.0*h - p - 90.0 # rho1 - arg[:,11] = t1 - 2.0*s + h - 90.0 # o1 - arg[:,12] = t1 - 2.0*s + 3.0*h + 90.0 # tau1 - arg[:,13] = t1 - s + h + 90.0 # M1 - arg[:,14] = t1 - s + 3.0*h - p + 90.0 # chi1 - arg[:,15] = t1 - 2.0*h + pp - 90.0 # pi1 - arg[:,16] = t1 - h - 90.0 # p1 - if kwargs['corrections'] in ('OTIS','ATLAS','TMD3','netcdf'): - arg[:,17] = t1 + 90.0 # s1 - elif kwargs['corrections'] in ('GOT','FES'): - arg[:,17] = t1 + 180.0 # s1 (Doodson's phase) - arg[:,18] = t1 + h + 90.0 # k1 - arg[:,19] = t1 + 2.0*h - pp + 90.0 # psi1 - arg[:,20] = t1 + 3.0*h + 90.0 # phi1 - arg[:,21] = t1 + s - h + p + 90.0 # theta1 - arg[:,22] = t1 + s + h - p + 90.0 # J1 - arg[:,23] = t1 + 2.0*s + h + 90.0 # OO1 - arg[:,24] = t2 - 4.0*s + 2.0*h + 2.0*p # 2N2 - arg[:,25] = t2 - 4.0*s + 4.0*h # mu2 - arg[:,26] = t2 - 3.0*s + 2.0*h + p # n2 - arg[:,27] = t2 - 3.0*s + 4.0*h - p # nu2 - arg[:,28] = t2 - 2.0*s + h + pp # M2a - arg[:,29] = t2 - 2.0*s + 2.0*h # M2 - arg[:,30] = t2 - 2.0*s + 3.0*h - pp # M2b - arg[:,31] = t2 - s + p + 180.0 # lambda2 - arg[:,32] = t2 - s + 2.0*h - p + 180.0 # L2 - arg[:,33] = t2 - h + pp # T2 - arg[:,34] = t2 # S2 - arg[:,35] = t2 + h - pp + 180.0 # R2 - arg[:,36] = t2 + 2.0*h # K2 - arg[:,37] = t2 + s + 2.0*h - pp # eta2 - arg[:,38] = t2 - 5.0*s + 4.0*h + p # MNS2 - arg[:,39] = t2 + 2.0*s - 2.0*h # 2SM2 - arg[:,40] = 1.5*arg[:,29] # M3 - arg[:,41] = arg[:,18] + arg[:,29] # MK3 - arg[:,42] = 3.0*t1 # S3 - arg[:,43] = arg[:,26] + arg[:,29] # MN4 - arg[:,44] = 2.0*arg[:,29] # M4 - arg[:,45] = arg[:,29] + arg[:,34] # MS4 - arg[:,46] = arg[:,29] + arg[:,36] # MK4 - arg[:,47] = 4.0*t1 # S4 - arg[:,48] = 5.0*t1 # S5 - arg[:,49] = 3.0*arg[:,29] # M6 - arg[:,50] = 3.0*t2 # S6 - arg[:,51] = 7.0*t1 # S7 - arg[:,52] = 4.0*t2 # S8 - # shallow water constituents - arg[:,53] = 4.0*arg[:,29] # m8 - arg[:,54] = arg[:,29] + arg[:,36] - arg[:,34] # mks2 - arg[:,55] = 4.0*s - 2.0*h # msqm - arg[:,56] = 3.0*s - p # mtm - arg[:,57] = 2.0*arg[:,26] # n4 - arg[:,58] = t2 - 5.0*s + 4.0*h + p # eps2 - # mean sea level - arg[:,59] = 0.0 # Z0 + # determine equilibrium arguments + fargs = np.c_[tau, s, h, p, n, pp, k] + arg = np.dot(fargs, _arguments_table(**kwargs)) # determine nodal corrections f and u - sinn = np.sin(omega*dtr) - cosn = np.cos(omega*dtr) - sin2n = np.sin(2.0*omega*dtr) - cos2n = np.cos(2.0*omega*dtr) - sin3n = np.sin(3.0*omega*dtr) + sinn = np.sin(n*dtr) + cosn = np.cos(n*dtr) + sin2n = np.sin(2.0*n*dtr) + cos2n = np.cos(2.0*n*dtr) + sin3n = np.sin(3.0*n*dtr) # set nodal corrections + # scale factor correction f = np.zeros((nt, 60)) + # phase correction u = np.zeros((nt, 60)) # determine nodal corrections f and u for each model type - if kwargs['corrections'] in ('OTIS','ATLAS','TMD3','netcdf'): + if kwargs['corrections'] in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): f[:,0] = 1.0 # Sa f[:,1] = 1.0 # Ssa f[:,2] = 1.0 - 0.130*cosn # Mm @@ -231,12 +183,14 @@ def arguments( temp2 = (0.189*sinn - 0.0058*sin2n)**2 f[:,11] = np.sqrt(temp1 + temp2) # O1 f[:,12] = 1.0 # tau1 - # Doodson's - # Mtmp1 = 2.0*np.cos(p*dtr) + 0.4*np.cos((p-omega)*dtr) - # Mtmp2 = np.sin(p*dtr) + 0.2*np.sin((p-omega)*dtr) - # Ray's - Mtmp1 = 1.36*np.cos(p*dtr) + 0.267*np.cos((p-omega)*dtr) - Mtmp2 = 0.64*np.sin(p*dtr) + 0.135*np.sin((p-omega)*dtr) + if (kwargs['M1'] == 'Doodson'): + # A. T. Doodson's coefficients for M1 tides + Mtmp1 = 2.0*np.cos(p*dtr) + 0.4*np.cos((p-n)*dtr) + Mtmp2 = np.sin(p*dtr) + 0.2*np.sin((p-n)*dtr) + elif (kwargs['M1'] == 'Ray'): + # R. Ray's coefficients for M1 tides + Mtmp1 = 1.36*np.cos(p*dtr) + 0.267*np.cos((p-n)*dtr) + Mtmp2 = 0.64*np.sin(p*dtr) + 0.135*np.sin((p-n)*dtr) f[:,13] = np.sqrt(Mtmp1**2 + Mtmp2**2) # M1 f[:,14] = np.sqrt((1.0+0.221*cosn)**2+(0.221*sinn)**2) # chi1 f[:,15] = 1.0 # pi1 @@ -262,8 +216,10 @@ def arguments( f[:,29] = f[:,24] # M2 f[:,30] = 1.0 # M2b f[:,31] = 1.0 # lambda2 - Ltmp1 = 1.0 - 0.25*np.cos(2*p*dtr) - 0.11*np.cos((2.0*p-omega)*dtr) - 0.04*cosn - Ltmp2 = 0.25*np.sin(2*p*dtr) + 0.11*np.sin((2.0*p-omega)*dtr) + 0.04*sinn + Ltmp1 = 1.0 - 0.25*np.cos(2*p*dtr) - \ + 0.11*np.cos((2.0*p - n)*dtr) - 0.04*cosn + Ltmp2 = 0.25*np.sin(2*p*dtr) + \ + 0.11*np.sin((2.0*p - n)*dtr) + 0.04*sinn f[:,32] = np.sqrt(Ltmp1**2 + Ltmp2**2) # L2 f[:,33] = 1.0 # T2 f[:,34] = 1.0 # S2 @@ -365,10 +321,10 @@ def arguments( elif kwargs['corrections'] in ('FES',): # additional astronomical terms for FES models - II = np.arccos(0.913694997 - 0.035692561*np.cos(omega*dtr)) - at1 = np.arctan(1.01883*np.tan(omega*dtr/2.0)) - at2 = np.arctan(0.64412*np.tan(omega*dtr/2.0)) - xi = -at1 - at2 + omega*dtr + II = np.arccos(0.913694997 - 0.035692561*np.cos(n*dtr)) + at1 = np.arctan(1.01883*np.tan(n*dtr/2.0)) + at2 = np.arctan(0.64412*np.tan(n*dtr/2.0)) + xi = -at1 - at2 + n*dtr xi[xi > np.pi] -= 2.0*np.pi nu = at1 - at2 I2 = np.tan(II/2.0) @@ -393,9 +349,14 @@ def arguments( f[:,9] = f[:,7] # q1 f[:,10] = f[:,7] # rho1 f[:,11] = f[:,7] # O1 - # Ray's - Mtmp1 = 1.36*np.cos(p*dtr) + 0.267*np.cos((p-omega)*dtr) - Mtmp2 = 0.64*np.sin(p*dtr) + 0.135*np.sin((p-omega)*dtr) + if (kwargs['M1'] == 'Doodson'): + # A. T. Doodson's coefficients for M1 tides + Mtmp1 = 2.0*np.cos(p*dtr) + 0.4*np.cos((p-n)*dtr) + Mtmp2 = np.sin(p*dtr) + 0.2*np.sin((p-n)*dtr) + elif (kwargs['M1'] == 'Ray'): + # R. Ray's coefficients for M1 tides + Mtmp1 = 1.36*np.cos(p*dtr) + 0.267*np.cos((p-n)*dtr) + Mtmp2 = 0.64*np.sin(p*dtr) + 0.135*np.sin((p-n)*dtr) f[:,13] = np.sqrt(Mtmp1**2 + Mtmp2**2) # M1 f[:,14] = np.sin(2.0*II) / 0.7214 # chi1 f[:,15] = 1.0 # pi1 @@ -533,15 +494,459 @@ def arguments( # take pu,pf,G for the set of given constituents nc = len(constituents) + # scale factor correction pu = np.zeros((nt,nc)) + # phase correction pf = np.zeros((nt,nc)) + # equilibrium arguments G = np.zeros((nt,nc)) for i,c in enumerate(constituents): # map between given constituents and supported in tidal program - j, = [j for j,val in enumerate(cindex) if (val == c)] + assert c.lower() in cindex, f'Unsupported constituent {c.lower()}' + j, = [j for j,val in enumerate(cindex) if (val == c.lower())] pu[:,i] = u[:,j]*dtr pf[:,i] = f[:,j] G[:,i] = arg[:,j] # return values as tuple return (pu, pf, G) + +def minor_arguments( + MJD: np.ndarray, + **kwargs + ): + """ + Calculates the nodal corrections for minor tidal constituents + in order to infer their values [1]_ [2]_ [3]_ [4]_ + + Parameters + ---------- + MJD: np.ndarray + modified Julian day of input date + deltat: float or np.ndarray, default 0.0 + time correction for converting to Ephemeris Time (days) + corrections: str, default 'OTIS' + use nodal corrections from OTIS/ATLAS or GOT models + + Returns + ------- + pu: np.ndarray + nodal correction for the constituent amplitude + pf: np.ndarray + nodal correction for the constituent phase + G: np.ndarray + phase correction in degrees + + References + ---------- + .. [1] A. T. Doodson and H. D. Warburg, "Admiralty Manual of Tides", + HMSO, London, (1941). + .. [2] P. Schureman, "Manual of Harmonic Analysis and Prediction of Tides," + *US Coast and Geodetic Survey*, Special Publication, 98, (1958). + .. [3] M. G. G. Foreman and R. F. Henry, "The harmonic analysis of tidal + model time series," *Advances in Water Resources*, 12(3), 109--120, + (1989). `doi: 10.1016/0309-1708(89)90017-1 + `_ + .. [4] G. D. Egbert and S. Y. Erofeeva, "Efficient Inverse Modeling of + Barotropic Ocean Tides," *Journal of Atmospheric and Oceanic + Technology*, 19(2), 183--204, (2002). + `doi: 10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2`__ + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + # set default keyword arguments + kwargs.setdefault('deltat', 0.0) + kwargs.setdefault('corrections', 'OTIS') + + # degrees to radians + dtr = np.pi/180.0 + # set function for astronomical longitudes + ASTRO5 = True if kwargs['corrections'] in ('GOT', 'FES') else False + # convert from Modified Julian Dates into Ephemeris Time + s, h, p, n, pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'], + ASTRO5=ASTRO5) + + # number of temporal values + nt = len(np.atleast_1d(MJD)) + # initial time conversions + hour = 24.0*np.mod(MJD, 1) + # convert from hours solar time into mean lunar time in degrees + tau = 15.0*hour - s + h + # variable for multiples of 90 degrees (Ray technical note 2017) + k = 90.0 + np.zeros((nt)) + + # determine equilibrium arguments + fargs = np.c_[tau, s, h, p, n, pp, k] + arg = np.dot(fargs, _minor_table()) + + # determine nodal corrections f and u + sinn = np.sin(n*dtr) + cosn = np.cos(n*dtr) + sin2n = np.sin(2.0*n*dtr) + cos2n = np.cos(2.0*n*dtr) + + # scale factor corrections for minor constituents + f = np.ones((nt, 20)) + f[:,0] = np.sqrt((1.0 + 0.189*cosn - 0.0058*cos2n)**2 + + (0.189*sinn - 0.0058*sin2n)**2)# 2Q1 + f[:,1] = f[:,0]# sigma1 + f[:,2] = f[:,0]# rho1 + f[:,3] = np.sqrt((1.0 + 0.185*cosn)**2 + (0.185*sinn)**2)# M12 + f[:,4] = np.sqrt((1.0 + 0.201*cosn)**2 + (0.201*sinn)**2)# M11 + f[:,5] = np.sqrt((1.0 + 0.221*cosn)**2 + (0.221*sinn)**2)# chi1 + f[:,9] = np.sqrt((1.0 + 0.198*cosn)**2 + (0.198*sinn)**2)# J1 + f[:,10] = np.sqrt((1.0 + 0.640*cosn + 0.134*cos2n)**2 + + (0.640*sinn + 0.134*sin2n)**2)# OO1 + f[:,11] = np.sqrt((1.0 - 0.0373*cosn)**2 + (0.0373*sinn)**2)# 2N2 + f[:,12] = f[:,11]# mu2 + f[:,13] = f[:,11]# nu2 + f[:,15] = f[:,11]# L2 + f[:,16] = np.sqrt((1.0 + 0.441*cosn)**2 + (0.441*sinn)**2)# L2 + + # phase corrections for minor constituents + u = np.zeros((nt, 20)) + u[:,0] = np.arctan2(0.189*sinn - 0.0058*sin2n, + 1.0 + 0.189*cosn - 0.0058*sin2n)/dtr# 2Q1 + u[:,1] = u[:,0]# sigma1 + u[:,2] = u[:,0]# rho1 + u[:,3] = np.arctan2( 0.185*sinn, 1.0 + 0.185*cosn)/dtr# M12 + u[:,4] = np.arctan2(-0.201*sinn, 1.0 + 0.201*cosn)/dtr# M11 + u[:,5] = np.arctan2(-0.221*sinn, 1.0 + 0.221*cosn)/dtr# chi1 + u[:,9] = np.arctan2(-0.198*sinn, 1.0 + 0.198*cosn)/dtr# J1 + u[:,10] = np.arctan2(-0.640*sinn - 0.134*sin2n, + 1.0 + 0.640*cosn + 0.134*cos2n)/dtr# OO1 + u[:,11] = np.arctan2(-0.0373*sinn, 1.0 - 0.0373*cosn)/dtr# 2N2 + u[:,12] = u[:,11]# mu2 + u[:,13] = u[:,11]# nu2 + u[:,15] = u[:,11]# L2 + u[:,16] = np.arctan2(-0.441*sinn, 1.0 + 0.441*cosn)/dtr# L2 + + if kwargs['corrections'] in ('FES',): + # additional astronomical terms for FES models + II = np.arccos(0.913694997 - 0.035692561*np.cos(n*dtr)) + at1 = np.arctan(1.01883*np.tan(n*dtr/2.0)) + at2 = np.arctan(0.64412*np.tan(n*dtr/2.0)) + xi = -at1 - at2 + n*dtr + xi[xi > np.pi] -= 2.0*np.pi + nu = at1 - at2 + I2 = np.tan(II/2.0) + Ra1 = np.sqrt(1.0 - 12.0*(I2**2)*np.cos(2.0*(p - xi)) + 36.0*(I2**4)) + P2 = np.sin(2.0*(p - xi)) + Q2 = 1.0/(6.0*(I2**2)) - np.cos(2.0*(p - xi)) + R = np.arctan(P2/Q2) + + # scale factor corrections for minor constituents + f[:,0] = np.sin(II)*(np.cos(II/2.0)**2)/0.38 # 2Q1 + f[:,1] = f[:,0] # sigma1 + f[:,2] = f[:,0] # rho1 + f[:,3] = f[:,0] # M12 + f[:,4] = np.sin(2.0*II)/0.7214 # M11 + f[:,5] = f[:,4] # chi1 + f[:,9] = f[:,5] # J1 + f[:,10] = np.sin(II)*np.power(np.sin(II/2.0),2.0)/0.01640 # OO1 + f[:,11] = np.power(np.cos(II/2.0),4.0)/0.9154 # 2N2 + f[:,12] = f[:,11] # mu2 + f[:,13] = f[:,11] # nu2 + f[:,14] = f[:,11] # lambda2 + f[:,15] = f[:,11]*Ra1 # L2 + f[:,18] = f[:,11] # eps2 + f[:,19] = np.power(np.sin(II),2.0)/0.1565 # eta2 + + # phase corrections for minor constituents + u[:,0] = (2.0*xi - nu)/dtr # 2Q1 + u[:,1] = u[:,0] # sigma1 + u[:,2] = u[:,0] # rho1 + u[:,3] = u[:,0] # M12 + u[:,4] = -nu/dtr # M11 + u[:,5] = u[:,4] # chi1 + u[:,9] = u[:,4] # J1 + u[:,10] = (-2.0*xi - nu)/dtr # OO1 + u[:,11] = (2.0*xi - 2.0*nu)/dtr # 2N2 + u[:,12] = u[:,11] # mu2 + u[:,13] = u[:,11] # nu2 + u[:,14] = (2.0*xi - 2.0*nu)/dtr # lambda2 + u[:,15] = (2.0*xi - 2.0*nu - R)/dtr# L2 + u[:,18] = u[:,12] # eps2 + u[:,19] = -2.0*nu/dtr # eta2 + + # return values as tuple + return (u, f, arg) + +def doodson_number( + constituents: str | list | np.ndarray, + **kwargs + ): + """ + Calculates the Doodson or Cartwright number for + tidal constituents [1]_ + + Parameters + ---------- + constituents: str, list or np.ndarray + tidal constituent ID(s) + corrections: str, default 'OTIS' + use arguments from OTIS/ATLAS or GOT models + formalism: str, default 'Doodson' + constituent identifier formalism + + - ``'Cartwright'`` + - ``'Doodson'`` + + Returns + ------- + numbers: float, np.ndarray or dict + Doodson or Cartwright number for each constituent + + References + ---------- + .. [1] A. T. Doodson and H. Lamb, "The harmonic development of + the tide-generating potential", *Proceedings of the Royal Society + of London. Series A, Containing Papers of a Mathematical and + Physical Character*, 100(704), 305--329, (1921). + `doi: 10.1098/rspa.1921.0088 `_ + """ + # set default keyword arguments + kwargs.setdefault('corrections', 'OTIS') + kwargs.setdefault('formalism', 'Doodson') + # validate inputs + assert kwargs['formalism'].title() in ('Cartwright', 'Doodson'), \ + f'Unknown formalism {kwargs["formalism"]}' + + # constituents array (not all are included in tidal program) + cindex = ['sa', 'ssa', 'mm', 'msf', 'mf', 'mt', 'alpha1', '2q1', 'sigma1', + 'q1', 'rho1', 'o1', 'tau1', 'm1', 'chi1', 'pi1', 'p1', 's1', 'k1', + 'psi1', 'phi1', 'theta1', 'j1', 'oo1', '2n2', 'mu2', 'n2', 'nu2', 'm2a', + 'm2', 'm2b', 'lambda2', 'l2', 't2', 's2', 'r2', 'k2', 'eta2', 'mns2', + '2sm2', 'm3', 'mk3', 's3', 'mn4', 'm4', 'ms4', 'mk4', 's4', 's5', 'm6', + 's6', 's7', 's8', 'm8', 'mks2', 'msqm', 'mtm', 'n4', 'eps2', 'z0'] + # get the table of coefficients + coefficients = _arguments_table(**kwargs) + if isinstance(constituents, str): + # map between given constituents and supported in tidal program + assert constituents.lower() in cindex, \ + f'Unsupported constituent {constituents}' + j, = [j for j,val in enumerate(cindex) if (val == constituents.lower())] + # extract identifier in formalism + if (kwargs['formalism'] == 'Cartwright'): + # extract Cartwright number + numbers = np.array(coefficients[:6,j]) + elif (kwargs['formalism'] == 'Doodson'): + # convert from coefficients to Doodson number + numbers = _to_doodson_number(coefficients[:,j]) + else: + # output dictionary with Doodson numbers + numbers = {} + # for each input constituent + for i,c in enumerate(constituents): + # map between given constituents and supported in tidal program + assert c.lower() in cindex, f'Unsupported constituent {c}' + j, = [j for j,val in enumerate(cindex) if (val == c.lower())] + # convert from coefficients to Doodson number + if (kwargs['formalism'] == 'Cartwright'): + # extract Cartwright number + numbers[c] = np.array(coefficients[:6,j]) + elif (kwargs['formalism'] == 'Doodson'): + # convert from coefficients to Doodson number + numbers[c] = _to_doodson_number(coefficients[:,j]) + # return the Doodson or Cartwright number + return numbers + +def _arguments_table(**kwargs): + """ + Arguments table for tidal constituents [1]_ [2]_ + + Parameters + ---------- + corrections: str, default 'OTIS' + use arguments from OTIS/ATLAS or GOT models + + Returns + ------- + coef: np.ndarray + Doodson coefficients (Cartwright numbers) for each constituent + + References + ---------- + .. [1] A. T. Doodson and H. Lamb, "The harmonic development of + the tide-generating potential", *Proceedings of the Royal Society + of London. Series A, Containing Papers of a Mathematical and + Physical Character*, 100(704), 305--329, (1921). + `doi: 10.1098/rspa.1921.0088 `_ + .. [2] A. T. Doodson and H. D. Warburg, "Admiralty Manual of Tides", + HMSO, London, (1941). + """ + # set default keyword arguments + kwargs.setdefault('corrections', 'OTIS') + + # modified Doodson coefficients for constituents + # using 7 index variables: tau, s, h, p, n, pp, k + # tau: mean lunar time + # s: mean longitude of moon + # h: mean longitude of sun + # p: mean longitude of lunar perigee + # n: mean longitude of ascending lunar node + # pp: mean longitude of solar perigee + # k: 90-degree phase + coef = np.zeros((7, 60)) + coef[:,0] = [0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0] # Sa + coef[:,1] = [0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0] # Ssa + coef[:,2] = [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0] # Mm + coef[:,3] = [0.0, 2.0, -2.0, 0.0, 0.0, 0.0, 0.0] # MSf + coef[:,4] = [0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0] # Mf + coef[:,5] = [0.0, 3.0, 0.0, -1.0, 0.0, 0.0, 0.0] # Mt + coef[:,6] = [1.0, -4.0, 2.0, 1.0, 0.0, 0.0, -1.0] # alpha1 + coef[:,7] = [1.0, -3.0, 0.0, 2.0, 0.0, 0.0, -1.0] # 2q1 + coef[:,8] = [1.0, -3.0, 2.0, 0.0, 0.0, 0.0, -1.0] # sigma1 + coef[:,9] = [1.0, -2.0, 0.0, 1.0, 0.0, 0.0, -1.0]# q1 + coef[:,10] = [1.0, -2.0, 2.0, -1.0, 0.0, 0.0, -1.0] # rho1 + coef[:,11] = [1.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0] # o1 + coef[:,12] = [1.0, -1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # tau1 + coef[:,13] = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0] # m1 + coef[:,14] = [1.0, 0.0, 2.0, -1.0, 0.0, 0.0, 1.0] # chi1 + coef[:,15] = [1.0, 1.0, -3.0, 0.0, 0.0, 1.0, -1.0] # pi1 + coef[:,16] = [1.0, 1.0, -2.0, 0.0, 0.0, 0.0, -1.0] # p1 + if kwargs['corrections'] in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): + coef[:,17] = [1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0] # s1 + elif kwargs['corrections'] in ('GOT', 'FES'): + coef[:,17] = [1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 2.0] # s1 (Doodson's phase) + coef[:,18] = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0] # k1 + coef[:,19] = [1.0, 1.0, 1.0, 0.0, 0.0, -1.0, 1.0] # psi1 + coef[:,20] = [1.0, 1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # phi1 + coef[:,21] = [1.0, 2.0, -2.0, 1.0, 0.0, 0.0, 1.0] # theta1 + coef[:,22] = [1.0, 2.0, 0.0, -1.0, 0.0, 0.0, 1.0] # j1 + coef[:,23] = [1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0] # oo1 + coef[:,24] = [2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # 2n2 + coef[:,25] = [2.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0] # mu2 + coef[:,26] = [2.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0] # n2 + coef[:,27] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 0.0] # nu2 + coef[:,28] = [2.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0] # m2a + coef[:,29] = [2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m2 + coef[:,30] = [2.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0] # m2b + coef[:,31] = [2.0, 1.0, -2.0, 1.0, 0.0, 0.0, 2.0] # lambda2 + coef[:,32] = [2.0, 1.0, 0.0, -1.0, 0.0, 0.0, 2.0] # l2 + coef[:,33] = [2.0, 2.0, -3.0, 0.0, 0.0, 1.0, 0.0] # t2 + coef[:,34] = [2.0, 2.0, -2.0, 0.0, 0.0, 0.0, 0.0] # s2 + coef[:,35] = [2.0, 2.0, -1.0, 0.0, 0.0, -1.0, 2.0] # r2 + coef[:,36] = [2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0] # k2 + coef[:,37] = [2.0, 3.0, 0.0, 0.0, 0.0, -1.0, 0.0] # eta2 + coef[:,38] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # mns2 + coef[:,39] = [2.0, 4.0, -4.0, 0.0, 0.0, 0.0, 0.0] # 2sm2 + coef[:,40] = [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m3 + coef[:,41] = coef[:,18] + coef[:,29] # mk3 + coef[:,42] = [3.0, 3.0, -3.0, 0.0, 0.0, 0.0, 0.0] # s3 + coef[:,43] = coef[:,26] + coef[:,29] # mn4 + coef[:,44] = [4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m4 + coef[:,45] = coef[:,29] + coef[:,34] # ms4 + coef[:,46] = coef[:,29] + coef[:,36] # mk4 + coef[:,47] = [4.0, 4.0, -4.0, 0.0, 0.0, 0.0, 0.0] # s4 + coef[:,48] = [5.0, 5.0, -5.0, 0.0, 0.0, 0.0, 0.0] # s5 + coef[:,49] = [6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m6 + coef[:,50] = [6.0, 6.0, -6.0, 0.0, 0.0, 0.0, 0.0] # s6 + coef[:,51] = [7.0, 7.0, -7.0, 0.0, 0.0, 0.0, 0.0] # s7 + coef[:,52] = [8.0, 8.0, -8.0, 0.0, 0.0, 0.0, 0.0] # s8 + # shallow water constituents + coef[:,53] = [8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # m8 + coef[:,54] = coef[:,29] + coef[:,36] - coef[:,34] # mks2 + coef[:,55] = [0.0, 4.0, -2.0, 0.0, 0.0, 0.0, 0.0] # msqm + coef[:,56] = [0.0, 3.0, 0.0, -1.0, 0.0, 0.0, 0.0] # mtm + coef[:,57] = [4.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # n4 + coef[:,58] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # eps2 + # mean sea level + coef[:,59] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # z0 + # return the coefficient table + return coef + +def _minor_table(**kwargs): + """ + Arguments table for minor tidal constituents [1]_ [2]_ + + Returns + ------- + coef: np.ndarray + Doodson coefficients (Cartwright numbers) for each constituent + + References + ---------- + .. [1] A. T. Doodson and H. Lamb, "The harmonic development of + the tide-generating potential", *Proceedings of the Royal Society + of London. Series A, Containing Papers of a Mathematical and + Physical Character*, 100(704), 305--329, (1921). + `doi: 10.1098/rspa.1921.0088 `_ + .. [2] A. T. Doodson and H. D. Warburg, "Admiralty Manual of Tides", + HMSO, London, (1941). + """ + # modified Doodson coefficients for constituents + # using 7 index variables: tau, s, h, p, n, pp, k + # tau: mean lunar time + # s: mean longitude of moon + # h: mean longitude of sun + # p: mean longitude of lunar perigee + # n: mean longitude of ascending lunar node + # pp: mean longitude of solar perigee + # k: 90-degree phase + coef = np.zeros((7, 20)) + coef[:,0] = [1.0, -3.0, 0.0, 2.0, 0.0, 0.0, -1.0] # 2q1 + coef[:,1] = [1.0, -3.0, 2.0, 0.0, 0.0, 0.0, -1.0] # sigma1 + coef[:,2] = [1.0, -2.0, 2.0, -1.0, 0.0, 0.0, -1.0] # rho1 + coef[:,3] = [1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0] # m1 + coef[:,4] = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0] # m1 + coef[:,5] = [1.0, 0.0, 2.0, -1.0, 0.0, 0.0, 1.0] # chi1 + coef[:,6] = [1.0, 1.0, -3.0, 0.0, 0.0, 1.0, -1.0] # pi1 + coef[:,7] = [1.0, 1.0, 2.0, 0.0, 0.0, 0.0, 1.0] # phi1 + coef[:,8] = [1.0, 2.0, -2.0, 1.0, 0.0, 0.0, 1.0] # theta1 + coef[:,9] = [1.0, 2.0, 0.0, -1.0, 0.0, 0.0, 1.0] # j1 + coef[:,10] = [1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0] # oo1 + coef[:,11] = [2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 0.0] # 2n2 + coef[:,12] = [2.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0] # mu2 + coef[:,13] = [2.0, -1.0, 2.0, -1.0, 0.0, 0.0, 0.0] # nu2 + coef[:,14] = [2.0, 1.0, -2.0, 1.0, 0.0, 0.0, 2.0]# lambda2 + coef[:,15] = [2.0, 1.0, 0.0, -1.0, 0.0, 0.0, 2.0] # l2 + coef[:,16] = [2.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0] # l2 + coef[:,17] = [2.0, 2.0, -3.0, 0.0, 0.0, 1.0, 0.0] # t2 + coef[:,18] = [2.0, -3.0, 2.0, 1.0, 0.0, 0.0, 0.0] # eps2 + coef[:,19] = [2.0, 3.0, 0.0, 0.0, 0.0, -1.0, 0.0] # eta2 + # return the coefficient table + return coef + +def _to_doodson_number(coef: list | np.ndarray): + """ + Converts Cartwright numbers into a Doodson number + + Parameters + ---------- + coef: list or np.ndarray + Doodson coefficients (Cartwright numbers) for constituent + + Returns + ------- + DO: float + Doodson number for constituent + """ + # assert length and verify array + coef = np.array(coef[:6]) + # add 5 to values following Doodson convention (prevent negatives) + coef[1:] += 5 + # convert to single number and round off floating point errors + DO = np.sum([v*10**(2-o) for o,v in enumerate(coef)]) + return np.round(DO, decimals=3) + +def _from_doodson_number(DO: float | np.ndarray): + """ + Converts Doodson numbers into Cartwright numbers + + Parameters + ---------- + DO: float or np.ndarray + Doodson number for constituent + + Returns + ------- + coef: np.ndarray + Doodson coefficients (Cartwright numbers) for constituent + """ + # convert from Doodson number to Cartwright numbers + # multiply by 1000 to prevent floating point errors + coef = np.array([np.mod(1e3*DO, 10**(6-o))//10**(5-o) for o in range(6)]) + # remove 5 from values following Doodson convention + coef[1:] -= 5 + return coef diff --git a/pyTMD/astro.py b/pyTMD/astro.py index f611229d..10c3a17c 100644 --- a/pyTMD/astro.py +++ b/pyTMD/astro.py @@ -1,18 +1,9 @@ #!/usr/bin/env python u""" astro.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (12/2023) Astronomical and nutation routines -mean_longitudes is a modification of the ASTROL fortran -subroutine by Richard Ray written in 03/1999 - -Computes the basic astronomical mean longitudes: s, h, p, N and PP -Note N is not N', i.e. N is decreasing with time. - -Formulae for the period 1990--2010 were derived by David Cartwright -MEEUS and ASTRO5 formulae are from versions of Meeus's Astronomical Algorithms - PYTHON DEPENDENCIES: numpy: Scientific Computing Tools For Python https://numpy.org @@ -25,6 +16,8 @@ Oliver Montenbruck, Practical Ephemeris Calculations, 1989. UPDATE HISTORY: + Updated 12/2023: refactored phase_angles function to doodson_arguments + added option to compute mean lunar time using equinox method Updated 05/2023: add wrapper function for nutation angles download JPL kernel file if not currently existing Updated 04/2023: added low resolution solar and lunar positions @@ -52,6 +45,7 @@ import logging import pathlib +import warnings import numpy as np from pyTMD.time import timescale from pyTMD.eop import iers_polar_motion @@ -127,9 +121,9 @@ def mean_longitudes( MEEUS: bool = False, ASTRO5: bool = False ): - """ + r""" Computes the basic astronomical mean longitudes: - `S`, `H`, `P`, `N` and `PP` [1]_ + `S`, `H`, `P`, `N` and `PP` [1]_ [2]_ Note `N` is not `N'`, i.e. `N` is decreasing with time. @@ -158,6 +152,14 @@ def mean_longitudes( References ---------- .. [1] J. Meeus, *Astronomical Algorithms*, 2nd edition, 477 pp., (1998). + .. [2] J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touz\ |eacute|\, + G. Francou, and J. Laskar "Numerical expressions for precession + formulae and mean elements for the Moon and the planets", + *Astronomy and Astrophysics*, 282(2), 663--683, (1994). + `bibcode: 1994A%26A...282..663S + `_ + + .. |eacute| unicode:: U+00E9 .. LATIN SMALL LETTER E WITH ACUTE """ circle = 360.0 if MEEUS: @@ -215,7 +217,7 @@ def mean_longitudes( # mean longitude of ascending lunar node N = 125.0445 - 0.05295377 * T # solar perigee at epoch 2000 - PP = 282.8 + PP = np.full_like(T, 282.8) # take the modulus of each S = np.mod(S, circle) H = np.mod(H, circle) @@ -226,28 +228,48 @@ def mean_longitudes( # PURPOSE: computes the phase angles of astronomical means def phase_angles(MJD: np.ndarray): + # raise warning for deprecated function call + warnings.warn(("Deprecated. Please use " + "pyTMD.astro.doodson_arguments instead"), + DeprecationWarning) + # call updated function to not break current workflows + TAU, S, H, P, ZNS, PS = doodson_arguments(MJD) + # return as tuple + return (S, H, P, TAU, ZNS, PS) + +# PURPOSE: computes the phase angles of astronomical means +def doodson_arguments( + MJD: np.ndarray, + equinox: bool = False, + apply_correction: bool = True, + ): """ - Computes astronomical phase angles for the sun and moon: - `S`, `H`, `P`, `TAU`, `ZNS` and `PS` [1]_ [2]_ + Computes astronomical phase angles for the six Doodson + Arguments: `tau`, `S`, `H`, `P`, and `N'`, and `Ps` + [1]_ [2]_ Parameters ---------- MJD: np.ndarray Modified Julian Day (MJD) of input date + equinox: bool, default False + use equinox method for calculating mean lunar time + apply_correction: bool, default True + Apply correction for mean lunar longitude Returns ------- + TAU: np.ndarray + mean lunar time (radians) S: np.ndarray - mean longitude of moon (radians) + mean longitude of the moon (radians) H: np.ndarray - mean longitude of sun (radians) + mean longitude of the sun (radians) P: np.ndarray mean longitude of lunar perigee (radians) - TAU: np.ndarray - time angle in lunar days (radians) - ZNS: np.ndarray - mean longitude of ascending lunar node `N'` (radians) - PS: np.ndarray + Np: np.ndarray + negative mean longitude of the ascending node (radians) + Ps: np.ndarray mean longitude of solar perigee (radians) References @@ -266,39 +288,48 @@ def phase_angles(MJD: np.ndarray): # 360 degrees circle = 360.0 # hour of the day - FHR = np.mod(MJD, 1)*24.0 - # calculate phase angles + hour = np.mod(MJD, 1)*24.0 + # calculate Doodson phase angles # mean longitude of moon (degrees) S = polynomial_sum(np.array([218.3164477, 481267.88123421, -1.5786e-3, 1.855835e-6, -1.53388e-8]), T) - # time angle in lunar days (degrees) - TAU = ((FHR*15.0) - S + polynomial_sum(np.array([280.4606184, - 36000.7700536, 3.8793e-4, -2.58e-8]), T)) + # mean lunar time (degrees) + if equinox: + # create timescale from Modified Julian Day (MJD) + ts = timescale(MJD=MJD) + # use Greenwich Mean Sidereal Time (GMST) from the + # Equinox method converted to degrees + TAU = 360.0*ts.st + 180.0 - S + else: + TAU = ((hour*15.0) - S + polynomial_sum(np.array([280.4606184, + 36000.7700536, 3.8793e-4, -2.58e-8]), T)) # calculate correction for mean lunar longitude (degrees) - PR = polynomial_sum(np.array([0.0, 1.396971278, - 3.08889e-4, 2.1e-8, 7.0e-9]), T) - S += PR + if apply_correction: + PR = polynomial_sum(np.array([0.0, 1.396971278, + 3.08889e-4, 2.1e-8, 7.0e-9]), T) + S += PR # mean longitude of sun (degrees) H = polynomial_sum(np.array([280.46645, 36000.7697489, 3.0322222e-4, 2.0e-8, -6.54e-9]), T) # mean longitude of lunar perigee (degrees) P = polynomial_sum(np.array([83.3532465, 4069.0137287, -1.032172222e-2, -1.24991e-5, 5.263e-8]), T) - # mean longitude of ascending lunar node (degrees) - ZNS = polynomial_sum(np.array([234.95544499, 1934.13626197, + # negative of the mean longitude of the ascending node + # of the moon (degrees) + Np = polynomial_sum(np.array([234.95544499, 1934.13626197, -2.07561111e-3, -2.13944e-6, 1.65e-8]), T) # mean longitude of solar perigee (degrees) - PS = polynomial_sum(np.array([282.93734098, 1.71945766667, + Ps = polynomial_sum(np.array([282.93734098, 1.71945766667, 4.5688889e-4, -1.778e-8, -3.34e-9]), T) # take the modulus of each and convert to radians S = dtr*np.mod(S, circle) H = dtr*np.mod(H, circle) P = dtr*np.mod(P, circle) TAU = dtr*np.mod(TAU, circle) - ZNS = dtr*np.mod(ZNS, circle) - PS = dtr*np.mod(PS, circle) + Np = dtr*np.mod(Np, circle) + Ps = dtr*np.mod(Ps, circle) # return as tuple - return (S, H, P, TAU, ZNS, PS) + return (TAU, S, H, P, Np, Ps) def delaunay_arguments(MJD: np.ndarray): """ diff --git a/pyTMD/check_points.py b/pyTMD/check_points.py index c8d5843e..2e3ca399 100644 --- a/pyTMD/check_points.py +++ b/pyTMD/check_points.py @@ -43,7 +43,7 @@ https://pypi.org/project/pyproj/ PROGRAM DEPENDENCIES: - convert_crs.py: convert points to and from Coordinates Reference Systems + crs.py: Coordinate Reference System (CRS) routines io/model.py: retrieves tide model parameters for named tide models io/OTIS.py: extract tidal harmonic constants from OTIS tide models io/ATLAS.py: extract tidal harmonic constants from ATLAS netcdf models @@ -52,6 +52,7 @@ interpolate.py: interpolation routines for spatial data UPDATE HISTORY: + Updated 12/2023: use new crs class for coordinate reprojection Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 04/2023: using pathlib to define and expand paths Updated 03/2023: add basic variable typing to function inputs @@ -74,9 +75,9 @@ import numpy as np import scipy.interpolate +import pyTMD.crs import pyTMD.io import pyTMD.io.model -import pyTMD.convert_crs import pyTMD.interpolate # attempt imports @@ -149,13 +150,7 @@ def check_points(x: np.ndarray, y: np.ndarray, # input shape of data idim = np.shape(x) # converting x,y from input coordinate reference system - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # convert to latitude and longitude + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon, lat = transformer.transform( @@ -171,7 +166,7 @@ def check_points(x: np.ndarray, y: np.ndarray, mz = np.logical_not(mz) # adjust dimensions of input coordinates to be iterable # run wrapper function to convert coordinate systems of input lat/lon - X, Y = pyTMD.convert_crs(lon, lat, model.projection, 'F') + X, Y = pyTMD.crs().convert(lon, lat, model.projection, 'F') elif (model.format == 'netcdf'): # if reading a netCDF OTIS atlas solution xi, yi, hz = pyTMD.io.ATLAS.read_netcdf_grid( diff --git a/pyTMD/compute_tide_corrections.py b/pyTMD/compute_tide_corrections.py index 6ab7006f..af4420b2 100644 --- a/pyTMD/compute_tide_corrections.py +++ b/pyTMD/compute_tide_corrections.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tide_corrections.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (12/2023) Calculates tidal elevations for correcting elevation or imagery data Ocean and Load Tides @@ -49,7 +49,7 @@ utilities.py: download and management utilities for syncing files arguments.py: load the nodal corrections for tidal constituents astro.py: computes the basic astronomical mean longitudes - convert_crs.py: convert points to and from Coordinates Reference Systems + crs.py: Coordinate Reference System (CRS) routines predict.py: predict tide values using harmonic constants io/model.py: retrieves tide model parameters for named tide models io/OTIS.py: extract tidal harmonic constants from OTIS tide models @@ -59,6 +59,7 @@ interpolate.py: interpolation routines for spatial data UPDATE HISTORY: + Updated 12/2023: use new crs class for coordinate reprojection Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 05/2023: use timescale class for time conversion operations use defaults from eop module for pole tide and EOP files @@ -106,6 +107,7 @@ import numpy as np import scipy.interpolate import pyTMD.constants +import pyTMD.crs import pyTMD.eop import pyTMD.io import pyTMD.time @@ -289,13 +291,7 @@ def compute_tide_corrections( y = np.atleast_1d(y) # converting x,y from EPSG to latitude/longitude - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # output coordinate reference system + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon, lat = transformer.transform(x.flatten(), y.flatten()) @@ -313,7 +309,7 @@ def compute_tide_corrections( nt = len(timescale) # read tidal constants and interpolate to grid points - if model.format in ('OTIS','ATLAS','TMD3'): + if model.format in ('OTIS', 'ATLAS', 'TMD3'): amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file, model.model_file, model.projection, type=model.type, method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, @@ -451,13 +447,7 @@ def compute_LPET_corrections( y = np.atleast_1d(y) # converting x,y from EPSG to latitude/longitude - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # output coordinate reference system + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon, lat = transformer.transform(x.flatten(), y.flatten()) @@ -579,13 +569,7 @@ def compute_LPT_corrections( y = np.atleast_1d(y) # converting x,y from EPSG to latitude/longitude - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # output coordinate reference system + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon,lat = transformer.transform(x.flatten(), y.flatten()) @@ -764,13 +748,7 @@ def compute_OPT_corrections( y = np.atleast_1d(y) # converting x,y from EPSG to latitude/longitude - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # output coordinate reference system + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon,lat = transformer.transform(x.flatten(), y.flatten()) @@ -963,13 +941,7 @@ def compute_SET_corrections( y = np.atleast_1d(y) # converting x,y from EPSG to latitude/longitude - try: - # EPSG projection code string or int - crs1 = pyproj.CRS.from_epsg(int(EPSG)) - except (ValueError,pyproj.exceptions.CRSError): - # Projection SRS string - crs1 = pyproj.CRS.from_string(EPSG) - # output coordinate reference system + crs1 = pyTMD.crs().from_input(EPSG) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) lon, lat = transformer.transform(x.flatten(), y.flatten()) diff --git a/pyTMD/constants.py b/pyTMD/constants.py index 28ea6093..ee655603 100755 --- a/pyTMD/constants.py +++ b/pyTMD/constants.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" constants.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Gravitational and ellipsoidal parameters @@ -37,6 +37,7 @@ https://iers-conventions.obspm.fr/content/tn36.pdf UPDATE HISTORY: + Updated 12/2023: add string representation of constants object Updated 10/2023: add value for WMO standard gravity Updated 03/2023: add basic variable typing set ellipsoid name and output units as constants attributes @@ -101,12 +102,12 @@ def __init__(self, ellipsoid: str = 'WGS84', units: str = 'MKS'): assert self.units in _units # set parameters for ellipsoid - if self.name in ('CLK66','NAD27'): + if self.name in ('CLK66', 'NAD27'): # Clarke 1866 self.a_axis = 6378206.4# [m] semimajor axis of the ellipsoid self.flat = 1.0/294.9786982# flattening of the ellipsoid - elif self.name in ('GRS80','NAD83'): + elif self.name in ('GRS80', 'NAD83'): # Geodetic Reference System 1980 # North American Datum 1983 self.a_axis = 6378135.0# [m] semimajor axis of the ellipsoid @@ -174,7 +175,7 @@ def __init__(self, ellipsoid: str = 'WGS84', units: str = 'MKS'): self.flat = 1.0/298.25642# flattening of the ellipsoid # set default parameters if not listed as part of ellipsoid - if self.name not in ('GRS80','GRS67','NAD83','TOPEX','EGM96'): + if self.name not in ('GRS80', 'GRS67', 'NAD83', 'TOPEX', 'EGM96'): # for ellipsoids not listing the Geocentric Gravitational Constant self.GM = 3.986004418e14# [m^3/s^2] @@ -388,3 +389,11 @@ def rho_e(self) -> float: """Average density """ return self.GM/(self.G*self.volume) + + def __str__(self): + """String representation of the ``constants`` object + """ + properties = ['pyTMD.constants'] + properties.append(f" name: {self.name}") + properties.append(f" units: {self.units}") + return '\n'.join(properties) diff --git a/pyTMD/convert_crs.py b/pyTMD/convert_crs.py index 8936fec2..48c46506 100644 --- a/pyTMD/convert_crs.py +++ b/pyTMD/convert_crs.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" convert_crs.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (12/2023) Converts points to and from Coordinates Reference Systems (CRS) CALLING SEQUENCE: @@ -30,6 +30,7 @@ https://pyproj4.github.io/pyproj/ UPDATE HISTORY: + Updated 03/2023: deprecated in favor of crs class Updated 03/2023: add basic variable typing to function inputs renamed coordinate reference system conversion functions Updated 02/2023: use named exception before passing to custom @@ -47,14 +48,9 @@ """ from __future__ import annotations -import logging +import warnings import numpy as np - -# attempt imports -try: - import pyproj -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - logging.critical("pyproj not available") +import pyTMD.crs def convert_crs( i1: np.ndarray, @@ -89,32 +85,10 @@ def convert_crs( o2: np.ndarray Projected y-coordinates (``'F'``) or latitude (``'B``') """ - # python dictionary with named conversion functions - conversion_functions = {} - conversion_functions['3031'] = _EPSG3031 - conversion_functions['3413'] = _EPSG3413 - conversion_functions['CATS2008'] = _CATS2008 - conversion_functions['3976'] = _EPSG3976 - conversion_functions['PSNorth'] = _PSNorth - conversion_functions['4326'] = _EPSG4326 - # check that PROJ for conversion was entered correctly - # run named conversion program and return values - try: - o1, o2 = conversion_functions[PROJ](i1, i2, BF, EPSG=EPSG) - except KeyError as exc: - pass - else: - return (o1, o2) - # try changing the projection using a custom projection - # run custom conversion program and return values - try: - o1, o2 = _custom(i1, i2, PROJ, BF, EPSG=EPSG) - except Exception as exc: - pass - else: - return (o1, o2) - # projection not found or available - raise Exception(f'PROJ: {PROJ} conversion function not found') + warnings.warn("Deprecated. Please use pyTMD.crs().convert instead", + DeprecationWarning) + # call updated function to not break current workflows + return pyTMD.crs().convert(i1, i2, PROJ, BF, EPSG=EPSG) # PURPOSE: try to get the projection information def crs_from_input(PROJECTION: int | str): @@ -126,343 +100,7 @@ def crs_from_input(PROJECTION: int | str): PROJECTION: int or str Coordinate Reference System code """ - # EPSG projection code - try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError, pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) - except (ValueError, pyproj.exceptions.CRSError): - pass - else: - return crs - # no projection can be made - raise pyproj.exceptions.CRSError - -# wrapper function for models in EPSG 3031 (Antarctic Polar Stereographic) -def _EPSG3031( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts models in EPSG 3031 (Antarctic Polar Stereographic) - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting from input EPSG (default latitude/longitude) - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_user_input({'proj':'stere','lat_0':-90,'lat_ts':-71, - 'lon_0':0,'x_0':0.,'y_0':0.,'ellps': 'WGS84','datum': 'WGS84', - 'units':'km'}) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to Polar-Stereographic x/y - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - # convert Polar-Stereographic x/y to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) - -# wrapper function for models in EPSG 3413 (Sea Ice Polar Stereographic North) -def _EPSG3413( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts models in EPSG 3413 (Sea Ice Polar Stereographic North) - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting from input EPSG (default latitude/longitude) - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_user_input({'proj':'stere','lat_0':90,'lat_ts':70, - 'lon_0':-45,'x_0':0.,'y_0':0.,'ellps': 'WGS84','datum': 'WGS84', - 'units':'km'}) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to Polar-Stereographic x/y - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - # convert Polar-Stereographic x/y to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) - -# wrapper function for CATS2008 tide models -def _CATS2008( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts Circum-Antarctic Tidal Simulation models - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting from input EPSG (default latitude/longitude) - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_user_input({'proj':'stere','lat_0':-90,'lat_ts':-71, - 'lon_0':-70,'x_0':0.,'y_0':0.,'ellps': 'WGS84','datum': 'WGS84', - 'units':'km'}) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to Polar-Stereographic x/y - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - # convert Polar-Stereographic x/y to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) - -# wrapper function for models in EPSG 3976 (NSIDC Sea Ice Stereographic South) -def _EPSG3976( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts models in EPSG 3976 (Sea Ice Polar Stereographic South) - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting from input EPSG (default latitude/longitude) - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_user_input({'proj':'stere','lat_0':-90,'lat_ts':-70, - 'lon_0':0,'x_0':0.,'y_0':0.,'ellps': 'WGS84','datum': 'WGS84', - 'units':'km'}) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to Polar-Stereographic x/y - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - # convert Polar-Stereographic x/y to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) - -# wrapper function for models in (idealized) PSNorth projection -def _PSNorth( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts idealized Arctic Polar Stereographic models - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting to and from input EPSG - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to (idealized) Polar-Stereographic x/y - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - lon, lat = transformer.transform(i1, i2, direction=direction) - o1 = (90.0-lat)*111.7*np.cos(lon/180.0*np.pi) - o2 = (90.0-lat)*111.7*np.sin(lon/180.0*np.pi) - # convert (idealized) Polar-Stereographic x/y to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - lon = np.arctan2(i2, i1)*180.0/np.pi - lat = 90.0 - np.sqrt(i1**2+i2**2)/111.7 - ii, = np.nonzero(lon < 0) - lon[ii] += 360.0 - o1, o2 = transformer.transform(lon, lat, direction=direction) - # return the output variables - return (o1, o2) - -# wrapper function to pass lat/lon values or convert if EPSG -def _EPSG4326( - i1: np.ndarray, - i2: np.ndarray, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts models in EPSG 4326 (WGS84 Latitude/Longitude) - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - crs1 = crs_from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) - -# wrapper function for using custom projections -def _custom( - i1: np.ndarray, - i2: np.ndarray, - PROJ: int | str, - BF: str, - EPSG: int | str = 4326 - ): - """ - Converts models in a custom projection - - Parameters - ---------- - i1: np.ndarray - Longitude (``'F'``) or projected x-coordinates (``'B'``) - i2: np.ndarray - Latitude (``'F'``) or projected y-coordinates (``'B'``) - PROJ: int or str - Spatial reference system code for coordinate transformations - BF: str - Direction of translation - - - ``'B'``: backwards - - ``'F'``: forwards - EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) - EPSG code for input (``'F'``) or output (``'B'``) coordinate system - - Returns - ------- - o1: np.ndarray - Projected x-coordinates (``'F'``) or longitude (``'B'``) - o2: np.ndarray - Projected y-coordinates (``'F'``) or latitude (``'B``') - """ - # projections for converting from input EPSG (default latitude/longitude) - crs1 = crs_from_input(EPSG) - crs2 = crs_from_input(PROJ) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # convert lat/lon to custom projection - if (BF.upper() == 'F'): - direction = pyproj.enums.TransformDirection.FORWARD - # convert custom projection to lat/lon - elif (BF.upper() == 'B'): - direction = pyproj.enums.TransformDirection.INVERSE - # return the output variables - return transformer.transform(i1, i2, direction=direction) + warnings.warn("Deprecated. Please use pyTMD.crs().from_input instead", + DeprecationWarning) + # call updated function to not break current workflows + return pyTMD.crs().from_input(PROJECTION) diff --git a/pyTMD/convert_ll_xy.py b/pyTMD/convert_ll_xy.py index 7c5b17f7..7f25e9b1 100644 --- a/pyTMD/convert_ll_xy.py +++ b/pyTMD/convert_ll_xy.py @@ -30,6 +30,7 @@ https://pyproj4.github.io/pyproj/ UPDATE HISTORY: + Updated 12/2023: updated function call to use crs class Updated 03/2023: deprecated in favor of convert_crs.py Updated 02/2023: use named exception before passing to custom Updated 11/2022: place some imports within try/except statements @@ -74,7 +75,7 @@ def convert_ll_xy(*args, **kwargs): o2: float Projected y-coordinates (``'F'``) or latitude (``'B``') """ - warnings.warn("Deprecated. Please use pyTMD.convert_crs instead", + warnings.warn("Deprecated. Please use pyTMD.crs().convert instead", DeprecationWarning) # call updated function to not break current workflows - return pyTMD.convert_crs(*args, **kwargs) + return pyTMD.crs().convert(*args, **kwargs) diff --git a/pyTMD/crs.py b/pyTMD/crs.py new file mode 100644 index 00000000..0d7924f6 --- /dev/null +++ b/pyTMD/crs.py @@ -0,0 +1,355 @@ +#!/usr/bin/env python +u""" +crs.py +Written by Tyler Sutterley (12/2023) +Coordinates Reference System (CRS) routines + +CALLING SEQUENCE: + x, y = pyTMD.crs().convert(lon, lat, PROJ, 'F') + lon, lat = pyTMD.crs().convert(x, y, PROJ, 'B') + +INPUTS: + i1: longitude ('F') or projection easting x ('B') + i2: latitude ('F') or projection northing y ('B') + PROJ: spatial reference system code for coordinate transformations + BF: backwards ('B') or forward ('F') translations + +OPTIONS: + EPSG: spatial reference system code for input (F) and output (B) coordinates + +OUTPUTS: + o1: projection easting x ('F') or longitude ('B') + o2: projection northing y ('F') or latitude ('B') + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python + https://numpy.org + https://numpy.org/doc/stable/user/numpy-for-matlab-users.html + pyproj: Python interface to PROJ library + https://pypi.org/project/pyproj/ + https://pyproj4.github.io/pyproj/ + +UPDATE HISTORY: + Updated 12/2023: converted conversion functions to class + Updated 03/2023: add basic variable typing to function inputs + renamed coordinate reference system conversion functions + Updated 02/2023: use named exception before passing to custom + Updated 11/2022: place some imports within try/except statements + use f-strings for formatting verbose or ascii output + Updated 04/2022: updated docstrings to numpy documentation format + Updated 09/2021: added function for using custom projections + Updated 06/2021: added 3413 for new 1km Greenland model from ESR + Updated 08/2020: using conversion protocols following pyproj-2 updates + https://pyproj4.github.io/pyproj/stable/gotchas.html + Updated 07/2020: added function docstrings. changed function name + Updated 03/2020: remove commented coordinate conversion functions + Updated 11/2019: using pyproj for coordinate conversions + Written 09/2017 +""" + +from __future__ import annotations + +import logging +import numpy as np + +# attempt imports +try: + import pyproj +except (AttributeError, ImportError, ModuleNotFoundError) as exc: + logging.critical("pyproj not available") + +class crs: + """Coordinate Reference System transformations for tide models + + Attributes + ---------- + name: str + Projection name + transformer: obj + ``pyproj`` transformer for changing coordinate reference system + """ + def __init__(self): + self.name = None + self.transformer = None + self._direction = None + + def convert(self, + i1: np.ndarray, + i2: np.ndarray, + PROJ: str, + BF: str, + EPSG: int | str = 4326 + ): + """ + Converts points to and from Coordinates Reference Systems (CRS) + + Parameters + ---------- + i1: np.ndarray + Input x-coordinates + i2: np.ndarray + Input y-coordinates + PROJ: str + Spatial reference system code for coordinate transformations + BF: str + Direction of transformation + + - ``'B'``: backwards + - ``'F'``: forwards + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + + Returns + ------- + o1: np.ndarray + Output transformed x-coordinates + o2: np.ndarray + Output transformed y-coordinates + """ + self.name = PROJ + # python dictionary with named conversion functions + transforms = {} + transforms['3031'] = self._EPSG3031 + transforms['3413'] = self._EPSG3413 + transforms['CATS2008'] = self._CATS2008 + transforms['3976'] = self._EPSG3976 + transforms['PSNorth'] = self._PSNorth + transforms['4326'] = self._EPSG4326 + # set the direction of the transform + self._direction = BF.upper() + # check that PROJ for conversion was entered correctly + # run named conversion program and return values + try: + transforms[PROJ](EPSG) + except KeyError as exc: + pass + else: + # return the output variables + return self.transform(i1, i2) + # try changing the projection using a custom projection + # run custom conversion program and return values + try: + self._custom(PROJ, EPSG=EPSG) + except Exception as exc: + pass + else: + return self.transform(i1, i2) + # projection not found or available + raise Exception(f'PROJ: {PROJ} conversion function not found') + + def transform(self, i1: np.ndarray, i2: np.ndarray): + """ + Performs Coordinates Reference System (CRS) transformations + + Parameters + ---------- + i1: np.ndarray + Input x-coordinates + i2: np.ndarray + Input y-coordinates + + Returns + ------- + o1: np.ndarray + Output transformed x-coordinates + o2: np.ndarray + Output transformed y-coordinates + """ + if (self.name == 'PSNorth') and (self.direction.name == 'FORWARD'): + # convert input coordinate reference system to lat/lon + lon, lat = self.transformer.transform(i1, i2, + direction=self.direction) + # convert lat/lon to (idealized) Polar-Stereographic x/y + o1 = (90.0 - lat)*111.7*np.cos(lon/180.0*np.pi) + o2 = (90.0 - lat)*111.7*np.sin(lon/180.0*np.pi) + elif (self.name == 'PSNorth') and (self.direction.name == 'INVERSE'): + # convert (idealized) Polar-Stereographic x/y to lat/lon + lon = np.arctan2(i2, i1)*180.0/np.pi + lat = 90.0 - np.sqrt(i1**2 + i2**2)/111.7 + # adjust longitudes to be -180:180 + ii, = np.nonzero(lon < 0) + lon[ii] += 360.0 + # convert to output coordinate reference system + o1, o2 = self.transformer.transform(lon, lat, + direction=self.direction) + else: + # convert coordinate reference system + o1, o2 = self.transformer.transform(i1, i2, + direction=self.direction) + # return the transformed coordinates + return (o1, o2) + + # PURPOSE: try to get the projection information + def from_input(self, PROJECTION: int | str): + """ + Attempt to retrieve the Coordinate Reference System + for an input code + + Parameters + ---------- + PROJECTION: int or str + Coordinate Reference System code + """ + # EPSG projection code + try: + CRS = pyproj.CRS.from_epsg(int(PROJECTION)) + except (ValueError, pyproj.exceptions.CRSError): + pass + else: + return CRS + # coordinate reference system string + try: + CRS = pyproj.CRS.from_string(PROJECTION) + except (ValueError, pyproj.exceptions.CRSError): + pass + else: + return CRS + # no projection can be made + raise pyproj.exceptions.CRSError + + def _EPSG3031(self, EPSG: int | str = 4326): + """ + Transform for models in EPSG:3031 (Antarctic Polar Stereographic) + + Parameters + ---------- + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # coordinate reference system information + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_user_input({'proj':'stere', + 'lat_0':-90, 'lat_ts':-71, 'lon_0':0, 'x_0':0., 'y_0':0., + 'ellps':'WGS84', 'datum':'WGS84', 'units':'km'}) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # wrapper function for models in EPSG 3413 (Sea Ice Polar Stereographic North) + def _EPSG3413(self, EPSG: int | str = 4326): + """ + Transform for models in EPSG:3413 (Sea Ice Polar Stereographic North) + + Parameters + ---------- + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # coordinate reference system information + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_user_input({'proj':'stere', + 'lat_0':90, 'lat_ts':70, 'lon_0':-45, 'x_0':0., 'y_0':0., + 'ellps':'WGS84', 'datum':'WGS84', 'units':'km'}) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # wrapper function for CATS2008 tide models + def _CATS2008(self, EPSG: int | str = 4326): + """ + Transform for Circum-Antarctic Tidal Simulation models + + Parameters + ---------- + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # coordinate reference system information + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_user_input({'proj':'stere', + 'lat_0':-90, 'lat_ts':-71, 'lon_0':-70, 'x_0':0., 'y_0':0., + 'ellps':'WGS84', 'datum':'WGS84', 'units':'km'}) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # wrapper function for models in EPSG 3976 (NSIDC Sea Ice Stereographic South) + def _EPSG3976(self, EPSG: int | str = 4326): + """ + Transform for models in EPSG:3976 (Sea Ice Polar Stereographic South) + + Parameters + ---------- + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # coordinate reference system information + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_user_input({'proj':'stere', + 'lat_0':-90, 'lat_ts':-70, 'lon_0':0, 'x_0':0., 'y_0':0., + 'ellps':'WGS84', 'datum':'WGS84', 'units':'km'}) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # function for models in (idealized) PSNorth projection + def _PSNorth(self, EPSG: int | str = 4326): + """ + Transform for idealized Arctic Polar Stereographic models + + Parameters + ---------- + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # projections for converting to and from input EPSG + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # wrapper function to pass lat/lon values or convert if EPSG + def _EPSG4326(self, EPSG: int | str = 4326): + """ + Transform for models in EPSG:4326 (WGS84 Latitude/Longitude) + + Parameters + ---------- + EPSG: int, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + crs1 = self.from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + # wrapper function for using custom projections + def _custom(self, PROJ: int | str, EPSG: int | str = 4326): + """ + Transform for models in a custom projection + + Parameters + ---------- + PROJ: int or str + Spatial reference system code for coordinate transformations + EPSG: int or str, default 4326 (WGS84 Latitude/Longitude) + input (``'F'``) or output (``'B'``) coordinate system + """ + # coordinate reference system information + crs1 = self.from_input(EPSG) + crs2 = self.from_input(PROJ) + self.transformer = pyproj.Transformer.from_crs(crs1, crs2, + always_xy=True) + return self + + @property + def direction(self): + """ + ``pyproj`` direction of the coordinate transform + """ + # convert from input coordinates to model coordinates + if (self._direction.upper() == 'F'): + return pyproj.enums.TransformDirection.FORWARD + # convert from model coordinates to coordinates + elif (self._direction.upper() == 'B'): + return pyproj.enums.TransformDirection.INVERSE + + def __str__(self): + """String representation of the ``crs`` object + """ + properties = ['pyTMD.crs'] + properties.append(f" name: {self.name}") + properties.append(f" direction: {self.direction.name}") + return '\n'.join(properties) diff --git a/pyTMD/io/FES.py b/pyTMD/io/FES.py index 9e0cf851..1dba8922 100644 --- a/pyTMD/io/FES.py +++ b/pyTMD/io/FES.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" FES.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (01/2024) Reads files for a tidal model and makes initial calculations to run tide program Includes functions to extract tidal harmonic constants from the @@ -56,6 +56,7 @@ interpolate.py: interpolation routines for spatial data UPDATE HISTORY: + Updated 01/2024: attempt to extract constituent IDs from filenames Updated 06/2023: extract ocean tide model variables for FES2012 Updated 04/2023: added global HAMTIDE11 model using pathlib to define and expand tide model paths @@ -353,6 +354,11 @@ def read_constants( model_file = pathlib.Path(model_file).expanduser() if not model_file.exists(): raise FileNotFoundError(str(model_file)) + # try to parse the constituent from the model file + try: + cons = pyTMD.io.model.parse_file(model_file, raise_error=True) + except ValueError as exc: + cons = str(i) # read constituent from elevation file if kwargs['version'] in ('FES1999','FES2004'): # FES ascii constituent files @@ -367,7 +373,7 @@ def read_constants( lon = extend_array(lon, dlon) hc = extend_matrix(hc) # append extended constituent - constituents.append(str(i), hc) + constituents.append(cons, hc) # set model coordinates setattr(constituents, 'longitude', lon) setattr(constituents, 'latitude', lat) diff --git a/pyTMD/io/OTIS.py b/pyTMD/io/OTIS.py index 59dd2d28..b7f5e122 100644 --- a/pyTMD/io/OTIS.py +++ b/pyTMD/io/OTIS.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" OTIS.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Reads files for a tidal model and makes initial calculations to run tide program Includes functions to extract tidal harmonic constants from OTIS tide models for @@ -55,10 +55,11 @@ https://unidata.github.io/netcdf4-python/netCDF4/index.html PROGRAM DEPENDENCIES: - convert_crs.py: converts lat/lon points to and from projected coordinates + crs.py: Coordinate Reference System (CRS) routines interpolate.py: interpolation routines for spatial data UPDATE HISTORY: + Updated 12/2023: use new crs class for coordinate reprojection Updated 10/2023: fix transport variable entry for TMD3 models Updated 09/2023: prevent overwriting ATLAS compact x and y coordinates Updated 08/2023: changed ESR netCDF4 format to TMD3 format @@ -120,9 +121,9 @@ import warnings import numpy as np import scipy.interpolate +import pyTMD.crs import pyTMD.interpolate import pyTMD.io.constituents -from pyTMD.convert_crs import convert_crs # attempt imports try: @@ -239,7 +240,7 @@ def extract_constants( ilon = np.atleast_1d(np.copy(ilon)) ilat = np.atleast_1d(np.copy(ilat)) # run wrapper function to convert coordinate systems of input lat/lon - x,y = convert_crs(ilon, ilat, EPSG, 'F') + x,y = pyTMD.crs().convert(ilon, ilat, EPSG, 'F') # grid step size of tide model dx = xi[1] - xi[0] dy = yi[1] - yi[0] @@ -704,7 +705,7 @@ def interpolate_constants( ilon = np.atleast_1d(np.copy(ilon)) ilat = np.atleast_1d(np.copy(ilat)) # run wrapper function to convert coordinate systems of input lat/lon - x,y = convert_crs(ilon, ilat, EPSG, 'F') + x,y = pyTMD.crs().convert(ilon, ilat, EPSG, 'F') # adjust longitudinal convention of input latitude and longitude # to fit tide model convention if (np.min(x) < np.min(xi)) & (EPSG == '4326'): diff --git a/pyTMD/io/constituents.py b/pyTMD/io/constituents.py index 5691158c..a47a659f 100644 --- a/pyTMD/io/constituents.py +++ b/pyTMD/io/constituents.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" constituents.py -Written by Tyler Sutterley (08/2023) +Written by Tyler Sutterley (01/2024) Basic tide model constituent class PYTHON DEPENDENCIES: @@ -10,6 +10,7 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 01/2024: added properties for Doodson and Cartwright numbers Updated 08/2023: added default for printing constituent class Updated 07/2023: output constituent from get and pop as copy Updated 03/2023: add basic variable typing to function inputs @@ -19,6 +20,7 @@ import copy import numpy as np +import pyTMD.arguments class constituents: """ @@ -152,6 +154,48 @@ def phase(self, field: str): ph.data[ph.mask] = ph.fill_value return ph + @property + def doodson_number(self): + """constituent Doodson number + """ + doodson_numbers = [] + # for each constituent ID + for f in self.fields: + try: + # try to get the Doodson number + n = pyTMD.arguments.doodson_number(f) + except (AssertionError, ValueError) as exc: + n = None + # add Doodson number to the combined list + doodson_numbers.append(n) + # return the list of Doodson numbers + return doodson_numbers + + @property + def cartwright_number(self): + """constituent Cartwright numbers + """ + cartwright_numbers = [] + # for each constituent ID + for f in self.fields: + try: + # try to get the Cartwright numbers + n = pyTMD.arguments.doodson_number(f, formalism='Cartwright') + except (AssertionError, ValueError) as exc: + n = None + # add Cartwright numbers to the combined list + cartwright_numbers.append(n) + # return the list of Cartwright numbers + return cartwright_numbers + + def __str__(self): + """String representation of the ``constituents`` object + """ + properties = ['pyTMD.constituents'] + fields = ', '.join(self.fields) + properties.append(f" constituents: {fields}") + return '\n'.join(properties) + def __len__(self): """Number of constituents """ @@ -174,8 +218,3 @@ def __next__(self): constituent = getattr(self, field) self.__index__ += 1 return (field, constituent) - - def __str__(self): - """Print the list of constituents - """ - return ', '.join(self.fields) diff --git a/pyTMD/io/model.py b/pyTMD/io/model.py index ec148d37..4b2554a0 100644 --- a/pyTMD/io/model.py +++ b/pyTMD/io/model.py @@ -1582,3 +1582,10 @@ def to_bool(self, val: str) -> bool: return False else: raise ValueError(f'Invalid boolean string {val}') + + def __str__(self): + """String representation of the ``io.model`` object + """ + properties = ['pyTMD.io.model'] + properties.append(f" name: {self.name}") + return '\n'.join(properties) diff --git a/pyTMD/load_nodal_corrections.py b/pyTMD/load_nodal_corrections.py index 8941316a..44eb2e2a 100755 --- a/pyTMD/load_nodal_corrections.py +++ b/pyTMD/load_nodal_corrections.py @@ -70,7 +70,7 @@ def load_nodal_corrections(*args, **kwargs): Parameters ---------- MJD: np.ndarray - modified julian day of input date + modified Julian day of input date constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 08c3cecf..14d8b9cc 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" predict.py -Written by Tyler Sutterley (09/2023) +Written by Tyler Sutterley (01/2024) Prediction routines for ocean, load, equilibrium and solid earth tides REFERENCES: @@ -20,6 +20,8 @@ spatial.py: utilities for working with geospatial data UPDATE HISTORY: + Updated 01/2024: moved minor arguments calculation into new function + Updated 12/2023: phase_angles function renamed to doodson_arguments Updated 09/2023: moved constituent parameters function within this module Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 04/2023: using renamed astro mean_longitudes function @@ -43,8 +45,8 @@ from __future__ import annotations import numpy as np +import pyTMD.arguments import pyTMD.astro -from pyTMD.arguments import arguments from pyTMD.constants import constants # PURPOSE: Predict tides at single times @@ -85,22 +87,25 @@ def map(t: float | np.ndarray, .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 """ # number of points and number of constituents - npts,nc = np.shape(hc) + npts, nc = np.shape(hc) # load the nodal corrections # convert time to Modified Julian Days (MJD) - pu,pf,G = arguments(t + 48622.0, constituents, - deltat=deltat, corrections=corrections) + pu, pf, G = pyTMD.arguments.arguments(t + 48622.0, + constituents, + deltat=deltat, + corrections=corrections + ) # allocate for output tidal elevation ht = np.ma.zeros((npts)) - ht.mask = np.zeros((npts),dtype=bool) + ht.mask = np.zeros((npts), dtype=bool) # for each constituent for k,c in enumerate(constituents): - if corrections in ('OTIS','ATLAS','TMD3','netcdf'): + if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = _constituent_parameters(c) # add component for constituent to output tidal elevation th = omega*t*86400.0 + ph + pu[0,k] - elif corrections in ('GOT','FES'): + elif corrections in ('GOT', 'FES'): th = G[0,k]*np.pi/180.0 + pu[0,k] # sum over all tides ht.data[:] += pf[0,k]*hc.real[:,k]*np.cos(th) - \ @@ -150,19 +155,22 @@ def drift(t: float | np.ndarray, nt = len(t) # load the nodal corrections # convert time to Modified Julian Days (MJD) - pu,pf,G = arguments(t + 48622.0, constituents, - deltat=deltat, corrections=corrections) + pu, pf, G = pyTMD.arguments.arguments(t + 48622.0, + constituents, + deltat=deltat, + corrections=corrections + ) # allocate for output time series ht = np.ma.zeros((nt)) - ht.mask = np.zeros((nt),dtype=bool) + ht.mask = np.zeros((nt), dtype=bool) # for each constituent for k,c in enumerate(constituents): - if corrections in ('OTIS','ATLAS','TMD3','netcdf'): + if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = _constituent_parameters(c) # add component for constituent to output tidal elevation th = omega*t*86400.0 + ph + pu[:,k] - elif corrections in ('GOT','FES'): + elif corrections in ('GOT', 'FES'): th = G[:,k]*np.pi/180.0 + pu[:,k] # sum over all tides ht.data[:] += pf[:,k]*hc.real[:,k]*np.cos(th) - \ @@ -212,19 +220,22 @@ def time_series(t: float | np.ndarray, nt = len(t) # load the nodal corrections # convert time to Modified Julian Days (MJD) - pu,pf,G = arguments(t + 48622.0, constituents, - deltat=deltat, corrections=corrections) + pu, pf, G = pyTMD.arguments.arguments(t + 48622.0, + constituents, + deltat=deltat, + corrections=corrections + ) # allocate for output time series ht = np.ma.zeros((nt)) - ht.mask = np.zeros((nt),dtype=bool) + ht.mask = np.zeros((nt), dtype=bool) # for each constituent for k,c in enumerate(constituents): - if corrections in ('OTIS','ATLAS','TMD3','netcdf'): + if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = _constituent_parameters(c) # add component for constituent to output tidal time series th = omega*t*86400.0 + ph + pu[:,k] - elif corrections in ('GOT','FES'): + elif corrections in ('GOT', 'FES'): th = G[:,k]*np.pi/180.0 + pu[:,k] # sum over all tides at location ht.data[:] += pf[:,k]*hc.real[0,k]*np.cos(th) - \ @@ -241,8 +252,8 @@ def infer_minor( **kwargs ): """ - Calculate the tidal corrections for minor constituents inferred using - major constituents [1]_ [2]_ [3]_ [4]_ + Infer the tidal values for minor constituents using their + relation with major constituents [1]_ [2]_ [3]_ [4]_ Parameters ---------- @@ -260,7 +271,7 @@ def infer_minor( Returns ------- dh: np.ndarray - Height from minor constituents + tidal time series for minor constituents References ---------- @@ -292,15 +303,13 @@ def infer_minor( n = nt if ((npts == 1) & (nt > 1)) else npts # allocate for output elevation correction dh = np.ma.zeros((n)) - # convert time from days relative to Jan 1, 1992 to Modified Julian Days - MJD = 48622.0 + t # major constituents used for inferring minor tides cindex = ['q1', 'o1', 'p1', 'k1', 'n2', 'm2', 's2', 'k2', '2n2'] # re-order major tides to correspond to order of cindex z = np.ma.zeros((n,9),dtype=np.complex64) nz = 0 for i,c in enumerate(cindex): - j = [j for j,val in enumerate(constituents) if (val == c)] + j = [j for j,val in enumerate(constituents) if (val.lower() == c)] if j: j1, = j z[:,i] = zmajor[:,j1] @@ -310,13 +319,14 @@ def infer_minor( raise Exception('Not enough constituents for inference') # list of minor constituents - minor = ['2q1','sigma1','rho1','m12','m11','chi1','pi1','phi1','theta1', - 'j1','oo1','2n2','mu2','nu2','lambda2','l2','l2','t2','eps2','eta2'] + minor = ['2q1', 'sigma1', 'rho1', 'm1', 'm1', 'chi1', 'pi1', + 'phi1', 'theta1', 'j1', 'oo1', '2n2', 'mu2', 'nu2', 'lambda2', + 'l2', 'l2', 't2', 'eps2', 'eta2'] # only add minor constituents that are not on the list of major values minor_indices = [i for i,m in enumerate(minor) if m not in constituents] # relationship between major and minor constituent amplitude and phase - zmin = np.zeros((n,20),dtype=np.complex64) + zmin = np.zeros((n, 20), dtype=np.complex64) zmin[:,0] = 0.263*z[:,0] - 0.0252*z[:,1]# 2Q1 zmin[:,1] = 0.297*z[:,0] - 0.0264*z[:,1]# sigma1 zmin[:,2] = 0.164*z[:,0] + 0.0048*z[:,1]# rho1 @@ -351,129 +361,19 @@ def infer_minor( zmin[:,18] = 0.53285*z[:,8] - 0.03304*z[:,4]# eps2 zmin[:,19] = -0.0034925*z[:,5] + 0.0831707*z[:,7]# eta2 - hour = (t % 1)*24.0 - t1 = 15.0*hour - t2 = 30.0*hour - # set function for astronomical longitudes - ASTRO5 = True if kwargs['corrections'] in ('GOT','FES') else False - # convert from Modified Julian Dates into Ephemeris Time - S,H,P,omega,pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'], - ASTRO5=ASTRO5) - - # determine equilibrium tidal arguments - arg = np.zeros((n,20)) - arg[:,0] = t1 - 4.0*S + H + 2.0*P - 90.0# 2Q1 - arg[:,1] = t1 - 4.0*S + 3.0*H - 90.0# sigma1 - arg[:,2] = t1 - 3.0*S + 3.0*H - P - 90.0# rho1 - arg[:,3] = t1 - S + H - P + 90.0# M12 - arg[:,4] = t1 - S + H + P + 90.0# M11 - arg[:,5] = t1 - S + 3.0*H - P + 90.0# chi1 - arg[:,6] = t1 - 2.0*H + pp - 90.0# pi1 - arg[:,7] = t1 + 3.0*H + 90.0# phi1 - arg[:,8] = t1 + S - H + P + 90.0# theta1 - arg[:,9] = t1 + S + H - P + 90.0# J1 - arg[:,10] = t1 + 2.0*S + H + 90.0# OO1 - arg[:,11] = t2 - 4.0*S + 2.0*H + 2.0*P# 2N2 - arg[:,12] = t2 - 4.0*S + 4.0*H# mu2 - arg[:,13] = t2 - 3.0*S + 4.0*H - P# nu2 - arg[:,14] = t2 - S + P + 180.0# lambda2 - arg[:,15] = t2 - S + 2.0*H - P + 180.0# L2 - arg[:,16] = t2 - S + 2.0*H + P# L2 - arg[:,17] = t2 - H + pp# t2 - arg[:,18] = t2 - 5.0*S + 4.0*H + P # eps2 - arg[:,19] = t2 + S + 2.0*H - pp # eta2 - - # determine nodal corrections f and u - sinn = np.sin(omega*dtr) - cosn = np.cos(omega*dtr) - sin2n = np.sin(2.0*omega*dtr) - cos2n = np.cos(2.0*omega*dtr) - - f = np.ones((n,20)) - f[:,0] = np.sqrt((1.0 + 0.189*cosn - 0.0058*cos2n)**2 + - (0.189*sinn - 0.0058*sin2n)**2)# 2Q1 - f[:,1] = f[:,0]# sigma1 - f[:,2] = f[:,0]# rho1 - f[:,3] = np.sqrt((1.0 + 0.185*cosn)**2 + (0.185*sinn)**2)# M12 - f[:,4] = np.sqrt((1.0 + 0.201*cosn)**2 + (0.201*sinn)**2)# M11 - f[:,5] = np.sqrt((1.0 + 0.221*cosn)**2 + (0.221*sinn)**2)# chi1 - f[:,9] = np.sqrt((1.0 + 0.198*cosn)**2 + (0.198*sinn)**2)# J1 - f[:,10] = np.sqrt((1.0 + 0.640*cosn + 0.134*cos2n)**2 + - (0.640*sinn + 0.134*sin2n)**2)# OO1 - f[:,11] = np.sqrt((1.0 - 0.0373*cosn)**2 + (0.0373*sinn)**2)# 2N2 - f[:,12] = f[:,11]# mu2 - f[:,13] = f[:,11]# nu2 - f[:,15] = f[:,11]# L2 - f[:,16] = np.sqrt((1.0 + 0.441*cosn)**2 + (0.441*sinn)**2)# L2 - - u = np.zeros((n,20)) - u[:,0] = np.arctan2(0.189*sinn - 0.0058*sin2n, - 1.0 + 0.189*cosn - 0.0058*sin2n)/dtr# 2Q1 - u[:,1] = u[:,0]# sigma1 - u[:,2] = u[:,0]# rho1 - u[:,3] = np.arctan2( 0.185*sinn, 1.0 + 0.185*cosn)/dtr# M12 - u[:,4] = np.arctan2(-0.201*sinn, 1.0 + 0.201*cosn)/dtr# M11 - u[:,5] = np.arctan2(-0.221*sinn, 1.0 + 0.221*cosn)/dtr# chi1 - u[:,9] = np.arctan2(-0.198*sinn, 1.0 + 0.198*cosn)/dtr# J1 - u[:,10] = np.arctan2(-0.640*sinn - 0.134*sin2n, - 1.0 + 0.640*cosn + 0.134*cos2n)/dtr# OO1 - u[:,11] = np.arctan2(-0.0373*sinn, 1.0 - 0.0373*cosn)/dtr# 2N2 - u[:,12] = u[:,11]# mu2 - u[:,13] = u[:,11]# nu2 - u[:,15] = u[:,11]# L2 - u[:,16] = np.arctan2(-0.441*sinn, 1.0 + 0.441*cosn)/dtr# L2 - - if kwargs['corrections'] in ('FES',): - # additional astronomical terms for FES models - II = np.arccos(0.913694997 - 0.035692561*np.cos(omega*dtr)) - at1 = np.arctan(1.01883*np.tan(omega*dtr/2.0)) - at2 = np.arctan(0.64412*np.tan(omega*dtr/2.0)) - xi = -at1 - at2 + omega*dtr - xi[xi > np.pi] -= 2.0*np.pi - nu = at1 - at2 - I2 = np.tan(II/2.0) - Ra1 = np.sqrt(1.0 - 12.0*(I2**2)*np.cos(2.0*(P - xi)) + 36.0*(I2**4)) - P2 = np.sin(2.0*(P - xi)) - Q2 = 1.0/(6.0*(I2**2)) - np.cos(2.0*(P - xi)) - R = np.arctan(P2/Q2) - - f[:,0] = np.sin(II)*(np.cos(II/2.0)**2)/0.38 # 2Q1 - f[:,1] = f[:,0] # sigma1 - f[:,2] = f[:,0] # rho1 - f[:,3] = f[:,0] # M12 - f[:,4] = np.sin(2.0*II)/0.7214 # M11 - f[:,5] = f[:,4] # chi1 - f[:,9] = f[:,5] # J1 - f[:,10] = np.sin(II)*np.power(np.sin(II/2.0),2.0)/0.01640 # OO1 - f[:,11] = np.power(np.cos(II/2.0),4.0)/0.9154 # 2N2 - f[:,12] = f[:,11] # mu2 - f[:,13] = f[:,11] # nu2 - f[:,14] = f[:,11] # lambda2 - f[:,15] = f[:,11]*Ra1 # L2 - f[:,18] = f[:,11] # eps2 - f[:,19] = np.power(np.sin(II),2.0)/0.1565 # eta2 - - u[:,0] = (2.0*xi - nu)/dtr # 2Q1 - u[:,1] = u[:,0] # sigma1 - u[:,2] = u[:,0] # rho1 - u[:,3] = u[:,0] # M12 - u[:,4] = -nu/dtr # M11 - u[:,5] = u[:,4] # chi1 - u[:,9] = u[:,4] # J1 - u[:,10] = (-2.0*xi - nu)/dtr # OO1 - u[:,11] = (2.0*xi - 2.0*nu)/dtr # 2N2 - u[:,12] = u[:,11] # mu2 - u[:,13] = u[:,11] # nu2 - u[:,14] = (2.0*xi - 2.0*nu)/dtr # lambda2 - u[:,15] = (2.0*xi - 2.0*nu - R)/dtr# L2 - u[:,18] = u[:,12] # eps2 - u[:,19] = -2.0*nu/dtr # eta2 + # load the nodal corrections for minor constituents + # convert time to Modified Julian Days (MJD) + pu, pf, G = pyTMD.arguments.minor_arguments(t + 48622.0, + deltat=kwargs['deltat'], + corrections=kwargs['corrections'] + ) # sum over the minor tidal constituents of interest for k in minor_indices: - th = (arg[:,k] + u[:,k])*dtr - dh += zmin.real[:,k]*f[:,k]*np.cos(th)-zmin.imag[:,k]*f[:,k]*np.sin(th) - # return the inferred elevation + th = (G[:,k] + pu[:,k])*dtr + dh += zmin.real[:,k]*pf[:,k]*np.cos(th) - \ + zmin.imag[:,k]*pf[:,k]*np.sin(th) + # return the inferred values return dh # PURPOSE: estimate long-period equilibrium tides @@ -714,51 +614,56 @@ def _constituent_parameters(c): .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 """ # constituents array that are included in tidal program - cindex = ['m2','s2','k1','o1','n2','p1','k2','q1','2n2','mu2','nu2','l2', - 't2','j1','m1','oo1','rho1','mf','mm','ssa','m4','ms4','mn4','m6','m8', - 'mk3','s6','2sm2','2mk3'] + cindex = ['m2', 's2', 'k1', 'o1', 'n2', 'p1', 'k2', 'q1', '2n2', 'mu2', + 'nu2', 'l2', 't2', 'j1', 'm1', 'oo1', 'rho1', 'mf', 'mm', 'ssa', + 'm4', 'ms4', 'mn4', 'm6', 'm8', 'mk3', 's6', '2sm2', '2mk3'] # species type (spherical harmonic dependence of quadrupole potential) - species_all = np.array([2,2,1,1,2,1,2,1,2,2,2,2,2,1,1,1,1,0,0,0,0,0,0, - 0,0,0,0,0,0]) - # loading love number + _species = np.array([2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) + # Load Love numbers # alpha = correction factor for first order load tides - alpha_all = np.array([0.693,0.693,0.736,0.695,0.693,0.706,0.693,0.695,0.693, - 0.693,0.693,0.693,0.693,0.695,0.695,0.695,0.695,0.693,0.693,0.693, - 0.693,0.693,0.693,0.693,0.693,0.693,0.693,0.693,0.693]) + _alpha = np.array([0.693, 0.693, 0.736, 0.695, 0.693, 0.706, 0.693, + 0.695, 0.693, 0.693, 0.693, 0.693, 0.693, 0.695, 0.695, 0.695, 0.695, + 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, + 0.693, 0.693]) # omega: angular frequency of constituent, in radians - omega_all = np.array([1.405189e-04,1.454441e-04,7.292117e-05,6.759774e-05, - 1.378797e-04,7.252295e-05,1.458423e-04,6.495854e-05,1.352405e-04, - 1.355937e-04,1.382329e-04,1.431581e-04,1.452450e-04,7.556036e-05, - 7.028195e-05,7.824458e-05,6.531174e-05,0.053234e-04,0.026392e-04, - 0.003982e-04,2.810377e-04,2.859630e-04,2.783984e-04,4.215566e-04, - 5.620755e-04,2.134402e-04,4.363323e-04,1.503693e-04,2.081166e-04]) + _omega = np.array([1.405189e-04, 1.454441e-04, 7.292117e-05, 6.759774e-05, + 1.378797e-04, 7.252295e-05, 1.458423e-04, 6.495854e-05, 1.352405e-04, + 1.355937e-04, 1.382329e-04, 1.431581e-04, 1.452450e-04, 7.556036e-05, + 7.028195e-05, 7.824458e-05, 6.531174e-05, 0.053234e-04, 0.026392e-04, + 0.003982e-04, 2.810377e-04, 2.859630e-04, 2.783984e-04, 4.215566e-04, + 5.620755e-04, 2.134402e-04, 4.363323e-04, 1.503693e-04, 2.081166e-04]) # Astronomical arguments (relative to t0 = 1 Jan 0:00 1992) # phases for each constituent are referred to the time when the phase of # the forcing for that constituent is zero on the Greenwich meridian - phase_all = np.array([1.731557546,0.000000000,0.173003674,1.558553872, - 6.050721243,6.110181633,3.487600001,5.877717569,4.086699633, - 3.463115091,5.427136701,0.553986502,0.052841931,2.137025284, - 2.436575100,1.929046130,5.254133027,1.756042456,1.964021610, - 3.487600001,3.463115091,1.731557546,1.499093481,5.194672637, - 6.926230184,1.904561220,0.000000000,4.551627762,3.809122439]) + _phase = np.array([1.731557546, 0.000000000, 0.173003674, 1.558553872, + 6.050721243, 6.110181633, 3.487600001, 5.877717569, 4.086699633, + 3.463115091, 5.427136701, 0.553986502, 0.052841931, 2.137025284, + 2.436575100, 1.929046130, 5.254133027, 1.756042456, 1.964021610, + 3.487600001, 3.463115091, 1.731557546, 1.499093481, 5.194672637, + 6.926230184, 1.904561220, 0.000000000, 4.551627762, 3.809122439]) # amplitudes of equilibrium tide in meters - # amplitude_all = np.array([0.242334,0.112743,0.141565,0.100661,0.046397, - amplitude_all = np.array([0.2441,0.112743,0.141565,0.100661,0.046397, - 0.046848,0.030684,0.019273,0.006141,0.007408,0.008811,0.006931,0.006608, - 0.007915,0.007915,0.004338,0.003661,0.042041,0.022191,0.019567,0.,0.,0., - 0.,0.,0.,0.,0.,0.]) + # _amplitude = np.array([0.242334,0.112743,0.141565,0.100661,0.046397, + _amplitude = np.array([0.2441, 0.112743, 0.141565, 0.100661, 0.046397, + 0.046848, 0.030684, 0.019273, 0.006141, 0.007408, 0.008811, 0.006931, + 0.006608, 0.007915, 0.007915, 0.004338, 0.003661, 0.042041, 0.022191, + 0.019567, 0., 0., 0., 0., 0., 0., 0., 0., 0.]) # map between input constituent and cindex j = [j for j,val in enumerate(cindex) if (val == c.lower())] # set the values for the constituent if j: - amplitude, = amplitude_all[j] - phase, = phase_all[j] - omega, = omega_all[j] - alpha, = alpha_all[j] - species, = species_all[j] + amplitude, = _amplitude[j] + phase, = _phase[j] + omega, = _omega[j] + alpha, = _alpha[j] + species, = _species[j] else: - amplitude = 0.0; phase = 0.0; omega = 0.0; alpha = 0.0; species = 0 + amplitude = 0.0 + phase = 0.0 + omega = 0.0 + alpha = 0.0 + species = 0 # return the values for the constituent return (amplitude, phase, omega, alpha, species) @@ -1014,8 +919,8 @@ def _frequency_dependence_diurnal( [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [3.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]) - # get phase angles - S, H, P, TAU, ZNS, PS = pyTMD.astro.phase_angles(MJD) + # get phase angles (Doodson arguments) + TAU, S, H, P, ZNS, PS = pyTMD.astro.doodson_arguments(MJD) # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius @@ -1070,8 +975,8 @@ def _frequency_dependence_long_period( [2.0, 0.0, 0.0, 0.0, 0.0, -0.13, -0.11, -0.15, -0.07], [2.0, 0.0, 0.0, 1.0, 0.0, -0.05, -0.05, -0.06, -0.03] ]) - # get phase angles - S, H, P, TAU, ZNS, PS = pyTMD.astro.phase_angles(MJD) + # get phase angles (Doodson arguments) + TAU, S, H, P, ZNS, PS = pyTMD.astro.doodson_arguments(MJD) # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius diff --git a/pyTMD/solve.py b/pyTMD/solve.py new file mode 100644 index 00000000..437748c9 --- /dev/null +++ b/pyTMD/solve.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python +u""" +solve.py +Written by Tyler Sutterley (12/2023) +Routines for estimating the harmonic constants for ocean tides + +REFERENCES: + G. D. Egbert and S. Erofeeva, "Efficient Inverse Modeling of Barotropic + Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002). + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python + https://numpy.org + https://numpy.org/doc/stable/user/numpy-for-matlab-users.html + scipy: Scientific Tools for Python + https://scipy.org + +PROGRAM DEPENDENCIES: + arguments.py: loads nodal corrections for tidal constituents + astro.py: computes the basic astronomical mean longitudes + predict.py: predict tidal values using harmonic constants + +UPDATE HISTORY: + Written 12/2023 +""" + +from __future__ import annotations + +import numpy as np +import scipy.linalg +import pyTMD.predict +from pyTMD.arguments import arguments + +def constants(t: float | np.ndarray, + ht: np.ndarray, + constituents: list | np.ndarray, + deltat: float | np.ndarray = 0.0, + corrections: str = 'OTIS', + solver: str = 'lstsq' + ): + """ + Estimate the harmonic constants for an elevation time series [1]_ + + Parameters + ---------- + t: float or np.ndarray + days relative to 1992-01-01T00:00:00 + ht: np.ndarray + elevation time series (meters) + constituents: list or np.ndarray + tidal constituent IDs + deltat: float or np.ndarray, default 0.0 + time correction for converting to Ephemeris Time (days) + corrections: str, default 'OTIS' + use nodal corrections from OTIS/ATLAS or GOT/FES models + solver: str, default 'lstsq' + least squares solver to use + + - ``'lstsq'``: least squares solution + - ``'gelsy'``: complete orthogonal factorization + - ``'gelss'``: singular value decomposition (SVD) + - ``'gelsd'``: SVD with divide and conquer method + + Returns + ------- + amp: np.ndarray + amplitude of each harmonic constant (meters) + phase: np.ndarray + phase of each harmonic constant (degrees) + + References + ---------- + .. [1] G. D. Egbert and S. Y. Erofeeva, "Efficient Inverse Modeling of + Barotropic Ocean Tides," *Journal of Atmospheric and Oceanic + Technology*, 19(2), 183--204, (2002). + `doi: 10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2`__ + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + # check that there are enough values for a time series fit + nt = len(np.atleast_1d(t)) + nc = len(constituents) + if (nt <= 2*nc): + raise ValueError('Not enough values for fit') + # check that the number of time values matches the number of height values + if (nt != len(np.atleast_1d(ht))): + raise ValueError('Dimension mismatch between input variables') + + # load the nodal corrections + # convert time to Modified Julian Days (MJD) + pu, pf, G = arguments(t + 48622.0, constituents, + deltat=deltat, corrections=corrections) + + # create design matrix + M = [] + # add constant term for mean + M.append(np.ones_like(t)) + # add constituent terms + for k,c in enumerate(constituents): + if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): + amp, ph, omega, alpha, species = \ + pyTMD.predict._constituent_parameters(c) + th = omega*t*86400.0 + ph + pu[:,k] + elif corrections in ('GOT', 'FES'): + th = G[:,k]*np.pi/180.0 + pu[:,k] + # add constituent to design matrix + M.append(pf[:,k]*np.cos(th)) + M.append(-pf[:,k]*np.sin(th)) + # take the transpose of the design matrix + M = np.transpose(M) + + # use a least-squares fit to solve for parameters + if (solver == 'lstsq'): + p, res, rnk, s = np.linalg.lstsq(M, ht, rcond=-1) + elif solver in ('gelsd', 'gelsy', 'gelss'): + p, res, rnk, s = scipy.linalg.lstsq(M, ht, + lapack_driver=solver) + + # calculate amplitude and phase for each constituent + amp = np.zeros((nc)) + ph = np.zeros((nc)) + # skip over the first indice in the fit (constant term) + for k,c in enumerate(constituents): + amp[k] = np.abs(1j*p[2*k+2] + p[2*k+1]) + ph[k] = np.arctan2(-p[2*k+2], p[2*k+1]) + # convert phase to degrees + phase = ph*180.0/np.pi + phase[phase < 0] += 360.0 + + # return the amplitude and phase + return (amp, phase) diff --git a/pyTMD/time.py b/pyTMD/time.py index 3f2c4b45..d4caf301 100644 --- a/pyTMD/time.py +++ b/pyTMD/time.py @@ -551,7 +551,7 @@ def convert_julian(JD: np.ndarray, **kwargs): else: single_value = False - # verify julian day + # verify Julian day JDO = np.floor(JD + 0.5) C = np.zeros_like(JD) # calculate C for dates before and after the switch to Gregorian @@ -903,6 +903,12 @@ def ndim(self): """ return np.ndim(self.MJD) + def __str__(self): + """String representation of the ``timescale`` object + """ + properties = ['pyTMD.time.timescale'] + return '\n'.join(properties) + def __len__(self): """Number of time values """ diff --git a/pyTMD/tools.py b/pyTMD/tools.py index 75c94abf..a90dc632 100644 --- a/pyTMD/tools.py +++ b/pyTMD/tools.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" tools.py -Written by Tyler Sutterley (08/2023) +Written by Tyler Sutterley (12/2023) Jupyter notebook, user interface and plotting tools PYTHON DEPENDENCIES: @@ -17,6 +17,7 @@ https://github.com/matplotlib/matplotlib UPDATE HISTORY: + Updated 12/2023: pass through VBox and HBox Updated 08/2023: place matplotlib within try/except statements Updated 05/2023: don't set a default directory for tide models Updated 04/2023: using pathlib to define and expand paths @@ -71,6 +72,9 @@ def __init__(self, **kwargs): kwargs.setdefault('style', {}) # set style self.style = copy.copy(kwargs['style']) + # pass through some ipywidgets objects + self.HBox = ipywidgets.HBox + self.VBox = ipywidgets.VBox # set the directory with tide models self.directory = ipywidgets.Text( diff --git a/scripts/compute_LPET_elevations.py b/scripts/compute_LPET_elevations.py index 1e7ecf87..b2d8642c 100644 --- a/scripts/compute_LPET_elevations.py +++ b/scripts/compute_LPET_elevations.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_LPET_elevations.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates long-period equilibrium tidal elevations for an input file INPUTS: @@ -61,12 +61,14 @@ https://pypi.org/project/pyproj/ PROGRAM DEPENDENCIES: + crs.py: Coordinate Reference System (CRS) routines time.py: utilities for calculating time operations spatial.py: utilities for reading and writing spatial data utilities.py: download and management utilities for syncing files predict.py: calculates long-period equilibrium ocean tides UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 05/2023: use timescale class for time conversion operations Updated 04/2023: check if datetime before converting to seconds @@ -110,21 +112,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/compute_LPT_displacements.py b/scripts/compute_LPT_displacements.py index 2ca3532a..dac94c3e 100644 --- a/scripts/compute_LPT_displacements.py +++ b/scripts/compute_LPT_displacements.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_LPT_displacements.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates radial load pole tide displacements for an input file following IERS Convention (2010) guidelines https://iers-conventions.obspm.fr/chapter7.php @@ -72,12 +72,14 @@ PROGRAM DEPENDENCIES: constants.py: gravitational and ellipsoidal parameters + crs.py: Coordinate Reference System (CRS) routines eop.py: utilities for calculating Earth Orientation Parameters (EOP) spatial: utilities for reading, writing and operating on spatial data time.py: utilities for calculating time operations utilities.py: download and management utilities for syncing files UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 05/2023: use timescale class for time conversion operations use defaults from eop module for pole tide and EOP files @@ -130,21 +132,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/compute_OPT_displacements.py b/scripts/compute_OPT_displacements.py index 36c1b99b..3550e710 100644 --- a/scripts/compute_OPT_displacements.py +++ b/scripts/compute_OPT_displacements.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_OPT_displacements.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates radial ocean pole load tide displacements for an input file following IERS Convention (2010) guidelines https://iers-conventions.obspm.fr/chapter7.php @@ -75,6 +75,7 @@ https://pypi.org/project/pyproj/ PROGRAM DEPENDENCIES: + crs.py: Coordinate Reference System (CRS) routines time.py: utilities for calculating time operations spatial.py: utilities for reading and writing spatial data utilities.py: download and management utilities for syncing files @@ -89,6 +90,7 @@ doi: 10.1007/s00190-015-0848-7 UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 05/2023: use timescale class for time conversion operations use defaults from eop module for pole tide and EOP files @@ -147,21 +149,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/compute_SET_displacements.py b/scripts/compute_SET_displacements.py index df355e20..c70fdfdf 100644 --- a/scripts/compute_SET_displacements.py +++ b/scripts/compute_SET_displacements.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_SET_displacements.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates radial solid earth tide displacements for an input file following IERS Convention (2010) guidelines https://iers-conventions.obspm.fr/chapter7.php @@ -72,6 +72,7 @@ https://pypi.org/project/pyproj/ PROGRAM DEPENDENCIES: + crs.py: Coordinate Reference System (CRS) routines time.py: utilities for calculating time operations spatial.py: utilities for reading and writing spatial data utilities.py: download and management utilities for syncing files @@ -96,6 +97,7 @@ doi: 10.1111/j.1365-246X.1981.tb02690.x UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 05/2023: use timescale class for time conversion operations add option for using higher resolution ephemerides from JPL @@ -125,21 +127,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/compute_tidal_currents.py b/scripts/compute_tidal_currents.py index 86a8cd65..254dbe4a 100755 --- a/scripts/compute_tidal_currents.py +++ b/scripts/compute_tidal_currents.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_currents.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates zonal and meridional tidal currents for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -85,7 +85,7 @@ utilities.py: download and management utilities for syncing files arguments.py: load the nodal corrections for tidal constituents astro.py: computes the basic astronomical mean longitudes - convert_crs.py: convert points to and from Coordinates Reference Systems + crs.py: Coordinate Reference System (CRS) routines io/model.py: retrieves tide model parameters for named tide models io/OTIS.py: extract tidal harmonic constants from OTIS tide models io/ATLAS.py: extract tidal harmonic constants from netcdf models @@ -94,6 +94,7 @@ predict.py: predict tidal values using harmonic constants UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 05/2023: use timescale class for time conversion operations @@ -160,21 +161,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/compute_tidal_elevations.py b/scripts/compute_tidal_elevations.py index 86e14b94..4e68bc29 100755 --- a/scripts/compute_tidal_elevations.py +++ b/scripts/compute_tidal_elevations.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_elevations.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (12/2023) Calculates tidal elevations for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -88,7 +88,7 @@ utilities.py: download and management utilities for syncing files arguments.py: load the nodal corrections for tidal constituents astro.py: computes the basic astronomical mean longitudes - convert_crs.py: convert points to and from Coordinates Reference Systems + crs.py: Coordinate Reference System (CRS) routines io/model.py: retrieves tide model parameters for named tide models io/OTIS.py: extract tidal harmonic constants from OTIS tide models io/ATLAS.py: extract tidal harmonic constants from netcdf models @@ -97,6 +97,7 @@ predict.py: predict tidal values using harmonic constants UPDATE HISTORY: + Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 05/2023: use timescale class for time conversion operations @@ -165,21 +166,14 @@ def get_projection(attributes, PROJECTION): # coordinate reference system string from file try: - crs = pyproj.CRS.from_string(attributes['projection']) + crs = pyTMD.crs().from_input(attributes['projection']) except (ValueError,KeyError,pyproj.exceptions.CRSError): pass else: return crs - # EPSG projection code + # coordinate reference system from input argument try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) + crs = pyTMD.crs().from_input(PROJECTION) except (ValueError,pyproj.exceptions.CRSError): pass else: diff --git a/scripts/reduce_OTIS_files.py b/scripts/reduce_OTIS_files.py index c9f4700f..4a98a178 100644 --- a/scripts/reduce_OTIS_files.py +++ b/scripts/reduce_OTIS_files.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" reduce_OTIS_files.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (12/2023) Read OTIS-format tidal files and reduce to a regional subset COMMAND LINE OPTIONS: @@ -25,9 +25,10 @@ io/model.py: retrieves tide model parameters for named tide models io/OTIS.py: extract tidal harmonic constants from OTIS tide models utilities.py: download and management utilities for syncing files - convert_crs.py: converts lat/lon points to and from projected coordinates + crs.py: Coordinate Reference System (CRS) routines UPDATE HISTORY: + Updated 12/2023: use new crs class for coordinate reprojection Updated 04/2023: using pathlib to define and expand paths Updated 03/2023: new function name for coordinate reference systems Updated 12/2022: refactored OTIS model input and output @@ -54,9 +55,9 @@ import pathlib import argparse import numpy as np +import pyTMD.crs import pyTMD.io import pyTMD.utilities -from pyTMD.convert_crs import convert_crs # attempt imports try: @@ -64,25 +65,6 @@ except (AttributeError, ImportError, ModuleNotFoundError) as exc: logging.critical("pyproj not available") -# PURPOSE: try to get the projection information for the input file -def get_projection(PROJECTION): - # EPSG projection code - try: - crs = pyproj.CRS.from_epsg(int(PROJECTION)) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # coordinate reference system string - try: - crs = pyproj.CRS.from_string(PROJECTION) - except (ValueError,pyproj.exceptions.CRSError): - pass - else: - return crs - # no projection can be made - raise pyproj.exceptions.CRSError - # PURPOSE: reads OTIS-format tidal files and reduces to a regional subset def make_regional_OTIS_files(tide_dir, TIDE_MODEL, BOUNDS=4*[None], PROJECTION='4326', MODE=0o775): @@ -106,7 +88,7 @@ def make_regional_OTIS_files(tide_dir, TIDE_MODEL, BOUNDS=4*[None], xi,yi,hz,mz,iob,dt = pyTMD.io.OTIS.read_otis_grid(model.grid_file) # converting bounds x,y from projection to latitude/longitude - crs1 = get_projection(PROJECTION) + crs1 = pyTMD.crs().from_input(PROJECTION) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) xbox = np.array([BOUNDS[0],BOUNDS[1],BOUNDS[1],BOUNDS[0],BOUNDS[0]]) @@ -114,7 +96,7 @@ def make_regional_OTIS_files(tide_dir, TIDE_MODEL, BOUNDS=4*[None], lon,lat = transformer.transform(xbox,ybox) # convert bounds from latitude/longitude to model coordinates - x,y = convert_crs(lon,lat,model.projection,'F') + x,y = pyTMD.crs().convert(lon,lat,model.projection,'F') # find indices to reduce to xmin,xmax,ymin,ymax gridx,gridy = np.meshgrid(xi,yi) diff --git a/test/test_arguments.py b/test/test_arguments.py new file mode 100644 index 00000000..6dc638ff --- /dev/null +++ b/test/test_arguments.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python +u""" +test_arguments.py (01/2024) +Verify arguments table matches prior arguments array +""" +import pytest +import numpy as np +import pyTMD.astro +from pyTMD.arguments import ( + doodson_number, + _arguments_table, + _minor_table, + _to_doodson_number, + _from_doodson_number +) + +@pytest.mark.parametrize("corrections", ['OTIS', 'GOT']) +def test_arguments(corrections): + """ + Tests the calculation of nodal corrections for tidal constituents + + Parameters + ---------- + corrections: str, default 'OTIS' + use nodal corrections from OTIS/ATLAS or GOT models + """ + # use a random set of modified Julian days + MJD = np.random.randint(58000, 61000, size=10) + # set function for astronomical longitudes + ASTRO5 = True if corrections in ('GOT', 'FES') else False + # convert from Modified Julian Dates into Ephemeris Time + s, h, p, omega, pp = pyTMD.astro.mean_longitudes(MJD, + ASTRO5=ASTRO5) + # number of temporal values + nt = len(np.atleast_1d(MJD)) + # initial time conversions + hour = 24.0*np.mod(MJD, 1) + # convert from hours into degrees + t1 = 15.0*hour + t2 = 30.0*hour + # convert from hours solar time into mean lunar time in degrees + tau = 15.0*hour - s + h + # variable for multiples of 90 degrees (Ray technical note 2017) + k = 90.0 + np.zeros((nt)) + + # determine equilibrium arguments + arg = np.zeros((nt, 60)) + arg[:,0] = h - pp # Sa + arg[:,1] = 2.0*h # Ssa + arg[:,2] = s - p # Mm + arg[:,3] = 2.0*s - 2.0*h # MSf + arg[:,4] = 2.0*s # Mf + arg[:,5] = 3.0*s - p # Mt + arg[:,6] = t1 - 5.0*s + 3.0*h + p - 90.0 # alpha1 + arg[:,7] = t1 - 4.0*s + h + 2.0*p - 90.0 # 2Q1 + arg[:,8] = t1 - 4.0*s + 3.0*h - 90.0 # sigma1 + arg[:,9] = t1 - 3.0*s + h + p - 90.0 # q1 + arg[:,10] = t1 - 3.0*s + 3.0*h - p - 90.0 # rho1 + arg[:,11] = t1 - 2.0*s + h - 90.0 # o1 + arg[:,12] = t1 - 2.0*s + 3.0*h + 90.0 # tau1 + arg[:,13] = t1 - s + h + 90.0 # M1 + arg[:,14] = t1 - s + 3.0*h - p + 90.0 # chi1 + arg[:,15] = t1 - 2.0*h + pp - 90.0 # pi1 + arg[:,16] = t1 - h - 90.0 # p1 + if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): + arg[:,17] = t1 + 90.0 # s1 + elif corrections in ('GOT', 'FES'): + arg[:,17] = t1 + 180.0 # s1 (Doodson's phase) + arg[:,18] = t1 + h + 90.0 # k1 + arg[:,19] = t1 + 2.0*h - pp + 90.0 # psi1 + arg[:,20] = t1 + 3.0*h + 90.0 # phi1 + arg[:,21] = t1 + s - h + p + 90.0 # theta1 + arg[:,22] = t1 + s + h - p + 90.0 # J1 + arg[:,23] = t1 + 2.0*s + h + 90.0 # OO1 + arg[:,24] = t2 - 4.0*s + 2.0*h + 2.0*p # 2N2 + arg[:,25] = t2 - 4.0*s + 4.0*h # mu2 + arg[:,26] = t2 - 3.0*s + 2.0*h + p # n2 + arg[:,27] = t2 - 3.0*s + 4.0*h - p # nu2 + arg[:,28] = t2 - 2.0*s + h + pp # M2a + arg[:,29] = t2 - 2.0*s + 2.0*h # M2 + arg[:,30] = t2 - 2.0*s + 3.0*h - pp # M2b + arg[:,31] = t2 - s + p + 180.0 # lambda2 + arg[:,32] = t2 - s + 2.0*h - p + 180.0 # L2 + arg[:,33] = t2 - h + pp # T2 + arg[:,34] = t2 # S2 + arg[:,35] = t2 + h - pp + 180.0 # R2 + arg[:,36] = t2 + 2.0*h # K2 + arg[:,37] = t2 + s + 2.0*h - pp # eta2 + arg[:,38] = t2 - 5.0*s + 4.0*h + p # MNS2 + arg[:,39] = t2 + 2.0*s - 2.0*h # 2SM2 + arg[:,40] = 1.5*arg[:,29] # M3 + arg[:,41] = arg[:,18] + arg[:,29] # MK3 + arg[:,42] = 3.0*t1 # S3 + arg[:,43] = arg[:,26] + arg[:,29] # MN4 + arg[:,44] = 2.0*arg[:,29] # M4 + arg[:,45] = arg[:,29] + arg[:,34] # MS4 + arg[:,46] = arg[:,29] + arg[:,36] # MK4 + arg[:,47] = 4.0*t1 # S4 + arg[:,48] = 5.0*t1 # S5 + arg[:,49] = 3.0*arg[:,29] # M6 + arg[:,50] = 3.0*t2 # S6 + arg[:,51] = 7.0*t1 # S7 + arg[:,52] = 4.0*t2 # S8 + # shallow water constituents + arg[:,53] = 4.0*arg[:,29] # m8 + arg[:,54] = arg[:,29] + arg[:,36] - arg[:,34] # mks2 + arg[:,55] = 4.0*s - 2.0*h # msqm + arg[:,56] = 3.0*s - p # mtm + arg[:,57] = 2.0*arg[:,26] # n4 + arg[:,58] = t2 - 5.0*s + 4.0*h + p # eps2 + # mean sea level + arg[:,59] = 0.0 # Z0 + + # determine equilibrium arguments using table + fargs = np.c_[tau, s, h, p, omega, pp, k] + test = np.dot(fargs, _arguments_table(corrections=corrections)) + + # validate arguments between methods + assert np.all(arg == test) + +def test_minor(): + """ + Tests the calculation of nodal corrections for minor tidal constituents + """ + # use a random set of modified Julian days + MJD = np.random.randint(58000, 61000, size=10) + # convert from Modified Julian Dates into Ephemeris Time + s, h, p, n, pp = pyTMD.astro.mean_longitudes(MJD) + # number of temporal values + nt = len(np.atleast_1d(MJD)) + # initial time conversions + hour = 24.0*np.mod(MJD, 1) + # convert from hours into degrees + t1 = 15.0*hour + t2 = 30.0*hour + # convert from hours solar time into mean lunar time in degrees + tau = 15.0*hour - s + h + # variable for multiples of 90 degrees (Ray technical note 2017) + k = 90.0 + np.zeros((nt)) + + # determine equilibrium tidal arguments + arg = np.zeros((nt, 20)) + arg[:,0] = t1 - 4.0*s + h + 2.0*p - 90.0# 2Q1 + arg[:,1] = t1 - 4.0*s + 3.0*h - 90.0# sigma1 + arg[:,2] = t1 - 3.0*s + 3.0*h - p - 90.0# rho1 + arg[:,3] = t1 - s + h - p + 90.0# M12 + arg[:,4] = t1 - s + h + p + 90.0# M11 + arg[:,5] = t1 - s + 3.0*h - p + 90.0# chi1 + arg[:,6] = t1 - 2.0*h + pp - 90.0# pi1 + arg[:,7] = t1 + 3.0*h + 90.0# phi1 + arg[:,8] = t1 + s - h + p + 90.0# theta1 + arg[:,9] = t1 + s + h - p + 90.0# J1 + arg[:,10] = t1 + 2.0*s + h + 90.0# OO1 + arg[:,11] = t2 - 4.0*s + 2.0*h + 2.0*p# 2N2 + arg[:,12] = t2 - 4.0*s + 4.0*h# mu2 + arg[:,13] = t2 - 3.0*s + 4.0*h - p# nu2 + arg[:,14] = t2 - s + p + 180.0# lambda2 + arg[:,15] = t2 - s + 2.0*h - p + 180.0# L2 + arg[:,16] = t2 - s + 2.0*h + p# L2 + arg[:,17] = t2 - h + pp# t2 + arg[:,18] = t2 - 5.0*s + 4.0*h + p # eps2 + arg[:,19] = t2 + s + 2.0*h - pp # eta2 + + # determine equilibrium arguments + fargs = np.c_[tau, s, h, p, n, pp, k] + test = np.dot(fargs, _minor_table()) + + # validate arguments between methods + assert np.all(arg == test) + +def test_doodson(): + """ + Tests the calculation of Doodson numbers + """ + # expected values + exp = {} + # semi-diurnal species + exp['m2'] = 255.555 + exp['s2'] = 273.555 + exp['n2'] = 245.655 + exp['nu2'] = 247.455 + exp['mu2'] = 237.555 + exp['2n2'] = 235.755 + exp['lambda2'] = 263.655 + exp['l2'] = 265.455 + exp['k2'] = 275.555 + # diurnal species + exp['m1'] = 155.555 + exp['s1'] = 164.555 + exp['o1'] = 145.555 + exp['oo1'] = 185.555 + exp['k1'] = 165.555 + exp['q1'] = 135.655 + exp['2q1'] = 125.755 + exp['p1'] = 163.555 + # long-period species + exp['mm'] = 65.455 + exp['ssa'] = 57.555 + exp['msf'] = 73.555 + exp['mf'] = 75.555 + exp['msqm'] = 93.555 + exp['mtm'] = 85.455 + # short-period species + exp['m3'] = 355.555 + exp['m4'] = 455.555 + exp['m6'] = 655.555 + exp['m8'] = 855.555 + # get observed values for constituents + obs = doodson_number(exp.keys()) + cartwright = doodson_number(exp.keys(), formalism='Cartwright') + # check values + for key,val in exp.items(): + assert val == obs[key] + # check values when entered as string + test = doodson_number(key) + assert val == test + # check conversion to and from Doodson numbers + doodson = _to_doodson_number(cartwright[key]) + # check values when entered as Cartwright + assert val == doodson + # check values when entered as Doodson + coefficients = _from_doodson_number(val) + assert np.all(cartwright[key] == coefficients) diff --git a/test/test_atlas_read.py b/test/test_atlas_read.py index c83b34e0..291ded8e 100644 --- a/test/test_atlas_read.py +++ b/test/test_atlas_read.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -test_atlas_read.py (04/2023) +test_atlas_read.py (01/2024) Tests that ATLAS compact and netCDF4 data can be downloaded from AWS S3 bucket Tests the read program to verify that constituents are being extracted @@ -16,6 +16,7 @@ https://boto3.amazonaws.com/v1/documentation/api/latest/index.html UPDATE HISTORY: + Updated 01/2024: test doodson and cartwright numbers of each constituent Updated 04/2023: using pathlib to define and expand paths Updated 12/2022: add check for read and interpolate constants Updated 11/2022: use f-strings for formatting verbose or ascii output @@ -38,6 +39,7 @@ import pyTMD.io import pyTMD.time import pyTMD.io.model +import pyTMD.arguments import pyTMD.utilities # current file path @@ -197,6 +199,13 @@ def test_compare_TPXO9_v2(METHOD): # calculate complex form of constituent oscillation hc2 = amp2*np.exp(-1j*ph2*np.pi/180.0) + # expected Doodson numbers for constituents + exp = {} + exp['m2'] = 255.555 + exp['s2'] = 273.555 + exp['o1'] = 145.555 + exp['k1'] = 165.555 + # will verify differences between model outputs are within tolerance eps = np.finfo(np.float16).eps # calculate differences between methods @@ -206,6 +215,11 @@ def test_compare_TPXO9_v2(METHOD): # calculate difference in amplitude and phase difference = hc1[:,i] - hc2[:,i] assert np.all(np.abs(difference) <= eps) + # verify doodson numbers + assert (constituents.doodson_number[i] == exp[cons]) + # verify cartwright numbers + assert np.all(constituents.cartwright_number[i] == + pyTMD.arguments._from_doodson_number(exp[cons])) # parameterize interpolation method @pytest.mark.parametrize("METHOD", ['bilinear']) diff --git a/test/test_coordinates.py b/test/test_coordinates.py index 8efe3dd7..f453f99c 100644 --- a/test/test_coordinates.py +++ b/test/test_coordinates.py @@ -1,11 +1,15 @@ #!/usr/bin/env python u""" -test_coordinates.py (08/2020) +test_coordinates.py (12/2023) Verify forward and backwards coordinate conversions + +UPDATE HISTORY: + Updated 12/2023: use new crs class for coordinate reprojection + Written 08/2020 """ import pytest import numpy as np -import pyTMD.convert_crs +import pyTMD.crs # parameterize projections @pytest.mark.parametrize("PROJ", ['3031','CATS2008','3976','PSNorth','4326']) @@ -16,8 +20,8 @@ def test_coordinates(PROJ): i1 = np.arange(-180,180+1,1) i2 = np.linspace(startlat[PROJ],endlat[PROJ],len(i1)) # convert latitude and longitude to and from projection - o1, o2 = pyTMD.convert_crs(i1,i2,PROJ,'F') - lon, lat = pyTMD.convert_crs(o1,o2,PROJ,'B') + o1, o2 = pyTMD.crs().convert(i1,i2,PROJ,'F') + lon, lat = pyTMD.crs().convert(o1,o2,PROJ,'B') # calculate great circle distance between inputs and outputs cdist = np.arccos(np.sin(i2*np.pi/180.0)*np.sin(lat*np.pi/180.0) + np.cos(i2*np.pi/180.0)*np.cos(lat*np.pi/180.0)* diff --git a/test/test_download_and_read.py b/test/test_download_and_read.py index 0855e332..041fa6f3 100644 --- a/test/test_download_and_read.py +++ b/test/test_download_and_read.py @@ -56,6 +56,7 @@ import pyTMD.utilities import pyTMD.check_points import pyTMD.ellipse +import pyTMD.solve from oct2py import octave # current file path @@ -546,6 +547,47 @@ def test_tidal_ellipse(self): if not np.all(difference.mask): assert np.all(np.abs(difference) < eps) + # PURPOSE: Tests solving for harmonic constants + def test_solve(self): + # get model parameters + model = pyTMD.io.model(filepath).elevation('CATS2008') + + # calculate a forecast every minute + minutes = np.arange(366*1440) + # convert time to days relative to Jan 1, 1992 (48622 MJD) + year, month, day = 2000, 1, 1 + tide_time = pyTMD.time.convert_calendar_dates(year, month, day, + minute=minutes) + + # read tidal constants and interpolate to coordinates + constituents = pyTMD.io.OTIS.read_constants( + model.grid_file, model.model_file, + model.projection, type=model.type, + grid=model.format) + c = constituents.fields + DELTAT = np.zeros_like(tide_time) + + # interpolate constants to a coordinate + LAT,LON = (-76.0, -40.0) + amp,ph,D = pyTMD.io.OTIS.interpolate_constants( + np.atleast_1d(LON), np.atleast_1d(LAT), + constituents, model.projection, type=model.type, + method='spline', extrapolate=True) + + # calculate complex form of constituent oscillation + hc = amp*np.exp(-1j*ph*np.pi/180.0) + # predict tidal elevations at times + TIDE = pyTMD.predict.time_series(tide_time, hc, c, + deltat=DELTAT, corrections=model.format) + # solve for amplitude and phase + famp, fph = pyTMD.solve.constants(tide_time, TIDE.data, c) + # calculate complex form of constituent oscillation + fhc = famp*np.exp(-1j*fph*np.pi/180.0) + # verify differences are within tolerance + eps = 5e-3 + for k,cons in enumerate(c): + assert np.isclose(hc[0][k], fhc[k], rtol=eps, atol=eps) + # parameterize type: heights versus currents # parameterize interpolation method parameters = [] diff --git a/test/test_fes_predict.py b/test/test_fes_predict.py index 3a4674a1..0954dc92 100644 --- a/test/test_fes_predict.py +++ b/test/test_fes_predict.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -test_fes_predict.py (04/2023) +test_fes_predict.py (01/2024) Tests that FES2014 data can be downloaded from AWS S3 bucket Tests the read program to verify that constituents are being extracted Tests that interpolated results are comparable to FES2014 program @@ -17,6 +17,7 @@ https://boto3.amazonaws.com/v1/documentation/api/latest/index.html UPDATE HISTORY: + Updated 01/2024: test doodson and cartwright numbers of each constituent Updated 04/2023: using pathlib to define and expand paths Updated 12/2022: add check for read and interpolate constants Updated 09/2021: update check tide points to add compression flags @@ -38,6 +39,7 @@ import pyTMD.io.model import pyTMD.utilities import pyTMD.predict +import pyTMD.arguments import pyTMD.check_points # current file path @@ -192,6 +194,24 @@ def test_compare_FES2014(METHOD): # calculate complex form of constituent oscillation hc2 = amp2*np.exp(-1j*ph2*np.pi/180.0) + # expected Doodson numbers for constituents + exp = {} + exp['2n2'] = 235.755 + exp['k1'] = 165.555 + exp['k2'] = 275.555 + exp['m2'] = 255.555 + exp['m4'] = 455.555 + exp['mf'] = 75.555 + exp['mm'] = 65.455 + exp['msqm'] = 93.555 + exp['mtm'] = 85.455 + exp['n2'] = 245.655 + exp['o1'] = 145.555 + exp['p1'] = 163.555 + exp['q1'] = 135.655 + exp['s1'] = 164.555 + exp['s2'] = 273.555 + # will verify differences between model outputs are within tolerance eps = np.finfo(np.float16).eps # calculate differences between methods @@ -199,6 +219,11 @@ def test_compare_FES2014(METHOD): # calculate difference in amplitude and phase difference = hc1[:,i] - hc2[:,i] assert np.all(np.abs(difference) <= eps) + # verify doodson numbers + assert (constituents.doodson_number[i] == exp[cons]) + # verify cartwright numbers + assert np.all(constituents.cartwright_number[i] == + pyTMD.arguments._from_doodson_number(exp[cons])) # validate iteration within constituents class for field, hc in constituents: diff --git a/test/test_solid_earth.py b/test/test_solid_earth.py index 7e27500d..5c5d77f9 100644 --- a/test/test_solid_earth.py +++ b/test/test_solid_earth.py @@ -1,5 +1,5 @@ """ -test_solid_earth.py (04/2023) +test_solid_earth.py (12/2023) Tests the steps for calculating the solid earth tides PYTHON DEPENDENCIES: @@ -8,6 +8,7 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 12/2023: phase_angles function renamed to doodson_arguments Updated 04/2023: added test for using JPL ephemerides for positions Written 04/2023 """ @@ -129,7 +130,7 @@ def test_phase_angles(): s, h, p, N, PP = pyTMD.astro.mean_longitudes(MJD, ASTRO5=True) PR = dtr*pyTMD.astro.polynomial_sum(np.array([0.0, 1.396971278, 3.08889e-4, 2.1e-8, 7.0e-9]), T) - S, H, P, TAU, ZNS, PS = pyTMD.astro.phase_angles(MJD) + TAU, S, H, P, ZNS, PS = pyTMD.astro.doodson_arguments(MJD) assert np.isclose(dtr*s + PR, S) assert np.isclose(dtr*h, H) assert np.isclose(dtr*p, P) diff --git a/test/test_spatial.py b/test/test_spatial.py index ecd105e4..6b596db5 100644 --- a/test/test_spatial.py +++ b/test/test_spatial.py @@ -400,7 +400,7 @@ def test_convert_geodetic(): assert np.isclose(longitude, ln1).all() assert np.isclose(latitude, lt1).all() assert np.isclose(height, h1).all() - # validate outputs for Bowring iterative method + # validate outputs for Bowring iterative method assert np.isclose(longitude, ln2).all() assert np.isclose(latitude, lt2).all() assert np.isclose(height, h2).all()