Merge pull request #21 from axiom-data-science/some_updates

Some updates

docs/whats_new.md

@@ -1,5 +1,18 @@
 # What's New
 
+## v0.8.0 (April 2, 2024)
+
+* `time_step_output` behavior has changed — 1 hour by default
+* `time_step` is now 5 min by default
+* added `Dcrit` parameter for accurately finding where drifters are stranded in tidal flats
+* `vertical_mixing` is True by default now
+* added seafloor_action option
+* fixed some Leeway/3D handling and log messaging
+* export_variables are specific to drift_model as needed
+* do not drop zeta anymore since used in opendrift
+* output_file is now an option
+
+
 ## v0.7.1 (February 21, 2024)
 
 * Small fix to some attributes to be less verbose

particle_tracking_manager/models/opendrift/config.json

@@ -7,14 +7,15 @@
     },
     "export_variables": {
         "default": [
-            "z"
+            "z",
+            "origin_marker"
         ],
         "ptm_level": 2,
         "type": "list",
         "description": "List of variables to export. Options available with `m.all_export_variables` for a given `drift_model`. ['lon', 'lat', 'ID', 'status'] will always be exported."
     },
     "radius": {
-        "default": 0.0,
+        "default": 1000.0,
         "ptm_level": 1,
         "type": "float",
         "min": 0.0,
@@ -58,6 +59,11 @@
         "od_mapping": "general:coastline_action",
         "ptm_level": 1
     },
+    "seafloor_action": {
+        "default": "previous",
+        "od_mapping": "general:seafloor_action",
+        "ptm_level": 1
+    },
     "max_speed": {
         "default": 2,
         "od_mapping": "drift:max_speed"
@@ -198,5 +204,10 @@
         "default": "low",
         "ptm_level": 2,
         "description": "Log verbosity"
+    },
+    "Dcrit": {
+        "default": 0.1,
+        "od_mapping": "general:seafloor_action_dcrit",
+        "ptm_level": 3
     }
 }

particle_tracking_manager/models/opendrift/opendrift.py

@@ -84,6 +84,9 @@ class OpenDriftModel(ParticleTrackingManager):
     coastline_action : str, optional
         Action to perform if a drifter hits the coastline, by default "previous". Options
         are 'stranding', 'previous'.
+    seafloor_action : str, optional
+        Action to perform if a drifter hits the seafloor, by default "deactivate". Options
+        are 'deactivate', 'previous', 'lift_to_seafloor'.
     max_speed : int
         Typical maximum speed of elements, used to estimate reader buffer size.
     wind_drift_factor : float
