Starting functions dispatch to sub-modules

vfrecovery/command_line_interface/utils.py

@@ -0,0 +1,37 @@
+import sys
+import logging
+
+log = logging.getLogger("vfrecovery.cli")
+
+
+PREF = "\033["
+RESET = f"{PREF}0m"
+class COLORS:
+    black = "30m"
+    red = "31m"
+    green = "32m"
+    yellow = "33m"
+    blue = "34m"
+    magenta = "35m"
+    cyan = "36m"
+    white = "37m"
+
+
+def puts(text, color=None, bold=False, file=sys.stdout):
+    """Alternative to print, uses no color by default but accepts any color from the COLORS class.
+
+    Parameters
+    ----------
+    text
+    color=None
+    bold=False
+    file=sys.stdout
+    """
+    if color is None:
+        txt = f'{PREF}{1 if bold else 0}m' + text + RESET
+        print(txt, file=file)
+    else:
+        txt = f'{PREF}{1 if bold else 0};{color}' + text + RESET
+        print(txt, file=file)
+    log.info(text)
+

vfrecovery/download_functions/__init__.py

@@ -0,0 +1,2 @@
+from armor3d import Armor3d
+from glorys import Glorys

vfrecovery/download_functions/armor3d.py

@@ -0,0 +1,116 @@
+import pandas as pd
+import copernicusmarine
+
+
+class Armor3d:
+    """Global Ocean 1/4° Multi Observation Product ARMOR3D
+
+    Product description:
+    https://data.marine.copernicus.eu/product/MULTIOBS_GLO_PHY_TSUV_3D_MYNRT_015_012
+
+    If start_date + n_days <= 2022-12-28:
+        Delivers the multi-year reprocessed (REP) weekly data
+
+    otherwise:
+        Delivers the near-real-time (NRT) weekly data
+
+    Examples
+    --------
+    >>> Armor3d([-25, -13, 6.5, 13], pd.to_datetime('20091130', utc=True)).to_xarray()
+    >>> Armor3d([-25, -13, 6.5, 13], pd.to_datetime('20231121', utc=True), n_days=10).to_xarray()
+
+    """
+
+    def __init__(self, box, start_date, n_days=1, max_depth=2500):
+        """
+        Parameters
+        ----------
+        box: list(float)
+            Define domain to load: [lon_min, lon_max, lat_min, lat_max]
+        start_date: :class:`pandas.Timestamp`
+            Starting date of the time series to load. Since ARMOR3D is weekly, the effective starting
+            date will be the first weekly period including the user-defined ``start_date``
+        n_days: int (default=1)
+            Number of days to load data for.
+        max_depth: float (default=2500)
+            Maximum depth levels to load data for.
+        """
+        self.box = box
+        self.start_date = start_date
+        self.n_days = n_days
+        self.max_depth = max_depth
+
+        dt = pd.Timedelta(n_days, 'D') if n_days > 1 else pd.Timedelta(0, 'D')
+        if start_date + dt <= pd.to_datetime('2022-12-28', utc=True):
+            self._loader = self._get_rep
+            self.dataset_id = "dataset-armor-3d-rep-weekly"
+            self.time_axis = pd.Series(pd.date_range('19930106', '20221228', freq='7D').tz_localize("UTC"))
+        else:
+            self._loader = self._get_nrt
+            self.dataset_id = "dataset-armor-3d-nrt-weekly"
+            self.time_axis = pd.Series(
+                pd.date_range('20190102', pd.to_datetime('now', utc=True).strftime("%Y%m%d"), freq='7D').tz_localize(
+                    "UTC")[0:-1])
+
+        if start_date < self.time_axis.iloc[0]:
+            raise ValueError('Date out of bounds')
+        elif start_date + dt > self.time_axis.iloc[-1]:
+            raise ValueError('Date out of bounds, %s > %s' % (
+                start_date + dt, self.time_axis.iloc[-1]))
+
+    def _get_this(self, dataset_id):
+        start_date = self.time_axis[self.time_axis <= self.start_date].iloc[-1]
+        if self.n_days == 1:
+            end_date = start_date
+        else:
+            end_date = \
+            self.time_axis[self.time_axis <= self.start_date + (self.n_days + 1) * pd.Timedelta(1, 'D')].iloc[-1]
+
+        ds = copernicusmarine.open_dataset(
+            dataset_id=dataset_id,
+            minimum_longitude=self.box[0],
+            maximum_longitude=self.box[1],
+            minimum_latitude=self.box[2],
+            maximum_latitude=self.box[3],
+            maximum_depth=self.max_depth,
+            start_datetime=start_date.strftime("%Y-%m-%dT%H:%M:%S"),
+            end_datetime=end_date.strftime("%Y-%m-%dT%H:%M:%S"),
+            variables=['ugo', 'vgo']
+        )
+        return ds
+
+    def _get_rep(self):
+        """multi-year reprocessed (REP) weekly data
+
+        Returns
+        -------
+        :class:xarray.dataset
+        """
+        return self._get_this(self.dataset_id)
+
+    def _get_nrt(self):
+        """near-real-time (NRT) weekly data
+
+        Returns
+        -------
+        :class:xarray.dataset
+        """
+        return self._get_this(self.dataset_id)
+
+    def to_xarray(self):
+        """Load and return data as a :class:`xarray.dataset`
+
+        Returns
+        -------
+        :class:xarray.dataset
+        """
+        return self._loader()
+
+    def __repr__(self):
+        summary = ["<CopernicusMarineData.Loader><Armor3D>"]
+        summary.append("dataset_id: %s" % self.dataset_id)
+        summary.append("First day: %s" % self.start_date)
+        summary.append("N days: %s" % self.n_days)
+        summary.append("Domain: %s" % self.box)
+        summary.append("Max depth (m): %s" % self.max_depth)
+        return "\n".join(summary)

