add ctd_bgc sampling capability to virtualship run execution

Author: j-atkins

src/virtualship/expedition/__init__.py

@@ -6,6 +6,7 @@
 from .ship_config import (
     ADCPConfig,
     ArgoFloatConfig,
+    CTD_BGCConfig,
     CTDConfig,
     DrifterConfig,
     ShipConfig,
@@ -17,6 +18,7 @@
     "ADCPConfig",
     "ArgoFloatConfig",
     "CTDConfig",
+    "CTD_BGCConfig",
     "DrifterConfig",
     "InputData",
     "InstrumentType",

src/virtualship/expedition/do_expedition.py

@@ -136,6 +136,7 @@ def _load_input_data(
         load_adcp=ship_config.adcp_config is not None,
         load_argo_float=ship_config.argo_float_config is not None,
         load_ctd=ship_config.ctd_config is not None,
+        load_ctd_bgc=ship_config.ctd_bgc_config is not None,
         load_drifter=ship_config.drifter_config is not None,
         load_xbt=ship_config.xbt_config is not None,
         load_ship_underwater_st=ship_config.ship_underwater_st_config is not None,

src/virtualship/expedition/input_data.py

@@ -15,6 +15,7 @@ class InputData:
     adcp_fieldset: FieldSet | None
     argo_float_fieldset: FieldSet | None
     ctd_fieldset: FieldSet | None
+    ctd_bgc_fieldset: FieldSet | None
     drifter_fieldset: FieldSet | None
     xbt_fieldset: FieldSet | None
     ship_underwater_st_fieldset: FieldSet | None
@@ -26,6 +27,7 @@ def load(
         load_adcp: bool,
         load_argo_float: bool,
         load_ctd: bool,
+        load_ctd_bgc: bool,
         load_drifter: bool,
         load_xbt: bool,
         load_ship_underwater_st: bool,
@@ -39,6 +41,7 @@ def load(
         :param load_adcp: Whether to load the ADCP fieldset.
         :param load_argo_float: Whether to load the argo float fieldset.
         :param load_ctd: Whether to load the CTD fieldset.
+        :param load_ctd_bgc: Whether to load the CTD BGC fieldset.
         :param load_drifter: Whether to load the drifter fieldset.
         :param load_ship_underwater_st: Whether to load the ship underwater ST fieldset.
         :returns: An instance of this class with loaded fieldsets.
@@ -51,6 +54,10 @@ def load(
             argo_float_fieldset = cls._load_argo_float_fieldset(directory)
         else:
             argo_float_fieldset = None
+        if load_ctd_bgc:
+            ctd_bgc_fieldset = cls._load_ctd_bgc_fieldset(directory)
+        else:
+            ctd_bgc_fieldset = None
         if load_adcp or load_ctd or load_ship_underwater_st or load_xbt:
             ship_fieldset = cls._load_ship_fieldset(directory)
         if load_adcp:
@@ -74,6 +81,7 @@ def load(
             adcp_fieldset=adcp_fieldset,
             argo_float_fieldset=argo_float_fieldset,
             ctd_fieldset=ctd_fieldset,
+            ctd_bgc_fieldset=ctd_bgc_fieldset,
             drifter_fieldset=drifter_fieldset,
             xbt_fieldset=xbt_fieldset,
             ship_underwater_st_fieldset=ship_underwater_st_fieldset,
@@ -122,6 +130,48 @@ def _load_ship_fieldset(cls, directory: Path) -> FieldSet:
 
         return fieldset
 
+    @classmethod
+    def _load_ctd_bgc_fieldset(cls, directory: Path) -> FieldSet:
+        filenames = {
+            "U": directory.joinpath("ship_uv.nc"),
+            "V": directory.joinpath("ship_uv.nc"),
+            "o2": directory.joinpath("ctd_bgc_o2.nc"),
+            "chl": directory.joinpath("ctd_bgc_chloro.nc"),
+        }
+        variables = {"U": "uo", "V": "vo", "o2": "o2", "chl": "chl"}
+        dimensions = {
+            "lon": "longitude",
+            "lat": "latitude",
+            "time": "time",
+            "depth": "depth",
+        }
+
+        fieldset = FieldSet.from_netcdf(
+            filenames, variables, dimensions, allow_time_extrapolation=True
+        )
+        fieldset.o2.interp_method = "linear_invdist_land_tracer"
+        fieldset.chl.interp_method = "linear_invdist_land_tracer"
+
+        # make depth negative
+        for g in fieldset.gridset.grids:
+            g.negate_depth()
+
+        # add bathymetry data
+        bathymetry_file = directory.joinpath("bathymetry.nc")
+        bathymetry_variables = ("bathymetry", "deptho")
+        bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"}
+        bathymetry_field = Field.from_netcdf(
+            bathymetry_file, bathymetry_variables, bathymetry_dimensions
+        )
+        # make depth negative
+        bathymetry_field.data = -bathymetry_field.data
+        fieldset.add_field(bathymetry_field)
+
+        # read in data already
+        fieldset.computeTimeChunk(0, 1)
+
+        return fieldset
+
     @classmethod
     def _load_drifter_fieldset(cls, directory: Path) -> FieldSet:
         filenames = {

src/virtualship/expedition/simulate_measurements.py

@@ -6,6 +6,7 @@
 from ..instruments.adcp import simulate_adcp
 from ..instruments.argo_float import simulate_argo_floats
 from ..instruments.ctd import simulate_ctd
+from ..instruments.ctd_bgc import simulate_ctd_bgc
 from ..instruments.drifter import simulate_drifters
 from ..instruments.ship_underwater_st import simulate_ship_underwater_st
 from ..instruments.xbt import simulate_xbt
@@ -75,6 +76,19 @@ def simulate_measurements(
             outputdt=timedelta(seconds=10),
         )
 
+    if len(measurements.ctd_bgcs) > 0:
+        print("Simulating BGC CTD casts.")
+        if ship_config.ctd_bgc_config is None:
+            raise RuntimeError("No configuration for CTD_BGC provided.")
+        if input_data.ctd_bgc_fieldset is None:
+            raise RuntimeError("No fieldset for CTD_BGC provided.")
+        simulate_ctd_bgc(
+            out_path=expedition_dir.joinpath("results", "ctd_bgc.zarr"),
+            fieldset=input_data.ctd_bgc_fieldset,
+            ctd_bgcs=measurements.ctd_bgcs,
+            outputdt=timedelta(seconds=10),
+        )
+
     if len(measurements.drifters) > 0:
         print("Simulating drifters")
         if ship_config.drifter_config is None:

src/virtualship/expedition/simulate_schedule.py

@@ -9,6 +9,7 @@
 
 from ..instruments.argo_float import ArgoFloat
 from ..instruments.ctd import CTD
+from ..instruments.ctd_bgc import CTD_BGC
 from ..instruments.drifter import Drifter
 from ..instruments.xbt import XBT
 from ..location import Location
@@ -42,6 +43,7 @@ class MeasurementsToSimulate:
     argo_floats: list[ArgoFloat] = field(default_factory=list, init=False)
     drifters: list[Drifter] = field(default_factory=list, init=False)
     ctds: list[CTD] = field(default_factory=list, init=False)
+    ctd_bgcs: list[CTD_BGC] = field(default_factory=list, init=False)
     xbts: list[XBT] = field(default_factory=list, init=False)
 
 
@@ -102,6 +104,7 @@ def simulate(self) -> ScheduleOk | ScheduleProblem:
             # check if waypoint was reached in time
             if waypoint.time is not None and self._time > waypoint.time:
                 print(
+                    # TODO: I think this should be wp_i + 1, not wp_i; otherwise it will be off by one
                     f"Waypoint {wp_i} could not be reached in time. Current time: {self._time}. Waypoint time: {waypoint.time}."
                 )
                 return ScheduleProblem(self._time, wp_i)
@@ -251,6 +254,15 @@ def _make_measurements(self, waypoint: Waypoint) -> timedelta:
                     )
                 )
                 time_costs.append(self._ship_config.ctd_config.stationkeeping_time)
+            elif instrument is InstrumentType.CTD_BGC:
+                self._measurements_to_simulate.ctd_bgcs.append(
+                    CTD_BGC(
+                        spacetime=Spacetime(self._location, self._time),
+                        min_depth=self._ship_config.ctd_bgc_config.min_depth_meter,
+                        max_depth=self._ship_config.ctd_bgc_config.max_depth_meter,
+                    )
+                )
+                time_costs.append(self._ship_config.ctd_bgc_config.stationkeeping_time)
             elif instrument is InstrumentType.DRIFTER:
                 self._measurements_to_simulate.drifters.append(
                     Drifter(