A new kernel to handle the specific use case of Argo floats recovery

docs/whats-new.rst

@@ -11,8 +11,20 @@ What's New
 Coming up next
 --------------
 
-- Virtual Argo float configuration manager :class:`FloatConfiguration` now uses a well documented `JSON schema <https://raw.githubusercontent.com/euroargodev/VirtualFleet/master/schemas/VF-ArgoFloat-Configuration.json>`_ to load/validate and save
-the set of parameters. (:pr:`29`) by `G. Maze <http://www.github.com/gmaze>`_.
+- A new kernel was developed to handle the specific use case of Argo floats recovery. With this new kernel it is possible
+to let the float free drift at the surface after a given number of simulated cycles. This kernel will automatically be
+selected when the ``free_surface_drift`` parameter is found in the :class:`FloatConfiguration` of the deployment plan. This
+can be used as follows:
+
+.. code-block:: python
+
+    cfg = FloatConfiguration([6903091, 120])  # Load data from the last cycle of the float to recover
+    cfg.add('free_surface_drift', 4)  #
+
+With this configuration, the virtual float will start to free drift at the surface after the 4th simulated cycle.
+
+- Virtual Argo float configuration manager :class:`FloatConfiguration` now uses a well documented `JSON schema <https://raw.githubusercontent.com/euroargodev/VirtualFleet/master/schemas/VF-ArgoFloat-Configuration.json>`_ to load/validate/save the set of parameters. (:pr:`29`) by `G. Maze <http://www.github.com/gmaze>`_.
+
 
 v0.4.0 (2 Feb. 2024)
 --------------------

environment.yml

@@ -7,7 +7,7 @@ dependencies:
   - parcels>=3.0.0
   - dask
   - distributed
-  - dask-kubernetes
+#  - dask-kubernetes
   - bottleneck
   - gcsfs
   - zarr

virtualargofleet/app_parcels.py

@@ -53,6 +53,10 @@ class ArgoParticle(JITParticle):
     life_expectancy = Variable('life_expectancy', dtype=np.int32, initial=200, to_write=False)
     """Float mission parameter life expectancy in cycle"""
 
+    # Kernel specifics parameters:
+    free_surface_drift = Variable('free_surface_drift', dtype=np.int8, initial=9999, to_write=False)
+    """Integer indicating when the virtual float will start to free drift at the surface"""
+
 
 def ArgoFloatKernel(particle, fieldset, time):
     """Default kernel to simulate an Argo float
@@ -247,7 +251,7 @@ class ArgoParticle_exp(ArgoParticle):
     :class:`parcels.particle.JITParticle`
     """
     in_area = Variable('in_area', dtype=np.float32, initial=0., to_write=False)
-    """Boolean indicating if the virtual float in the experiment area (1) or not (0)"""
+    """Boolean indicating if the virtual float is in the experiment area (1) or not (0)"""
 
 
 def ArgoFloatKernel_exp(particle, fieldset, time):
@@ -372,6 +376,20 @@ def ArgoFloatKernel_exp(particle, fieldset, time):
         particle.cycle_age += particle.dt  # update cycle_age
 
 
+def ArgoFloatKernel_recovery(particle, fieldset, time):
+    """
+    To be chained after the standard ArgoFloatKernel
+    """
+    if particle.cycle_number >= particle.free_surface_drift:
+        if fieldset.verbose_events:
+            print("This float will start to free drift at the surface")
+
+        particle_ddepth = 0  # Reset change in depth
+
+        if particle.cycle_phase == 0:
+            particle.cycle_phase = 5
+
+
 def PeriodicBoundaryConditionKernel(particle, fieldset, time):
     """Define periodic Boundary Conditions."""
     if particle.lon < fieldset.halo_west:
@@ -411,9 +429,10 @@ def KeepInWater(particle, fieldset, time):
 #             # Go throught KeepInDomain
 #             pass 
 
+
 def KeepInDomain(particle, fieldset, time):
     # out of geographical area : here we can delete the particle
     if particle.state == StatusCode.ErrorOutOfBounds:
-        if fieldset.verbose_events == 1:            
-            print("Field warning : Float out of the horizontal geographical domain OR interpolation error --> deleted")
-        particle.delete()
+        if fieldset.verbose_events == 1:
+            print("Field warning : Float out of the horizontal geographical domain --> deleted")
+        particle.delete()

virtualargofleet/assets/FloatConfiguration_recovery.json