vfrecovery/download_functions/glorys.py

@@ -0,0 +1,154 @@
+import os
+import numpy as np
+import glob
+import pandas as pd
+import xarray as xr
+import copernicusmarine
+
+
+def get_glorys_forecast_from_datarmor(a_box, a_start_date, n_days=1):
+    """Load Datarmor Global Ocean 1/12° Physics Analysis and Forecast updated Daily
+
+    Fields: daily, from 2020-11-25T12:00 to 'now' + 5 days
+    Src: /home/ref-ocean-model-public/multiparameter/physic/global/cmems/global-analysis-forecast-phy-001-024
+    Info: https://resources.marine.copernicus.eu/product-detail/GLOBAL_ANALYSISFORECAST_PHY_001_024/INFORMATION
+
+    Parameters
+    ----------
+    a_box
+    a_start_date
+    n_days
+    """
+    def get_forecast_files(a_date, n_days=1):
+        file_list = []
+        for n in range(0, n_days):
+            t = a_date + pd.Timedelta(n, 'D')
+            p = os.path.join(src, "%i" % t.year, "%0.3d" % t.day_of_year)
+            # print(p, os.path.exists(p))
+            if os.path.exists(p):
+                file_list.append(sorted(glob.glob(os.path.join(p, "*.nc")))[0])
+        return file_list
+
+    def preprocess(this_ds):
+        idpt = np.argwhere(this_ds['depth'].values > 2000)[0][0]
+        ilon = np.argwhere(this_ds['longitude'].values >= a_box[0])[0][0], \
+               np.argwhere(this_ds['longitude'].values >= a_box[1])[0][0]
+        ilat = np.argwhere(this_ds['latitude'].values >= a_box[2])[0][0], \
+               np.argwhere(this_ds['latitude'].values >= a_box[3])[0][0]
+        this_ds = this_ds.isel({'depth': range(0, idpt),
+                                'longitude': range(ilon[0], ilon[1]),
+                                'latitude': range(ilat[0], ilat[1])})
+        return this_ds
+
+    root = "/home/ref-ocean-model-public" if not os.uname()[0] == 'Darwin' else "/Volumes/MODEL-PUBLIC/"
+    src = os.path.join(root, "multiparameter/physic/global/cmems/global-analysis-forecast-phy-001-024")
+    # puts("\t%s" % src, color=COLORS.green)
+    flist = get_forecast_files(a_start_date, n_days=n_days)
+    if len(flist) == 0:
+        raise ValueError("This float cycle is too old for this velocity field.")
+    glorys = xr.open_mfdataset(flist, preprocess=preprocess, combine='nested', concat_dim='time', parallel=True)
+    #
+    return glorys
+
+
+class Glorys:
+    """Global Ocean 1/12° Physics Re-Analysis or Forecast
+
+    If start_date + n_days <= 2021-01-09:
+        Delivers the multi-year reprocessed (REP) daily data
+        https://resources.marine.copernicus.eu/product-detail/GLOBAL_MULTIYEAR_PHY_001_030
+
+    otherwise:
+        Delivers the near-real-time (NRT) Analysis and Forecast daily data
+        https://resources.marine.copernicus.eu/product-detail/GLOBAL_ANALYSISFORECAST_PHY_001_024
+
+    Examples
+    --------
+    >>> Glorys([-25, -13, 6.5, 13], pd.to_datetime('20091130', utc=True)).to_xarray()
+    >>> Glorys([-25, -13, 6.5, 13], pd.to_datetime('20231121', utc=True), n_days=10).to_xarray()
+
+    """
+
+    def __init__(self, box, start_date, n_days=1, max_depth=2500):
+        """
+        Parameters
+        ----------
+        box: list(float)
+            Define domain to load: [lon_min, lon_max, lat_min, lat_max]
+        start_date: :class:`pandas.Timestamp`
+            Starting date of the time series to load.
+        n_days: int (default=1)
+            Number of days to load data for.
+        max_depth: float (default=2500)
+            Maximum depth levels to load data for.
+        """
+        self.box = box
+        self.start_date = start_date
+        self.n_days = n_days
+        self.max_depth = max_depth
+
+        dt = pd.Timedelta(n_days, 'D') if n_days > 1 else pd.Timedelta(0, 'D')
+        if start_date + dt <= pd.to_datetime('2021-01-09', utc=True):
+            self._loader = self._get_reanalysis
+            self.dataset_id = "cmems_mod_glo_phy_my_0.083_P1D-m"
+        else:
+            self._loader = self._get_forecast
+            self.dataset_id = "cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m"
+
+    def _get_this(self, dataset_id, dates):
+        ds = copernicusmarine.open_dataset(
+            dataset_id=dataset_id,
+            minimum_longitude=self.box[0],
+            maximum_longitude=self.box[1],
+            minimum_latitude=self.box[2],
+            maximum_latitude=self.box[3],
+            maximum_depth=self.max_depth,
+            start_datetime=dates[0].strftime("%Y-%m-%dT%H:%M:%S"),
+            end_datetime=dates[1].strftime("%Y-%m-%dT%H:%M:%S"),
+            variables=['uo', 'vo']
+        )
+        return ds
+
+    def _get_forecast(self):
+        """
+        Returns
+        -------
+        :class:`xarray.dataset`
+        """
+        start_date = self.start_date
+        if self.n_days == 1:
+            end_date = start_date
+        else:
+            end_date = start_date + pd.Timedelta(self.n_days - 1, 'D')
+        return self._get_this(self.dataset_id, [start_date, end_date])
+
+    def _get_reanalysis(self):
+        """
+        Returns
+        -------
+        :class:`xarray.dataset`
+        """
+        start_date = self.start_date
+        if self.n_days == 1:
+            end_date = start_date
+        else:
+            end_date = self.start_date + pd.Timedelta(self.n_days - 1, 'D')
+        return self._get_this(self.dataset_id, [start_date, end_date])
+
+    def to_xarray(self):
+        """ Load and return data as a :class:`xarray.dataset`
+        Returns
+        -------
+        :class:`xarray.dataset`
+        """
+        return self._loader()
+
+    def __repr__(self):
+        summary = ["<CopernicusMarineData.Loader><Glorys>"]
+        summary.append("dataset_id: %s" % self.dataset_id)
+        summary.append("First day: %s" % self.start_date)
+        summary.append("N days: %s" % self.n_days)
+        summary.append("Domain: %s" % self.box)
+        summary.append("Max depth (m): %s" % self.max_depth)
+        return "\n".join(summary)
+

