Biogeochemistry features for expedition execution

src/virtualship/cli/_fetch.py

@@ -87,7 +87,7 @@ def _fetch(path: str | Path, username: str | None, password: str | None) -> None
 
     if (
         (
-            {"XBT", "CTD", "SHIP_UNDERWATER_ST"}
+            {"XBT", "CTD", "CDT_BGC", "SHIP_UNDERWATER_ST"}
             & set(instrument.name for instrument in instruments_in_schedule)
         )
         or ship_config.ship_underwater_st_config is not None
@@ -144,7 +144,6 @@ def _fetch(path: str | Path, username: str | None, password: str | None) -> None
             shutil.rmtree(download_folder)
             raise e
 
-        complete_download(download_folder)
         click.echo("Ship data download based on space-time region completed.")
 
     if InstrumentType.DRIFTER in instruments_in_schedule:
@@ -187,7 +186,6 @@ def _fetch(path: str | Path, username: str | None, password: str | None) -> None
             shutil.rmtree(download_folder)
             raise e
 
-        complete_download(download_folder)
         click.echo("Drifter data download based on space-time region completed.")
 
     if InstrumentType.ARGO_FLOAT in instruments_in_schedule:
@@ -235,9 +233,53 @@ def _fetch(path: str | Path, username: str | None, password: str | None) -> None
             shutil.rmtree(download_folder)
             raise e
 
-        complete_download(download_folder)
         click.echo("Argo_float data download based on space-time region completed.")
 
+    if InstrumentType.CTD_BGC in instruments_in_schedule:
+        print("CTD_BGC data will be downloaded. Please wait...")
+
+        ctd_bgc_download_dict = {
+            "o2data": {
+                "dataset_id": "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1D-m",
+                "variables": ["o2"],
+                "output_filename": "ctd_bgc_o2.nc",
+            },
+            "chlorodata": {
+                "dataset_id": "cmems_mod_glo_bgc-pft_anfc_0.25deg_P1D-m",
+                "variables": ["chl"],
+                "output_filename": "ctd_bgc_chloro.nc",
+            },
+        }
+
+        # Iterate over all datasets and download each based on space_time_region
+        try:
+            for dataset in ctd_bgc_download_dict.values():
+                copernicusmarine.subset(
+                    dataset_id=dataset["dataset_id"],
+                    variables=dataset["variables"],
+                    minimum_longitude=spatial_range.minimum_longitude - 3.0,
+                    maximum_longitude=spatial_range.maximum_longitude + 3.0,
+                    minimum_latitude=spatial_range.minimum_latitude - 3.0,
+                    maximum_latitude=spatial_range.maximum_latitude + 3.0,
+                    start_datetime=start_datetime,
+                    end_datetime=end_datetime + timedelta(days=21),
+                    minimum_depth=abs(1),
+                    maximum_depth=abs(spatial_range.maximum_depth),
+                    output_filename=dataset["output_filename"],
+                    output_directory=download_folder,
+                    username=username,
+                    password=password,
+                    overwrite=True,
+                    coordinates_selection_method="outside",
+                )
+        except InvalidUsernameOrPassword as e:
+            shutil.rmtree(download_folder)
+            raise e
+
+        click.echo("CTD_BGC data download based on space-time region completed.")
+
+    complete_download(download_folder)
+
 
 def _hash(s: str, *, length: int) -> str:
     """Create a hash of a string."""

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/ship_config.py

@@ -17,9 +17,10 @@
 
 
 class InstrumentType(Enum):
-    """Types of instruments."""
+    """Types of the instruments."""
 
     CTD = "CTD"
+    CTD_BGC = "CTD_BGC"
     DRIFTER = "DRIFTER"
     ARGO_FLOAT = "ARGO_FLOAT"
     XBT = "XBT"
@@ -80,6 +81,28 @@ def _validate_stationkeeping_time(cls, value: int | float | timedelta) -> timede
         return _validate_numeric_mins_to_timedelta(value)
 
 
+class CTD_BGCConfig(pydantic.BaseModel):
+    """Configuration for CTD_BGC instrument."""
+
+    stationkeeping_time: timedelta = pydantic.Field(
+        serialization_alias="stationkeeping_time_minutes",
+        validation_alias="stationkeeping_time_minutes",
+        gt=timedelta(),
+    )
+    min_depth_meter: float = pydantic.Field(le=0.0)
+    max_depth_meter: float = pydantic.Field(le=0.0)
+
+    model_config = pydantic.ConfigDict(populate_by_name=True)
+
+    @pydantic.field_serializer("stationkeeping_time")
+    def _serialize_stationkeeping_time(self, value: timedelta, _info):
+        return value.total_seconds() / 60.0
+
+    @pydantic.field_validator("stationkeeping_time", mode="before")
+    def _validate_stationkeeping_time(cls, value: int | float | timedelta) -> timedelta:
+        return _validate_numeric_mins_to_timedelta(value)
+
+
 class ShipUnderwaterSTConfig(pydantic.BaseModel):
     """Configuration for underwater ST."""
 
@@ -159,6 +182,13 @@ class ShipConfig(pydantic.BaseModel):
     If None, no CTDs can be cast.
     """
 
+    ctd_bgc_config: CTD_BGCConfig | None = None
+    """
+    CTD_BGC configuration.
+
+    If None, no BGC CTDs can be cast.
+    """
+
     ship_underwater_st_config: ShipUnderwaterSTConfig | None = None
     """
     Ship underwater salinity temperature measurementconfiguration.
@@ -239,6 +269,7 @@ def verify(self, schedule: Schedule) -> None:
             "DRIFTER",
             "XBT",
             "CTD",
+            "CTD_BGC",
         ]:  # TODO make instrument names consistent capitals or lowercase throughout codebase
             if hasattr(self, instrument.lower() + "_config") and not any(
                 instrument == schedule_instrument.name
@@ -248,6 +279,7 @@ def verify(self, schedule: Schedule) -> None:
                 setattr(self, instrument.lower() + "_config", None)
 
         # verify instruments in schedule have configuration
+        # TODO: the ConfigError message could be improved to explain that the **schedule** file has X instrument but the **ship_config** file does not
         for instrument in instruments_in_schedule:
             try:
                 InstrumentType(instrument)
@@ -266,6 +298,12 @@ def verify(self, schedule: Schedule) -> None:
                 raise ConfigError(
                     "Planning has a waypoint with CTD instrument, but configuration does not configure CTDs."
                 )
+            if instrument == InstrumentType.CTD_BGC and (
+                not hasattr(self, "ctd_bgc_config") or self.ctd_bgc_config is None
+            ):
+                raise ConfigError(
+                    "Planning has a waypoint with CTD_BGC instrument, but configuration does not configure CTD_BGCs."
+                )
             if instrument == InstrumentType.DRIFTER and (
                 not hasattr(self, "drifter_config") or self.drifter_config is None
             ):

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(