@@ -0,0 +1,68 @@
+{
+    "created": "2024-03-07T08:44:43.560904+00:00",
+    "version": "2.0",
+    "name": "recovery",
+    "parameters": [
+        {
+            "name": "cycle_duration",
+            "value": 240.0,
+            "description": "Maximum length of float complete cycle",
+            "meta": {
+                "unit": "hours",
+                "dtype": "float",
+                "techkey": "CONFIG_CycleTime_hours"
+            }
+        },
+        {
+            "name": "life_expectancy",
+            "value": 200,
+            "description": "Maximum number of completed cycle",
+            "meta": {
+                "unit": "cycle",
+                "dtype": "int",
+                "techkey": "CONFIG_MaxCycles_NUMBER"
+            }
+        },
+        {
+            "name": "parking_depth",
+            "value": 1000.0,
+            "description": "Drifting depth",
+            "meta": {
+                "unit": "m",
+                "dtype": "float",
+                "techkey": "CONFIG_ParkPressure_dbar"
+            }
+        },
+        {
+            "name": "profile_depth",
+            "value": 2000.0,
+            "description": "Maximum profile depth",
+            "meta": {
+                "unit": "m",
+                "dtype": "float",
+                "techkey": "CONFIG_ProfilePressure_dbar"
+            }
+        },
+        {
+            "name": "vertical_speed",
+            "value": 0.09,
+            "description": "Vertical profiling speed",
+            "meta": {
+                "unit": "m/s",
+                "dtype": "float",
+                "techkey": ""
+            }
+        },
+        {
+            "name": "reco_free_surface_drift",
+            "value": 9999,
+            "description": "Virtual cycle number to start free surface drift, inclusive",
+            "meta": {
+                "unit": "",
+                "dtype": "int",
+                "techkey": ""
+            }
+        }
+    ],
+    "$schema": "https://raw.githubusercontent.com/euroargodev/VirtualFleet/json-schemas-FloatConfiguration/schemas/VF-ArgoFloat-Configuration.json"
+}

virtualargofleet/utilities.py

@@ -252,6 +252,7 @@ class FloatConfiguration:
     --------
     >>> cfg = FloatConfiguration('default')  # Internally defined
     >>> cfg = FloatConfiguration('local-change')  # Internally defined
+    >>> cfg = FloatConfiguration('recovery')  # Internally defined
     >>> cfg = FloatConfiguration('cfg_file.json')  # From any json file
     >>> cfg = FloatConfiguration([6902919, 132])  # From Euro-Argo Fleet API
     >>> cfg.update('parking_depth', 500)  # Update one parameter value
@@ -336,6 +337,9 @@ def load_from_json(name):
         elif name == 'gse-experiment' or name == "gulf-stream":
             name, data = load_from_json(os.path.join(path2data, 'FloatConfiguration_gulf_stream.json'))
 
+        elif name == 'recovery':
+            name, data = load_from_json(os.path.join(path2data, 'FloatConfiguration_recovery.json'))
+
         elif isinstance(name, str) and os.path.splitext(name)[-1] == ".json":
             name, data = load_from_json(name)
 
@@ -873,3 +877,22 @@ def getSystemInfo():
         return info
     except Exception as e:
         logging.exception(e)
+
+
+def save_figurefile(this_fig, a_name, folder: Path = Path('.')):
+    """
+    Parameters
+    ----------
+    this_fig
+    a_name
+
+    Returns
+    -------
+    path
+    """
+    # figname = os.path.join(folder, "%s.png" % a_name)
+    figname = folder.joinpath("%s.png" % a_name)
+    # log.debug("Saving %s ..." % figname)
+    this_fig.savefig(figname)
+    return figname
+

virtualargofleet/velocity_helpers.py

@@ -1,16 +1,38 @@
-"""
-Velocity Field Helper
-"""
-
 from parcels import FieldSet, ParticleSet, Field
 import xarray as xr
+import pandas as pd
+import numpy as np
 import glob
 from abc import ABC
 import logging
+from pathlib import Path
+
 from .app_parcels import ArgoParticle
+from .utilities import save_figurefile
+
+logger = logging.getLogger("virtualfleet.velocity")
+
+class default_logger:
 
+    def __init__(self, txt, log_level):
+        """Log text"""
+        getattr(logger, log_level.lower())(txt)
 
-log = logging.getLogger("virtualfleet.velocity")
+    @staticmethod
+    def info(txt) -> 'default_logger':
+        return default_logger(txt, 'INFO')
+
+    @staticmethod
+    def debug(txt) -> 'default_logger':
+        return default_logger(txt, 'DEBUG')
+
+    @staticmethod
+    def warning(txt) -> 'default_logger':
+        return default_logger(txt, 'WARNING')
+
+    @staticmethod
+    def error(txt) -> 'default_logger':
+        return default_logger(txt, 'ERROR')
 
 
 class VelocityField(ABC):