vfrecovery/utilities/formatters.py

@@ -0,0 +1,18 @@
+
+def strfdelta(tdelta, fmt):
+    """
+
+    Parameters
+    ----------
+    tdelta
+    fmt
+
+    Returns
+    -------
+
+    """
+    d = {"days": tdelta.days}
+    d["hours"], rem = divmod(tdelta.seconds, 3600)
+    d["minutes"], d["seconds"] = divmod(rem, 60)
+    return fmt.format(**d)
+

vfrecovery/utilities/geo.py

@@ -0,0 +1,63 @@
+from math import radians, cos, sin, asin, sqrt
+import pyproj
+
+
+def haversine(lon1, lat1, lon2, lat2):
+    """
+    Calculate the great circle distance (in [km]) between two points
+    on the earth (specified in decimal degrees)
+
+    see: https://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points
+
+    Parameters
+    ----------
+    lon1
+    lat1
+    lon2
+    lat2
+
+    Returns
+    -------
+    km
+    """
+    # convert decimal degrees to radians
+    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
+
+    # haversine formula
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
+    c = 2 * asin(sqrt(a))
+    r = 6371  # Radius of earth in kilometers.
+    return c * r
+
+
+def bearing(lon1, lat1, lon2, lat2):
+    """
+
+    Parameters
+    ----------
+    lon1
+    lat1
+    lon2
+    lat2
+
+    Returns
+    -------
+
+    """
+    # from math import cos, sin, atan2, degrees
+    # b = atan2(cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1), sin(lon2 - lon1) * cos(lat2))
+    # b = degrees(b)
+    # return b
+
+    geodesic = pyproj.Geod(ellps='WGS84')
+    fwd_azimuth, back_azimuth, distance = geodesic.inv(lon1, lat1, lon2, lat2)
+    return fwd_azimuth
+
+
+def fixLON(x):
+    """Ensure a 0-360 longitude"""
+    if x < 0:
+        x = 360 + x
+    return x