diff --git a/src/virtualship/cli/_fetch.py b/src/virtualship/cli/_fetch.py index 2e09dec2..f07c5b92 100644 --- a/src/virtualship/cli/_fetch.py +++ b/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.""" diff --git a/src/virtualship/expedition/__init__.py b/src/virtualship/expedition/__init__.py index 74403732..65b7cf11 100644 --- a/src/virtualship/expedition/__init__.py +++ b/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", diff --git a/src/virtualship/expedition/do_expedition.py b/src/virtualship/expedition/do_expedition.py index 97a66f7f..a4243361 100644 --- a/src/virtualship/expedition/do_expedition.py +++ b/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, diff --git a/src/virtualship/expedition/input_data.py b/src/virtualship/expedition/input_data.py index 554efb78..ca99f7e8 100644 --- a/src/virtualship/expedition/input_data.py +++ b/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 = { diff --git a/src/virtualship/expedition/ship_config.py b/src/virtualship/expedition/ship_config.py index 4f238b36..c0814f5b 100644 --- a/src/virtualship/expedition/ship_config.py +++ b/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 ): diff --git a/src/virtualship/expedition/simulate_measurements.py b/src/virtualship/expedition/simulate_measurements.py index bf28c989..c7cd8746 100644 --- a/src/virtualship/expedition/simulate_measurements.py +++ b/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: diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index df36b8b2..596d16a2 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/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( diff --git a/src/virtualship/instruments/ctd_bgc.py b/src/virtualship/instruments/ctd_bgc.py new file mode 100644 index 00000000..cb218e3a --- /dev/null +++ b/src/virtualship/instruments/ctd_bgc.py @@ -0,0 +1,143 @@ +"""CTD_BGC instrument.""" + +from dataclasses import dataclass +from datetime import timedelta +from pathlib import Path + +import numpy as np +from parcels import FieldSet, JITParticle, ParticleSet, Variable + +from ..spacetime import Spacetime + + +@dataclass +class CTD_BGC: + """Configuration for a single BGC CTD.""" + + spacetime: Spacetime + min_depth: float + max_depth: float + + +_CTD_BGCParticle = JITParticle.add_variables( + [ + Variable("o2", dtype=np.float32, initial=np.nan), + Variable("chl", dtype=np.float32, initial=np.nan), + Variable("raising", dtype=np.int8, initial=0.0), # bool. 0 is False, 1 is True. + Variable("max_depth", dtype=np.float32), + Variable("min_depth", dtype=np.float32), + Variable("winch_speed", dtype=np.float32), + ] +) + + +def _sample_o2(particle, fieldset, time): + particle.o2 = fieldset.o2[time, particle.depth, particle.lat, particle.lon] + + +def _sample_chlorophyll(particle, fieldset, time): + particle.chl = fieldset.chl[time, particle.depth, particle.lat, particle.lon] + + +def _ctd_bgc_cast(particle, fieldset, time): + # lowering + if particle.raising == 0: + particle_ddepth = -particle.winch_speed * particle.dt + if particle.depth + particle_ddepth < particle.max_depth: + particle.raising = 1 + particle_ddepth = -particle_ddepth + # raising + else: + particle_ddepth = particle.winch_speed * particle.dt + if particle.depth + particle_ddepth > particle.min_depth: + particle.delete() + + +def simulate_ctd_bgc( + fieldset: FieldSet, + out_path: str | Path, + ctd_bgcs: list[CTD_BGC], + outputdt: timedelta, +) -> None: + """ + Use Parcels to simulate a set of BGC CTDs in a fieldset. + + :param fieldset: The fieldset to simulate the BGC CTDs in. + :param out_path: The path to write the results to. + :param ctds: A list of BGC CTDs to simulate. + :param outputdt: Interval which dictates the update frequency of file output during simulation + :raises ValueError: Whenever provided BGC CTDs, fieldset, are not compatible with this function. + """ + WINCH_SPEED = 1.0 # sink and rise speed in m/s + DT = 10.0 # dt of CTD simulation integrator + + if len(ctd_bgcs) == 0: + print( + "No BGC CTDs provided. Parcels currently crashes when providing an empty particle set, so no BGC CTD simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + + fieldset_starttime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[0]) + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + + # deploy time for all ctds should be later than fieldset start time + if not all( + [ + np.datetime64(ctd_bgc.spacetime.time) >= fieldset_starttime + for ctd_bgc in ctd_bgcs + ] + ): + raise ValueError("BGC CTD deployed before fieldset starts.") + + # depth the bgc ctd will go to. shallowest between bgc ctd max depth and bathymetry. + max_depths = [ + max( + ctd_bgc.max_depth, + fieldset.bathymetry.eval( + z=0, + y=ctd_bgc.spacetime.location.lat, + x=ctd_bgc.spacetime.location.lon, + time=0, + ), + ) + for ctd_bgc in ctd_bgcs + ] + + # CTD depth can not be too shallow, because kernel would break. + # This shallow is not useful anyway, no need to support. + if not all([max_depth <= -DT * WINCH_SPEED for max_depth in max_depths]): + raise ValueError( + f"BGC CTD max_depth or bathymetry shallower than maximum {-DT * WINCH_SPEED}" + ) + + # define parcel particles + ctd_bgc_particleset = ParticleSet( + fieldset=fieldset, + pclass=_CTD_BGCParticle, + lon=[ctd_bgc.spacetime.location.lon for ctd_bgc in ctd_bgcs], + lat=[ctd_bgc.spacetime.location.lat for ctd_bgc in ctd_bgcs], + depth=[ctd_bgc.min_depth for ctd_bgc in ctd_bgcs], + time=[ctd_bgc.spacetime.time for ctd_bgc in ctd_bgcs], + max_depth=max_depths, + min_depth=[ctd_bgc.min_depth for ctd_bgc in ctd_bgcs], + winch_speed=[WINCH_SPEED for _ in ctd_bgcs], + ) + + # define output file for the simulation + out_file = ctd_bgc_particleset.ParticleFile(name=out_path, outputdt=outputdt) + + # execute simulation + ctd_bgc_particleset.execute( + [_sample_o2, _sample_chlorophyll, _ctd_bgc_cast], + endtime=fieldset_endtime, + dt=DT, + verbose_progress=False, + output_file=out_file, + ) + + # there should be no particles left, as they delete themselves when they resurface + if len(ctd_bgc_particleset.particledata) != 0: + raise ValueError( + "Simulation ended before BGC CTD resurfaced. This most likely means the field time dimension did not match the simulation time span." + ) diff --git a/src/virtualship/static/schedule.yaml b/src/virtualship/static/schedule.yaml index d15736d0..7cb39423 100644 --- a/src/virtualship/static/schedule.yaml +++ b/src/virtualship/static/schedule.yaml @@ -12,6 +12,7 @@ space_time_region: waypoints: - instrument: - CTD + - CTD_BGC location: latitude: 0 longitude: 0 diff --git a/src/virtualship/static/ship_config.yaml b/src/virtualship/static/ship_config.yaml index 5066e38d..34d6c6ea 100644 --- a/src/virtualship/static/ship_config.yaml +++ b/src/virtualship/static/ship_config.yaml @@ -14,6 +14,10 @@ ctd_config: max_depth_meter: -2000.0 min_depth_meter: -11.0 stationkeeping_time_minutes: 20.0 +ctd_bgc_config: + max_depth_meter: -2000.0 + min_depth_meter: -11.0 + stationkeeping_time_minutes: 20.0 drifter_config: depth_meter: 0.0 lifetime_minutes: 60480.0 diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py index 1d3203de..6dbbb49c 100644 --- a/src/virtualship/utils.py +++ b/src/virtualship/utils.py @@ -156,6 +156,7 @@ def mfp_to_yaml(coordinates_file_path: str, yaml_output_path: str): # noqa: D41 instrument_max_depths = { "XBT": 2000, "CTD": 5000, + "CTD_BGC": 5000, "DRIFTER": 1, "ARGO_FLOAT": 2000, } diff --git a/tests/expedition/expedition_dir/input_data/ctd_bgc_chloro.nc b/tests/expedition/expedition_dir/input_data/ctd_bgc_chloro.nc new file mode 100644 index 00000000..efdc51d3 Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/ctd_bgc_chloro.nc differ diff --git a/tests/expedition/expedition_dir/input_data/ctd_bgc_o2.nc b/tests/expedition/expedition_dir/input_data/ctd_bgc_o2.nc new file mode 100644 index 00000000..d552f8bb Binary files /dev/null and b/tests/expedition/expedition_dir/input_data/ctd_bgc_o2.nc differ diff --git a/tests/expedition/expedition_dir/ship_config.yaml b/tests/expedition/expedition_dir/ship_config.yaml index 09f40c0a..1bae9d1d 100644 --- a/tests/expedition/expedition_dir/ship_config.yaml +++ b/tests/expedition/expedition_dir/ship_config.yaml @@ -14,6 +14,10 @@ ctd_config: max_depth_meter: -2000.0 min_depth_meter: -11.0 stationkeeping_time_minutes: 20.0 +ctd_bgc_config: + max_depth_meter: -2000.0 + min_depth_meter: -11.0 + stationkeeping_time_minutes: 20.0 drifter_config: depth_meter: 0.0 lifetime_minutes: 40320.0 diff --git a/tests/expedition/test_ship_config.py b/tests/expedition/test_ship_config.py index 44bfd524..920058d2 100644 --- a/tests/expedition/test_ship_config.py +++ b/tests/expedition/test_ship_config.py @@ -50,6 +50,12 @@ def ship_config_no_ctd(ship_config): return ship_config +@pytest.fixture +def ship_config_no_ctd_bgc(ship_config): + delattr(ship_config, "ctd_bgc_config") + return ship_config + + @pytest.fixture def ship_config_no_argo_float(ship_config): delattr(ship_config, "argo_float_config") @@ -91,6 +97,12 @@ def test_verify_ship_config_no_instrument(ship_config, schedule_no_xbt) -> None: "Planning has a waypoint with CTD instrument, but configuration does not configure CTD.", id="ShipConfigNoCTD", ), + pytest.param( + "ship_config_no_ctd_bgc", + ConfigError, + "Planning has a waypoint with CTD_BGC instrument, but configuration does not configure CTD_BGCs.", + id="ShipConfigNoCTD_BGC", + ), pytest.param( "ship_config_no_argo_float", ConfigError, diff --git a/tests/expedition/test_simulate_schedule.py b/tests/expedition/test_simulate_schedule.py index 6381b805..8c42097b 100644 --- a/tests/expedition/test_simulate_schedule.py +++ b/tests/expedition/test_simulate_schedule.py @@ -53,4 +53,5 @@ def test_time_in_minutes_in_ship_schedule() -> None: ship_config = ShipConfig.from_yaml("expedition_dir/ship_config.yaml") assert ship_config.adcp_config.period == timedelta(minutes=5) assert ship_config.ctd_config.stationkeeping_time == timedelta(minutes=20) + assert ship_config.ctd_bgc_config.stationkeeping_time == timedelta(minutes=20) assert ship_config.ship_underwater_st_config.period == timedelta(minutes=5) diff --git a/tests/instruments/test_ctd_bgc.py b/tests/instruments/test_ctd_bgc.py new file mode 100644 index 00000000..8083c580 --- /dev/null +++ b/tests/instruments/test_ctd_bgc.py @@ -0,0 +1,137 @@ +""" +Test the simulation of CTD_BGC instruments. + +Fields are kept static over time and time component of CTD_BGC measurements is not tested because it's tricky to provide expected measurements. +""" + +import datetime +from datetime import timedelta + +import numpy as np +import xarray as xr +from parcels import Field, FieldSet + +from virtualship import Location, Spacetime +from virtualship.instruments.ctd_bgc import CTD_BGC, simulate_ctd_bgc + + +def test_simulate_ctd_bgcs(tmpdir) -> None: + # arbitrary time offset for the dummy fieldset + base_time = datetime.datetime.strptime("1950-01-01", "%Y-%m-%d") + + # where to cast CTD_BGCs + ctd_bgcs = [ + CTD_BGC( + spacetime=Spacetime( + location=Location(latitude=0, longitude=1), + time=base_time + datetime.timedelta(hours=0), + ), + min_depth=0, + max_depth=float("-inf"), + ), + CTD_BGC( + spacetime=Spacetime( + location=Location(latitude=1, longitude=0), + time=base_time, + ), + min_depth=0, + max_depth=float("-inf"), + ), + ] + + # expected observations for ctd_bgcs at surface and at maximum depth + ctd_bgc_exp = [ + { + "surface": { + "o2": 9, + "chl": 10, + "lat": ctd_bgcs[0].spacetime.location.lat, + "lon": ctd_bgcs[0].spacetime.location.lon, + }, + "maxdepth": { + "o2": 11, + "chl": 12, + "lat": ctd_bgcs[0].spacetime.location.lat, + "lon": ctd_bgcs[0].spacetime.location.lon, + }, + }, + { + "surface": { + "o2": 9, + "chl": 10, + "lat": ctd_bgcs[1].spacetime.location.lat, + "lon": ctd_bgcs[1].spacetime.location.lon, + }, + "maxdepth": { + "o2": 11, + "chl": 12, + "lat": ctd_bgcs[1].spacetime.location.lat, + "lon": ctd_bgcs[1].spacetime.location.lon, + }, + }, + ] + + # create fieldset based on the expected observations + # indices are time, depth, latitude, longitude + u = np.zeros((2, 2, 2, 2)) + v = np.zeros((2, 2, 2, 2)) + o2 = np.zeros((2, 2, 2, 2)) # o2 field + chl = np.zeros((2, 2, 2, 2)) # chl field + + o2[:, 1, 0, 1] = ctd_bgc_exp[0]["surface"]["o2"] + o2[:, 0, 0, 1] = ctd_bgc_exp[0]["maxdepth"]["o2"] + o2[:, 1, 1, 0] = ctd_bgc_exp[1]["surface"]["o2"] + o2[:, 0, 1, 0] = ctd_bgc_exp[1]["maxdepth"]["o2"] + + chl[:, 1, 0, 1] = ctd_bgc_exp[0]["surface"]["chl"] + chl[:, 0, 0, 1] = ctd_bgc_exp[0]["maxdepth"]["chl"] + chl[:, 1, 1, 0] = ctd_bgc_exp[1]["surface"]["chl"] + chl[:, 0, 1, 0] = ctd_bgc_exp[1]["maxdepth"]["chl"] + + fieldset = FieldSet.from_data( + {"V": v, "U": u, "o2": o2, "chl": chl}, + { + "time": [ + np.datetime64(base_time + datetime.timedelta(hours=0)), + np.datetime64(base_time + datetime.timedelta(hours=1)), + ], + "depth": [-1000, 0], + "lat": [0, 1], + "lon": [0, 1], + }, + ) + fieldset.add_field(Field("bathymetry", [-1000], lon=0, lat=0)) + + # perform simulation + out_path = tmpdir.join("out.zarr") + + simulate_ctd_bgc( + ctd_bgcs=ctd_bgcs, + fieldset=fieldset, + out_path=out_path, + outputdt=timedelta(seconds=10), + ) + + # test if output is as expected + results = xr.open_zarr(out_path) + + assert len(results.trajectory) == len(ctd_bgcs) + + for ctd_i, (traj, exp_bothloc) in enumerate( + zip(results.trajectory, ctd_bgc_exp, strict=True) + ): + obs_surface = results.sel(trajectory=traj, obs=0) + min_index = np.argmin(results.sel(trajectory=traj)["z"].data) + obs_maxdepth = results.sel(trajectory=traj, obs=min_index) + + for obs, loc in [ + (obs_surface, "surface"), + (obs_maxdepth, "maxdepth"), + ]: + exp = exp_bothloc[loc] + for var in ["o2", "chl", "lat", "lon"]: + obs_value = obs[var].values.item() + exp_value = exp[var] + assert np.isclose(obs_value, exp_value), ( + f"Observation incorrect {ctd_i=} {loc=} {var=} {obs_value=} {exp_value=}." + )