@@ -62,6 +84,7 @@ def add_mask(self):
             - ``self.var`` with ``U`` and ``V`` keys
         """
         if self.fieldset:
+            self.logger.debug("Create mask for grounding management from velocity dataset")
             if isinstance(self.field, xr.core.dataset.Dataset):
                 ds = self.field[{self.dim['time']: 0}]
                 ds = ds[[self.var['U'], self.var['V']]].squeeze()
@@ -112,12 +135,13 @@ def __init__(self,
                  isglobal: bool = False,
                  **kwargs):
         """Create a custom VelocityField for known products"""
-
         if 'U' not in variables:
             raise ValueError("'variables' dictionary must have a 'U' key")
         if 'V' not in variables:
             raise ValueError("'variables' dictionary must have a 'V' key")
 
+        self.logger = kwargs['logger'] if 'logger' in kwargs else default_logger
+
         if isinstance(src, xr.core.dataset.Dataset):
             # If the source data is a Xarray dataset, we ensure to have all the required variables and coordinates:
             for v in variables:
@@ -145,13 +169,15 @@ def __init__(self,
 
         # Define parcels fieldset
         if not isinstance(src, xr.core.dataset.Dataset):
+            self.logger.debug("Define parcels fieldset from_netcdf")
             self.field = src  # Dictionary with 'U' and 'V' as keys and list of corresponding files as values
             self.fieldset = FieldSet.from_netcdf(
                 src, self.var, self.dim,
                 allow_time_extrapolation=True,
                 time_periodic=False,
                 deferred_load=True)
         else:
+            self.logger.debug("Define parcels fieldset from_xarray_dataset")
             self.field = src  # Xarray dataset
             self.fieldset = FieldSet.from_xarray_dataset(
                 src, self.var, self.dim,
@@ -164,6 +190,42 @@ def __init__(self,
         # Create mask to manage grounding:
         self.add_mask()
 
+    def plot(self, it=0, iz=0, save: bool = False, workdir: Path = Path('.')):
+        """
+
+        Parameters
+        ----------
+        save
+        workdir
+
+        Returns
+        -------
+        fig, ax
+        """
+        import matplotlib.pyplot as plt
+        import argopy.plot as argoplot
+        import cartopy.crs as ccrs
+
+        fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 20), dpi=100,
+                               subplot_kw={'projection': ccrs.PlateCarree()})
+
+        # ax.set_extent(box)
+        argoplot.utils.latlongrid(ax)
+        ax.add_feature(argoplot.utils.land_feature, edgecolor="black")
+
+        self.field.isel(time=it, depth=iz).plot.quiver(x="longitude", y="latitude",
+                                                    u=self.var['U'], v=self.var['V'], ax=ax, color='grey', alpha=0.5,
+                                                    add_guide=False)
+
+        plt.tight_layout()
+        if save:
+            t = pd.to_datetime(self.field.isel(time=it, depth=iz)['time'].data).strftime("%Y%m%d")
+            z = str(np.abs(self.field.isel(time=it, depth=iz)['depth'].data)).split('.')[0]
+            figname = save_figurefile(fig, 'velocity_%s_t%s_z%s' % (self.name, t, z), workdir)
+            return fig, ax, figname
+        else:
+            return fig, ax
+
 
 def VelocityFieldFacade(model: str = 'GLOBAL_ANALYSIS_FORECAST_PHY_001_024', *args: object, **kwargs: object) -> object:
     """Function to return a :class:`VelocityField` instance for known products
@@ -248,14 +310,17 @@ def VelocityField_PSY4QV3R1(**kwargs):
     else:
         src = {'U': kwargs['src'],
                'V': kwargs['src']}
+    kwargs.pop('src', None)
+
     variables = {'U': 'uo',
                  'V': 'vo'}
     dimensions = {'time': 'time',
                   'depth': 'depth',
                   'lat': 'latitude',
                   'lon': 'longitude'}
     isglobal = kwargs['isglobal'] if 'isglobal' in kwargs else False
-    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=isglobal)
+    kwargs.pop('isglobal', None)
+    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=isglobal, **kwargs)
     V.name = 'PSY4QV3R1'
     return V
 
@@ -275,13 +340,14 @@ def VelocityField_MEDSEA_ANALYSISFORECAST_PHY_006_013(**kwargs):
     else:
         src = {'U': kwargs['src'],
                'V': kwargs['src']}
+    kwargs.pop('src', None)
     variables = {'U': 'uo',
                  'V': 'vo'}
     dimensions = {'time': 'time',
                   'depth': 'depth',
                   'lat': 'lat',
                   'lon': 'lon'}
-    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=False)
+    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=False, **kwargs)
     V.name = 'MEDSEA_ANALYSISFORECAST_PHY_006_013'
     return V
 
@@ -301,13 +367,15 @@ def VelocityField_ARMOR3D(**kwargs):
     else:
         src = {'U': kwargs['src'],
                'V': kwargs['src']}
+    kwargs.pop('src', None)
     variables = {'U': 'ugo',
                  'V': 'vgo'}
     dimensions = {'time': 'time',
                   'depth': 'depth',
                   'lat': 'latitude',
                   'lon': 'longitude'}
     isglobal = kwargs['isglobal'] if 'isglobal' in kwargs else False
-    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=isglobal)
+    kwargs.pop('isglobal', None)
+    V = VelocityField_CUSTOM(src=src, variables=variables, dimensions=dimensions, isglobal=isglobal, **kwargs)
     V.name = 'ARMOR3D'
     return V