@@ -173,6 +176,7 @@ def __init__(
         stokes_drift: bool = config_model["stokes_drift"]["default"],
         mixed_layer_depth: float = config_model["mixed_layer_depth"]["default"],
         coastline_action: str = config_model["coastline_action"]["default"],
+        seafloor_action: str = config_model["seafloor_action"]["default"],
         max_speed: int = config_model["max_speed"]["default"],
         wind_drift_factor: float = config_model["wind_drift_factor"]["default"],
         wind_drift_depth: float = config_model["wind_drift_depth"]["default"],
@@ -319,6 +323,15 @@ def __setattr_model__(self, name: str, value) -> None:
             self.config_model[name]["value"] = value
         self._update_config()
 
+        if name == "ocean_model":
+            if value == "NWGOA":
+                self.Dcrit = 0.5
+            elif "CIOFS" in value:
+                self.Dcrit = 0.3
+            else:
+                self.Dcrit = 0.1
+            self.logger.info(f"For ocean_model {value}, setting Dcrit to {self.Dcrit}.")
+
         if name in ["ocean_model", "horizontal_diffusivity"]:
 
             # just set the value and move on if purposely setting a non-None value
@@ -332,12 +345,11 @@ def __setattr_model__(self, name: str, value) -> None:
             # in all other cases that ocean_model is a known model, want to use the
             # grid-dependent value
             elif self.ocean_model in _KNOWN_MODELS:
-                print(name, value)
-                self.logger.info(
-                    "Setting horizontal_diffusivity parameter to one tuned to reader model"
-                )
 
                 hdiff = self.calc_known_horizontal_diffusivity()
+                self.logger.info(
+                    f"Setting horizontal_diffusivity parameter to one tuned to reader model of value {hdiff}."
+                )
                 # when editing the __dict__ directly have to also update config_model
                 self.__dict__["horizontal_diffusivity"] = hdiff
                 self.config_model["horizontal_diffusivity"]["value"] = hdiff
@@ -381,9 +393,21 @@ def __setattr_model__(self, name: str, value) -> None:
 
         # Leeway doesn't have this option available
         if name == "do3D" and not value and self.drift_model != "Leeway":
+            self.logger.info("do3D is False so disabling vertical motion.")
             self.o.disable_vertical_motion()
-        elif name == "do3D" and value:
+        elif name == "do3D" and not value and self.drift_model == "Leeway":
+            self.logger.info(
+                "do3D is False but drift_model is Leeway so doing nothing."
+            )
+
+        if name == "do3D" and value and self.drift_model != "Leeway":
+            self.logger.info("do3D is True so turning on vertical advection.")
             self.o.set_config("drift:vertical_advection", True)
+        elif name == "do3D" and value and self.drift_model == "Leeway":
+            self.logger.info(
+                "do3D is True but drift_model is Leeway so " "changing do3D to False."
+            )
+            self.do3D = False
 
         # Make sure vertical_mixing_timestep equals default value if vertical_mixing False
         if name in ["vertical_mixing", "vertical_mixing_timestep"]:
@@ -464,6 +488,25 @@ def __setattr_model__(self, name: str, value) -> None:
                 self.__dict__["stokes_drift"] = False
                 self.config_model["stokes_drift"]["value"] = False
 
+        # Add export variables for certain drift_model values
+        # drift_model is always set initially only
+        if name == "export_variables" and self.drift_model == "OpenOil":
+            oil_vars = [
+                "mass_oil",
+                "density",
+                "mass_evaporated",
+                "mass_dispersed",
+                "mass_biodegraded",
+                "viscosity",
+                "water_fraction",
+            ]
+            self.__dict__["export_variables"] += oil_vars
+            self.config_model["export_variables"]["value"] += oil_vars
+        elif name == "export_variables" and self.drift_model == "Leeway":
+            vars = ["object_type"]
+            self.__dict__["export_variables"] += vars
+            self.config_model["export_variables"]["value"] += vars
+
         self._update_config()
 
     def run_add_reader(
@@ -617,14 +660,13 @@ def run_add_reader(
 
                 drop_vars += [
                     "wetdry_mask_psi",
-                    "zeta",
                 ]
                 if self.ocean_model == "CIOFS":
 
                     loc_local = "/mnt/vault/ciofs/HINDCAST/ciofs_kerchunk.parq"
                     loc_remote = "http://xpublish-ciofs.srv.axds.co/datasets/ciofs_hindcast/zarr/"
 
-                elif self.ocean_model.upper() == "CIOFSOP":
+                elif self.ocean_model == "CIOFSOP":
 
                     standard_name_mapping.update(
                         {
@@ -674,7 +716,6 @@ def run_add_reader(
             # For NWGOA, need to calculate wetdry mask from a variable
             if self.ocean_model == "NWGOA" and not self.use_static_masks:
                 ds["wetdry_mask_rho"] = (~ds.zeta.isnull()).astype(int)
-                ds.drop_vars("zeta", inplace=True)
 
             # For CIOFSOP need to rename u/v to have "East" and "North" in the variable names
             # so they aren't rotated in the ROMS reader (the standard names have to be x/y not east/north)
@@ -683,7 +724,14 @@ def run_add_reader(
                 # grid = xr.open_dataset("/mnt/vault/ciofs/HINDCAST/nos.ciofs.romsgrid.nc")
                 # ds["angle"] = grid["angle"]
 
-            units_date = pd.Timestamp(ds.ocean_time.attrs["units"].split("since ")[1])
+            try:
+                units_date = pd.Timestamp(
+                    ds.ocean_time.attrs["units"].split("since ")[1]
+                )
+            except KeyError:  # for remote
+                units_date = pd.Timestamp(
+                    ds.ocean_time.encoding["units"].split("since ")[0]
+                )
             # use reader start time if not otherwise input
             if self.start_time is None:
                 self.logger.info("setting reader start_time as simulation start_time")
@@ -787,11 +835,18 @@ def run_drifters(self):
         }
 
         self.o._config = config_input_to_opendrift  # only OpenDrift config
+
+        output_file = (
+            self.output_file
+            or f"output-results_{datetime.datetime.utcnow():%Y-%m-%dT%H%M:%SZ}.nc"
+        )
+
         self.o.run(
             time_step=timedir * self.time_step,
+            time_step_output=self.time_step_output,
             steps=self.steps,
             export_variables=self.export_variables,
-            outfile=f"output-results_{datetime.datetime.utcnow():%Y-%m-%dT%H%M:%SZ}.nc",
+            outfile=output_file,
         )
 
         self.o._config = full_config  # reinstate config
@@ -902,6 +957,30 @@ def export_variables(self):
 
         return self.o.export_variables
 
+    def drift_model_config(self, ptm_level=[1, 2, 3], prefix=""):
+        """Show config for this drift model selection.
+
+        This shows all PTM-controlled parameters for the OpenDrift
+        drift model selected and their current values, at the selected ptm_level
+        of importance.
+
+        Parameters
+        ----------
+        ptm_level : int, list, optional
+            Options are 1, 2, 3, or lists of combinations. Use [1,2,3] for all.
+            Default is 1.
+        prefix : str, optional
+            prefix to search config for, only for OpenDrift parameters (not PTM).
+        """
+
+        return [
+            (key, value_dict["value"])
+            for key, value_dict in self.show_config(
+                substring=":", ptm_level=ptm_level, level=[1, 2, 3], prefix=prefix
+            ).items()
+            if "value" in value_dict
+        ]
+
     def get_configspec(self, prefix, substring, excludestring, level, ptm_level):
         """Copied from OpenDrift, then modified."""
 

particle_tracking_manager/plotting.py

@@ -3,7 +3,7 @@
 import opendrift
 
 
-def plot_dest(o):
+def plot_dest(o, filename):
     """This is copied from an opendrift example."""
 
     import cmocean
@@ -24,4 +24,6 @@ def plot_dest(o):
         vmin=0,
         vmax=vmax,
         cmap=cmap,
+        fast=True,
+        filename=filename,
     )

particle_tracking_manager/the_manager.py

@@ -71,10 +71,10 @@ class ParticleTrackingManager:
     run_forward : bool, optional
         True to run forward in time, False to run backward, by default True
     time_step : int, optional
-        Time step in seconds, options >0, <86400 (1 day in seconds), by default 3600
-    time_step_output : int, optional
-        How often to output model output, in seconds. Should be a multiple of time_step.
-        By default will take the value of time_step (this change occurs in the model).
+        Time step in seconds, options >0, <86400 (1 day in seconds), by default 300.
+    time_step_output : int, Timedelta, optional
+        How often to output model output. Should be a multiple of time_step.
+        By default 3600.
     steps : int, optional
         Number of time steps to run in simulation. Options >0.
         steps, end_time, or duration must be input by user. By default steps is 3 and
@@ -129,6 +129,8 @@ class ParticleTrackingManager:
         with a user-input ocean_model, you can drop the wetdry_mask_rho etc variables from the
         dataset before inputting to PTM. Setting this to True may save computation time but
         will be less accurate, especially in the tidal flat regions of the model.
+    output_file : Optional[str], optional
+        Name of output file to save, by default None. If None, default is set in the model.
 
     Notes
     -----
@@ -176,6 +178,7 @@ def __init__(
         do3D: bool = config_ptm["do3D"]["default"],
         vertical_mixing: bool = config_ptm["vertical_mixing"]["default"],
         use_static_masks: bool = config_ptm["use_static_masks"]["default"],
+        output_file: Optional[str] = config_ptm["output_file"]["default"],
         **kw,
     ) -> None:
         """Inputs necessary for any particle tracking."""
@@ -328,7 +331,7 @@ def __setattr__(self, name: str, value) -> None:
                     # this is not a user-defined option
                     if -180 < self.lon < 0:
                         self.__dict__["lon"] += 360
-                        self.config_ptm["lon"]["value"] = False
+                        self.config_ptm["lon"]["value"] += 360  # this isn't really used
 
             if name == "surface_only" and value:
                 self.logger.info(

particle_tracking_manager/the_manager_config.json

@@ -67,19 +67,19 @@
             "float",
             "datetime.timedelta"
         ],
-        "default": 3600,
+        "default": 300,
         "min": 1,
         "max": 86400,
         "units": "seconds",
         "description": "Interval between particles updates, in seconds or as timedelta.",
-        "ptm_level": 1
+        "ptm_level": 2
     },
     "time_step_output": {
         "type": [
             "float",
             "datetime.timedelta"
         ],
-        "default": "None",
+        "default": 3600,
         "min": 1,
         "max": 604800,
         "units": "seconds",
@@ -143,7 +143,7 @@
         "ptm_level": 1
     },
     "vertical_mixing": {
-        "default": false,
+        "default": true,
         "od_mapping": "drift:vertical_mixing",
         "ptm_level": 1
     },
@@ -162,5 +162,11 @@
         "default": false,
         "ptm_level": 2,
         "description": "Set to True to use static masks for known models instead of wetdry masks. If False, the masks are change in time."
+    },
+    "output_file": {
+        "type": "str",
+        "default": "None",
+        "description": "Name of file to write output to. If None, default name is used.",
+        "ptm_level": 3
     }
 }

tests/test_manager.py

@@ -97,10 +97,11 @@ def test_ocean_model_not_None(mock_reader_metadata):
         m.has_added_reader = True
 
 
+@pytest.mark.slow
 def test_parameter_passing():
     """make sure parameters passed into package make it to simulation runtime."""
 
-    ts = 2 * 3600
+    ts = 5
     diffmodel = "windspeed_Sundby1983"
     use_auto_landmask = True
     vertical_mixing = True
@@ -223,21 +224,26 @@ def test_input_too_many_end_of_simulation():
 
 
 def test_changing_end_of_simulation():
-    """change end_time, steps, and duration and make sure others are updated accordingly."""
+    """change end_time, steps, and duration
+
+    and make sure others are updated accordingly.
+    This accounts for the default time_step of 300 seconds.
+
+    """
 
     m = ptm.OpenDriftModel(start_time=pd.Timestamp("2000-1-1"))
     m.start_time = pd.Timestamp("2000-1-2")
     m.end_time = pd.Timestamp("2000-1-3")
-    assert m.steps == 24
+    assert m.steps == 288
     assert m.duration == pd.Timedelta("1 days 00:00:00")
 
     m.steps = 48
-    assert m.end_time == pd.Timestamp("2000-1-4")
-    assert m.duration == pd.Timedelta("2 days 00:00:00")
+    assert m.end_time == pd.Timestamp("2000-01-02 04:00:00")
+    assert m.duration == pd.Timedelta("0 days 04:00:00")
 
     m.duration = pd.Timedelta("2 days 12:00:00")
     assert m.end_time == pd.Timestamp("2000-01-04 12:00:00")
-    assert m.steps == 60
+    assert m.steps == 720
 
 
 class TestTheManager(unittest.TestCase):