diff --git a/docs/conf.py b/docs/conf.py index 2460915d..27353cbb 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -75,6 +75,10 @@ "user-guide/assignments/Research_Proposal_only": "user-guide/_images/MFP.jpg", "user-guide/assignments/Virtualship_research_proposal": "user-guide/_images/AnnaWeber.jpeg", "user-guide/assignments/sciencecommunication_assignment": "user-guide/_images/marine_ss.jpg", - "user-guide/assignments/Sail_the_ship": "user-guide/_images/vessel.jpg", + "user-guide/assignments/Sail_the_ship": "user-guide/_images/freepik_research_vessel.jpg", + "user-guide/assignments/Code_of_conduct": "user-guide/_images/freepik_code_of_conduct.jpg", } + +sphinx_gallery_conf = {"default_thumb_file": "_static/virtual_ship_logo.png"} + nbsphinx_execute = "never" diff --git a/docs/user-guide/_images/ILOs.jpg b/docs/user-guide/_images/ILOs.jpg new file mode 100644 index 00000000..81c09523 Binary files /dev/null and b/docs/user-guide/_images/ILOs.jpg differ diff --git a/docs/user-guide/_images/freepik_code_of_conduct.jpg b/docs/user-guide/_images/freepik_code_of_conduct.jpg new file mode 100644 index 00000000..ae1707bd Binary files /dev/null and b/docs/user-guide/_images/freepik_code_of_conduct.jpg differ diff --git a/docs/user-guide/_images/freepik_research_vessel.jpg b/docs/user-guide/_images/freepik_research_vessel.jpg new file mode 100644 index 00000000..71163e06 Binary files /dev/null and b/docs/user-guide/_images/freepik_research_vessel.jpg differ diff --git a/docs/user-guide/_images/vessel.jpg b/docs/user-guide/_images/vessel.jpg deleted file mode 100644 index 7c0d54f0..00000000 Binary files a/docs/user-guide/_images/vessel.jpg and /dev/null differ diff --git a/docs/user-guide/assignments/Code_of_conduct.ipynb b/docs/user-guide/assignments/Code_of_conduct.ipynb new file mode 100644 index 00000000..470add25 --- /dev/null +++ b/docs/user-guide/assignments/Code_of_conduct.ipynb @@ -0,0 +1,104 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Code of conduct\n", + "_As used during the Dynamical Oceanography 2024/25 course at Utrecht University_" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Please read this code of conduct, fill out the [parts between brackets] after consultation with your group, remove these italic sections, and then all individually upload a copy._ \n", + "\n", + "_Note that the Code of Conduct applies to the actual collaboration within your group in the Virtual Ship Assignment; so is not about how you would work on the ship._" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This code of conduct has been decided by [NAMES] for their group project during the Virtual Ship Classroom and is enacted from [DD-MM-YYYY]. The procedure is based on [this exercise by Aurelia Moser](http://aureliamoser.com/aaas-guides/conduct/index.html)\n", + "\n", + "Everyone taking part in the course and group discussions (mentors, helpers, coordinators, and learners) is required to conform to the following Code of Conduct. Coordinators will oversee adherence to this code throughout the course.\n", + "\n", + "Characteristics we value: [FILL IN]\n", + "\n", + "Behaviors we encourage: [FILL IN]\n", + "\n", + "Behaviors we discourage: [FILL IN]\n", + "\n", + "The way we redistribute the grade: [FILL IN]\n", + "\n", + "Choose for example “Everyone in the team gets the same grade”, “Everyone can propose a bonus point to one other person, which is then subtracted from other team member’s grades”, “We decide as a group who gets up to x bonus/penalty points”, or any other way to redistribute the grade, as long as the average of the group is not affected. \n", + "\n", + "We make each others feel safe and supported by: [FILL IN]\n", + "\n", + "How to report an issue, should someone violate the code?\n", + "1.\tContact the course instructor by private message, or in person. All communication will be treated as confidential.\n", + "2.\tIf for any reason you don’t want to do 1, you can contact the university Academic Integrity Counsellor." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Community Participation Guidelines" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Apart from the code of conduct, all participants in this course must follow these general guidelines.\n", + "Participation \n", + "\n", + "When participating in Dynamical Oceanography, respect the Utrecht University guideline for academic and scientific integrity and code of conduct. These guidelines cover our behaviour as participants, mentors, experts, staff, volunteers, and anyone else involved in making this course possible.\n", + "\n", + "How to treat each other\n", + "\n", + "To create a collaborative and inviting learning environment, we also emphasise certain values in how we treat each other:\n", + "*\tBe respectful and value each other’s ideas, styles and viewpoints.\n", + "*\tBe direct but professional; we cannot withhold hard truths.\n", + "*\tBe inclusive and help new perspectives be heard.\n", + "*\tAppreciate and accommodate our many cultural practices, attitudes and beliefs.\n", + "*\tBe open to learning from others.\n", + "*\tLead by example and match your actions with your words.\n", + "\n", + "The following will not be tolerated during the activities related to this course:\n", + "*\tviolence and threats of violence.\n", + "*\tpersonal attacks; derogatory language.\n", + "*\tunwelcome sexual attention or physical contact.\n", + "*\tdisruptive behavior.\n", + "*\tinfluencing unacceptable behavior.\n", + "\n", + "### Inclusion and Diversity\n", + "\n", + "We welcome contributions from everyone as long as they interact constructively with our community, including, but not limited to people of varied age, culture, ethnicity, gender, gender-identity, language, race, sexual orientation, geographical location and religious views.\n", + "\n", + "Raising Issues\n", + "\n", + "If you believe you‘re experiencing practices which don‘t meet the above policies, or if you feel you are being harassed in any way, please immediately contact the course coordinator, or the university Academic Integrity Counsellor.\n", + "\n", + "The course coordinator reserves the right to refuse admission to anyone violating these policies, and/or take further action including reporting to the responsible figure for academic integrity at the department.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Thumbnail image](_images/freepik_code_of_conduct.jpg)" + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/user-guide/assignments/Sail_the_ship.ipynb b/docs/user-guide/assignments/Sail_the_ship.ipynb index 3f62d06a..8bfcd463 100644 --- a/docs/user-guide/assignments/Sail_the_ship.ipynb +++ b/docs/user-guide/assignments/Sail_the_ship.ipynb @@ -177,6 +177,13 @@ "![Our ocean provides](https://seewhatgrows.org/wp-content/uploads/2019/06/our-ocean.jpg)\n", "![10 ways to help our ocean](https://oceanservice.noaa.gov/ocean/help-ocean.jpg)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Thumbnail image](_images/freepik_research_vessel.jpg)" + ] } ], "metadata": { diff --git a/docs/user-guide/assignments/index.md b/docs/user-guide/assignments/index.md index 2b37970f..0e0db109 100644 --- a/docs/user-guide/assignments/index.md +++ b/docs/user-guide/assignments/index.md @@ -5,9 +5,9 @@ maxdepth: 1 caption: Assignments --- +Virtualship_research_proposal.ipynb Research_proposal_intro.ipynb Research_Proposal_only.ipynb -Virtualship_research_proposal.ipynb sciencecommunication_assignment.ipynb ``` diff --git a/docs/user-guide/index.md b/docs/user-guide/index.md index 6d950781..2578fbea 100644 --- a/docs/user-guide/index.md +++ b/docs/user-guide/index.md @@ -4,6 +4,7 @@ :maxdepth: 1 quickstart +teacher-content/index tutorials/index assignments/index ``` diff --git a/docs/user-guide/teacher-content/ILOs.ipynb b/docs/user-guide/teacher-content/ILOs.ipynb new file mode 100644 index 00000000..f17cb126 --- /dev/null +++ b/docs/user-guide/teacher-content/ILOs.ipynb @@ -0,0 +1,69 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "84d1c60b", + "metadata": {}, + "source": [ + "# Aligned learning outcomes\n", + "\n", + "_Aligned using the revised version of Blooms taxonomy (Anderson et al. 2001)_\n", + "\n", + "1. Identify practical challenges of planning sea-based research, such as:\n", + "\n", + "- Availability and capabilities of research vessels and equipment, e.g. ice breaker vessel needed, infrastructure for equipment available\n", + "- General logistics, e.g. transportation of container labs and other equipment\n", + "- Requesting permits\n", + " - Travel permits for staff\n", + " - Diplomatic clearance (EEZ)\n", + " - Customs documents\n", + "- Safety trainings (beforehand and onboard)\n", + "- Preparing, testing and transporting equipment\n", + "\n", + "2. Be aware of (describe) matters during sea-based research, such as:\n", + "\n", + "- Contingency due to weather conditions and other unforeseen circumstances\n", + "- The need to be flexible and patient\n", + "- Taking notes, recording coordinates, stations, casts\n", + "- Measurement protocols\n", + "- Communication with crew, e.g. captain, winch driver\n", + "- Living at sea, e.g. sea legs, seasickness, relaxation, water/energy/waste footprint\n", + "\n", + "3. Give examples of uncertainty and variability of sea-based observations, such as:\n", + "\n", + "- Likely errors or wrong values, and seeing the data through the noise\n", + "- Time and date coordinates, e.g. UTC, local time, ship time\n", + "- Variability and limitations of measurements\n", + "- Different scales of variability, e.g. tidal, seasonal\n", + "- ‘Synopticity’ of measurements (i.e. not possible to measure with one ship at multiple locations at the same time; therefore mixing time and space variability)\n", + "- Instrument specific qualities, e.g. ADCP reflection\n", + "\n", + "4. Learn how to analyse and interpret sea-based observations. Be able to:\n", + "\n", + "- Read in different data formats, e.g. nc, csv, zarr\n", + "- Manipulate arrays and other data structures in python\n", + "- Create and interpret density/TS profiles\n", + "- Create spatial temperature/salinity/velocity maps, do GIS things\n", + "- Interpret plots and deduce ocean circulation\n", + "- Report findings\n", + "- Contribute to a cruise report with pictures or plots\n", + "\n", + "5. Plan a research cruise. Be able to:\n", + "\n", + "- Formulate research questions\n", + "- Select sampling sites\n", + "- Consider the duration of deployments, e.g. time a CTD station takes\n", + "- Decide on on-site measurement order\n", + "- Plan rest/working shifts for scientific personnel\n", + "- Deal with contingency, e.g. consider what sites/measurements are essential, which are ok to skip" + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/user-guide/teacher-content/index.md b/docs/user-guide/teacher-content/index.md new file mode 100644 index 00000000..ba1596b6 --- /dev/null +++ b/docs/user-guide/teacher-content/index.md @@ -0,0 +1,26 @@ +# Teacher content + +VirtualShip is used as part of the VirtualShip Classroom, that combines authentic tools with VR to create a virtual fieldwork experience and allows you to teach about sea-based research from your regular classroom. + +All VirtualShip Classroom (VSC) material is open under an MIT licence and freely available! You can use the VSC to teach anything from a 4 hour masterclass up to an open assignment of more than 40 hours. Example assignments are available below and please feel free to customize anything offline or [contribute](../../contributing/index.md) to the assignments provided here. + +Our core learning outcome is for students to appreciate the difficulty of measuring the ocean, making the VSC a valuable tool for students that will go into fieldwork, modelling or anything else. Additional Intended Learning Outcomes [(ILOs)](ILOs.md) were formulated in collaboration with ocean scientists from the Royal Netherlands Institute for Sea Research [(NIOZ)](https://www.nioz.nl/en) and Utrecht University [(UU)](https://www.uu.nl/en). + +The 360 videos are available on our YouTube channel [**_@VirtualShip Classroom_**](https://www.youtube.com/@VirtualShipClassroom) and in high resolution upon request. The videos can be viewed in the classroom or at home using: + +- VR Headsets – Wear a headset and look around naturally by moving your head. +- Smartphones/Tablets – Move your device or swipe the screen to explore. +- PC/Mac Browsers – Click and drag with your mouse to look around. + +The VSC design focuses on creating didactically sound, authentic learning experiences grounded in established learning theories in science education, such as constructivism [(Piaget 1954)](https://doi.org/10.4324/9781315009650) and constructionism [(Papert 1980)](https://worrydream.com/refs/Papert_1980_-_Mindstorms,_1st_ed.pdf). By integrating realistic tasks and a gamified narrative approach within Jupyter notebooks, students learn within a digital replica of the real world, constructing knowledge through ‘learning by doing’ and ‘trial and error’ as they explore oceanography concepts, research methods, and analysis tools. + +We evaluated in several (under)graduate courses and find that the VirtualShip Classroom is highly engaging, and students report on enhanced confidence and knowledge [(Daniels et al. 2025)](https://current-journal.com/articles/10.5334/cjme.121). + +```{nbgallery} +--- +maxdepth: 1 +caption: Teaching material +--- + +ILOs.md +``` diff --git a/src/virtualship/cli/_creds.py b/src/virtualship/cli/_creds.py index 9f1d2435..d9169ec4 100644 --- a/src/virtualship/cli/_creds.py +++ b/src/virtualship/cli/_creds.py @@ -1,3 +1,8 @@ +# + + +# TODO: TO DELETE?! + from __future__ import annotations from pathlib import Path diff --git a/src/virtualship/cli/_fetch.py b/src/virtualship/cli/_fetch.py deleted file mode 100644 index 19f000a7..00000000 --- a/src/virtualship/cli/_fetch.py +++ /dev/null @@ -1,563 +0,0 @@ -from __future__ import annotations - -import hashlib -import shutil -from datetime import datetime, timedelta -from pathlib import Path -from typing import TYPE_CHECKING - -import numpy as np -from pydantic import BaseModel - -from virtualship.errors import CopernicusCatalogueError, IncompleteDownloadError -from virtualship.utils import ( - _dump_yaml, - _generic_load_yaml, - _get_expedition, -) - -if TYPE_CHECKING: - from virtualship.models import SpaceTimeRegion - -import click -import copernicusmarine -from copernicusmarine.core_functions.credentials_utils import InvalidUsernameOrPassword - -import virtualship.cli._creds as creds -from virtualship.utils import EXPEDITION - -DOWNLOAD_METADATA = "download_metadata.yaml" - - -def _fetch(path: str | Path, username: str | None, password: str | None) -> None: - """ - Download input data for an expedition. - - Entrypoint for the tool to download data based on space-time region provided in the - schedule file. Data is downloaded from Copernicus Marine, credentials for which can be - obtained via registration: https://data.marine.copernicus.eu/register . Credentials can - be provided on prompt, via command line arguments, or via a YAML config file. Run - `virtualship fetch` on an expedition for more info. - """ - from virtualship.models import InstrumentType - - if sum([username is None, password is None]) == 1: - raise ValueError("Both username and password must be provided when using CLI.") - - path = Path(path) - - data_folder = path / "data" - data_folder.mkdir(exist_ok=True) - - expedition = _get_expedition(path) - - expedition.schedule.verify( - expedition.ship_config.ship_speed_knots, - input_data=None, - check_space_time_region=True, - ignore_missing_fieldsets=True, - ) - - space_time_region_hash = get_space_time_region_hash( - expedition.schedule.space_time_region - ) - - existing_download = get_existing_download(data_folder, space_time_region_hash) - if existing_download is not None: - click.echo( - f"Data download for space-time region already completed ('{existing_download}')." - ) - return - - creds_path = path / creds.CREDENTIALS_FILE - username, password = creds.get_credentials_flow(username, password, creds_path) - - # Extract space_time_region details from the schedule - spatial_range = expedition.schedule.space_time_region.spatial_range - time_range = expedition.schedule.space_time_region.time_range - start_datetime = time_range.start_time - end_datetime = time_range.end_time - instruments_in_schedule = expedition.schedule.get_instruments() - - # Create download folder and set download metadata - download_folder = data_folder / hash_to_filename(space_time_region_hash) - download_folder.mkdir() - DownloadMetadata(download_complete=False).to_yaml( - download_folder / DOWNLOAD_METADATA - ) - shutil.copyfile(path / EXPEDITION, download_folder / EXPEDITION) - - # data download - if ( - ( - {"XBT", "CTD", "CDT_BGC", "SHIP_UNDERWATER_ST"} - & set(instrument.name for instrument in instruments_in_schedule) - ) - or expedition.instruments_config.ship_underwater_st_config is not None - or expedition.instruments_config.adcp_config is not None - ): - print("Ship data will be downloaded. Please wait...") - - phys_product_id = select_product_id( - physical=True, - schedule_start=start_datetime, - schedule_end=end_datetime, - username=username, - password=password, - ) - - # Define all ship datasets to download, including bathymetry - download_dict = { - "Bathymetry": { - "dataset_id": "cmems_mod_glo_phy_my_0.083deg_static", - "variables": ["deptho"], - "output_filename": "bathymetry.nc", - }, - "UVdata": { - "dataset_id": phys_product_id, - "variables": ["uo", "vo"], - "output_filename": "ship_uv.nc", - }, - "Sdata": { - "dataset_id": phys_product_id, - "variables": ["so"], - "output_filename": "ship_s.nc", - }, - "Tdata": { - "dataset_id": phys_product_id, - "variables": ["thetao"], - "output_filename": "ship_t.nc", - }, - } - - # Iterate over all datasets and download each based on space_time_region - try: - for dataset in download_dict.values(): - copernicusmarine.subset( - dataset_id=dataset["dataset_id"], - variables=dataset["variables"], - minimum_longitude=spatial_range.minimum_longitude, - maximum_longitude=spatial_range.maximum_longitude, - minimum_latitude=spatial_range.minimum_latitude, - maximum_latitude=spatial_range.maximum_latitude, - start_datetime=start_datetime, - end_datetime=end_datetime, - minimum_depth=abs(spatial_range.minimum_depth), - 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("Ship data download based on space-time region completed.") - - if InstrumentType.DRIFTER in instruments_in_schedule: - print("Drifter data will be downloaded. Please wait...") - drifter_download_dict = { - "UVdata": { - "dataset_id": phys_product_id, - "variables": ["uo", "vo"], - "output_filename": "drifter_uv.nc", - }, - "Tdata": { - "dataset_id": phys_product_id, - "variables": ["thetao"], - "output_filename": "drifter_t.nc", - }, - } - - # Iterate over all datasets and download each based on space_time_region - try: - for dataset in drifter_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(1), - 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("Drifter data download based on space-time region completed.") - - if InstrumentType.ARGO_FLOAT in instruments_in_schedule: - print("Argo float data will be downloaded. Please wait...") - argo_download_dict = { - "UVdata": { - "dataset_id": phys_product_id, - "variables": ["uo", "vo"], - "output_filename": "argo_float_uv.nc", - }, - "Sdata": { - "dataset_id": phys_product_id, - "variables": ["so"], - "output_filename": "argo_float_s.nc", - }, - "Tdata": { - "dataset_id": phys_product_id, - "variables": ["thetao"], - "output_filename": "argo_float_t.nc", - }, - } - - # Iterate over all datasets and download each based on space_time_region - try: - for dataset in argo_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("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...") - - bgc_args = { - "physical": False, - "schedule_start": start_datetime, - "schedule_end": end_datetime, - "username": username, - "password": password, - } - - ctd_bgc_download_dict = { - "o2data": { - "dataset_id": select_product_id(**{**bgc_args, "variable": "o2"}), - "variables": ["o2"], - "output_filename": "ctd_bgc_o2.nc", - }, - "chlorodata": { - "dataset_id": select_product_id(**{**bgc_args, "variable": "chl"}), - "variables": ["chl"], - "output_filename": "ctd_bgc_chl.nc", - }, - "nitratedata": { - "dataset_id": select_product_id(**{**bgc_args, "variable": "no3"}), - "variables": ["no3"], - "output_filename": "ctd_bgc_no3.nc", - }, - "phosphatedata": { - "dataset_id": select_product_id(**{**bgc_args, "variable": "po4"}), - "variables": ["po4"], - "output_filename": "ctd_bgc_po4.nc", - }, - "phdata": { - "dataset_id": select_product_id( - **{**bgc_args, "variable": "ph"} - ), # this will be monthly resolution if reanalysis(_interim) period - "variables": ["ph"], - "output_filename": "ctd_bgc_ph.nc", - }, - "phytoplanktondata": { - "dataset_id": select_product_id( - **{**bgc_args, "variable": "phyc"} - ), # this will be monthly resolution if reanalysis(_interim) period, - "variables": ["phyc"], - "output_filename": "ctd_bgc_phyc.nc", - }, - "primaryproductiondata": { - "dataset_id": select_product_id(**{**bgc_args, "variable": "nppv"}), - "variables": ["nppv"], - "output_filename": "ctd_bgc_nppv.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.""" - assert length % 2 == 0, "Length must be even." - half_length = length // 2 - - return hashlib.shake_128(s.encode("utf-8")).hexdigest(half_length) - - -def create_hash(s: str) -> str: - """Create an 8 digit hash of a string.""" - return _hash(s, length=8) - - -def hash_model(model: BaseModel, salt: int = 0) -> str: - """ - Hash a Pydantic model. - - :param region: The region to hash. - :param salt: Salt to add to the hash. - :returns: The hash. - """ - return create_hash(model.model_dump_json() + str(salt)) - - -def get_space_time_region_hash(space_time_region: SpaceTimeRegion) -> str: - # Increment salt in the event of breaking data fetching changes with prior versions - # of virtualship where you want to force new hashes (i.e., new data downloads) - salt = 0 - return hash_model(space_time_region, salt=salt) - - -def filename_to_hash(filename: str) -> str: - """Extract hash from filename of the format YYYYMMDD_HHMMSS_{hash}.""" - parts = filename.split("_") - if len(parts) != 3: - raise ValueError( - f"Filename '{filename}' must have 3 parts delimited with underscores." - ) - return parts[-1] - - -def hash_to_filename(hash: str) -> str: - """Return a filename of the format YYYYMMDD_HHMMSS_{hash}.""" - if "_" in hash: - raise ValueError("Hash cannot contain underscores.") - return f"{datetime.now().strftime('%Y%m%d_%H%M%S')}_{hash}" - - -class DownloadMetadata(BaseModel): - """Metadata for a data download.""" - - download_complete: bool - download_date: datetime | None = None - - def to_yaml(self, file_path: str | Path) -> None: - with open(file_path, "w") as file: - _dump_yaml(self, file) - - @classmethod - def from_yaml(cls, file_path: str | Path) -> DownloadMetadata: - return _generic_load_yaml(file_path, cls) - - -def get_existing_download( - data_folder: Path, space_time_region_hash: str -) -> Path | None: - """Check if a download has already been completed. If so, return the path for existing download.""" - for download_path in data_folder.rglob("*"): - try: - hash = filename_to_hash(download_path.name) - except ValueError: - continue - - if hash == space_time_region_hash: - assert_complete_download(download_path) - return download_path - - return None - - -def assert_complete_download(download_path: Path) -> None: - download_metadata = download_path / DOWNLOAD_METADATA - try: - with open(download_metadata) as file: - assert DownloadMetadata.from_yaml(file).download_complete - except (FileNotFoundError, AssertionError) as e: - raise IncompleteDownloadError( - f"Download at {download_path} was found, but looks to be incomplete " - f"(likely due to interupting it mid-download). Please delete this folder and retry." - ) from e - return - - -def complete_download(download_path: Path) -> None: - """Mark a download as complete.""" - download_metadata = download_path / DOWNLOAD_METADATA - metadata = DownloadMetadata(download_complete=True, download_date=datetime.now()) - metadata.to_yaml(download_metadata) - return - - -def select_product_id( - physical: bool, - schedule_start: datetime, - schedule_end: datetime, - username: str, - password: str, - variable: str | None = None, # only needed for BGC datasets -) -> str: - """ - Determine which copernicus product id should be selected (reanalysis, reanalysis-interim, analysis & forecast), for prescribed schedule and physical vs. BGC. - - BGC is more complicated than physical products. Often (re)analysis period and variable dependent, hence more custom logic here. - """ - product_ids = { - "phys": { - "reanalysis": "cmems_mod_glo_phy_my_0.083deg_P1D-m", - "reanalysis_interim": "cmems_mod_glo_phy_myint_0.083deg_P1D-m", - "analysis": "cmems_mod_glo_phy_anfc_0.083deg_P1D-m", - }, - "bgc": { - "reanalysis": "cmems_mod_glo_bgc_my_0.25deg_P1D-m", - "reanalysis_interim": "cmems_mod_glo_bgc_myint_0.25deg_P1D-m", - "analysis": None, # will be set per variable - }, - } - - bgc_analysis_ids = { - "o2": "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1D-m", - "chl": "cmems_mod_glo_bgc-pft_anfc_0.25deg_P1D-m", - "no3": "cmems_mod_glo_bgc-nut_anfc_0.25deg_P1D-m", - "po4": "cmems_mod_glo_bgc-nut_anfc_0.25deg_P1D-m", - "ph": "cmems_mod_glo_bgc-car_anfc_0.25deg_P1D-m", - "phyc": "cmems_mod_glo_bgc-pft_anfc_0.25deg_P1D-m", - "nppv": "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1D-m", - } - - # pH and phytoplankton variables are available as *monthly* products only in renalysis(_interim) period - monthly_bgc_reanalysis_ids = { - "ph": "cmems_mod_glo_bgc_my_0.25deg_P1M-m", - "phyc": "cmems_mod_glo_bgc_my_0.25deg_P1M-m", - } - monthly_bgc_reanalysis_interim_ids = { - "ph": "cmems_mod_glo_bgc_myint_0.25deg_P1M-m", - "phyc": "cmems_mod_glo_bgc_myint_0.25deg_P1M-m", - } - - key = "phys" if physical else "bgc" - selected_id = None - - for period, pid in product_ids[key].items(): - # for BGC analysis, set pid per variable - if key == "bgc" and period == "analysis": - if variable is None or variable not in bgc_analysis_ids: - continue - pid = bgc_analysis_ids[variable] - # for BGC reanalysis, check if requires monthly product - if ( - key == "bgc" - and period == "reanalysis" - and variable in monthly_bgc_reanalysis_ids - ): - monthly_pid = monthly_bgc_reanalysis_ids[variable] - ds_monthly = copernicusmarine.open_dataset( - monthly_pid, - username=username, - password=password, - ) - time_end_monthly = ds_monthly["time"][-1].values - if np.datetime64(schedule_end) <= time_end_monthly: - pid = monthly_pid - # for BGC reanalysis_interim, check if requires monthly product - if ( - key == "bgc" - and period == "reanalysis_interim" - and variable in monthly_bgc_reanalysis_interim_ids - ): - monthly_pid = monthly_bgc_reanalysis_interim_ids[variable] - ds_monthly = copernicusmarine.open_dataset( - monthly_pid, username=username, password=password - ) - time_end_monthly = ds_monthly["time"][-1].values - if np.datetime64(schedule_end) <= time_end_monthly: - pid = monthly_pid - if pid is None: - continue - ds = copernicusmarine.open_dataset(pid, username=username, password=password) - time_end = ds["time"][-1].values - if np.datetime64(schedule_end) <= time_end: - selected_id = pid - break - - if selected_id is None: - raise CopernicusCatalogueError( - "No suitable product found in the Copernicus Marine Catalogue for the scheduled time and variable." - ) - - # handle the rare situation where start time and end time span different products, which is possible for reanalysis and reanalysis_interim - # in this case, return the analysis product which spans far back enough - if start_end_in_product_timerange( - selected_id, schedule_start, schedule_end, username, password - ): - return selected_id - - else: - return ( - product_ids["phys"]["analysis"] if physical else bgc_analysis_ids[variable] - ) - - -def start_end_in_product_timerange( - selected_id: str, - schedule_start: datetime, - schedule_end: datetime, - username: str, - password: str, -) -> bool: - """Check schedule_start and schedule_end are both within a selected Copernicus product's time range.""" - ds_selected = copernicusmarine.open_dataset( - selected_id, username=username, password=password - ) - time_values = ds_selected["time"].values - time_min, time_max = np.min(time_values), np.max(time_values) - - if ( - np.datetime64(schedule_start) >= time_min - and np.datetime64(schedule_end) <= time_max - ): - return True - - else: - return False diff --git a/src/virtualship/cli/_plan.py b/src/virtualship/cli/_plan.py index eff3ac56..6aa6ff28 100644 --- a/src/virtualship/cli/_plan.py +++ b/src/virtualship/cli/_plan.py @@ -29,6 +29,7 @@ type_to_textual, ) from virtualship.errors import UnexpectedError, UserError +from virtualship.instruments.types import InstrumentType from virtualship.models import ( ADCPConfig, ArgoFloatConfig, @@ -36,7 +37,6 @@ CTDConfig, DrifterConfig, Expedition, - InstrumentType, Location, ShipConfig, ShipUnderwaterSTConfig, @@ -630,7 +630,7 @@ def _update_schedule(self): 0, ) wp.instrument = [] - for instrument in InstrumentType: + for instrument in [inst for inst in InstrumentType if not inst.is_underway]: switch_on = self.query_one(f"#wp{i}_{instrument.value}").value if instrument.value == "DRIFTER" and switch_on: count_str = self.query_one(f"#wp{i}_drifter_count").value @@ -847,7 +847,6 @@ def compose(self) -> ComposeResult: yield Select( [ (str(year), year) - # TODO: change from hard coding? ...flexibility for different datasets... for year in range( 1993, datetime.datetime.now().year + 1, @@ -902,7 +901,7 @@ def compose(self) -> ComposeResult: ) yield Label("Instruments:") - for instrument in InstrumentType: + for instrument in [i for i in InstrumentType if not i.is_underway]: is_selected = instrument in (self.waypoint.instrument or []) with Horizontal(): yield Label(instrument.value) @@ -952,7 +951,9 @@ def copy_from_previous(self) -> None: if prev and curr: curr.value = prev.value - for instrument in InstrumentType: + for instrument in [ + inst for inst in InstrumentType if not inst.is_underway + ]: prev_switch = schedule_editor.query_one( f"#wp{self.index - 1}_{instrument.value}" ) @@ -1045,9 +1046,9 @@ def save_pressed(self) -> None: # verify schedule expedition_editor.expedition.schedule.verify( ship_speed_value, - input_data=None, + bathy_data_dir=None, check_space_time_region=True, - ignore_missing_fieldsets=True, + ignore_missing_bathymetry=True, ) expedition_saved = expedition_editor.save_changes() diff --git a/src/virtualship/expedition/do_expedition.py b/src/virtualship/cli/_run.py similarity index 59% rename from src/virtualship/expedition/do_expedition.py rename to src/virtualship/cli/_run.py index 5c46d2eb..df8a24fd 100644 --- a/src/virtualship/expedition/do_expedition.py +++ b/src/virtualship/cli/_run.py @@ -1,35 +1,51 @@ """do_expedition function.""" +import logging import os import shutil from pathlib import Path import pyproj -from virtualship.cli._fetch import get_existing_download, get_space_time_region_hash -from virtualship.models import Expedition, Schedule +from virtualship.expedition.simulate_schedule import ( + MeasurementsToSimulate, + ScheduleProblem, + simulate_schedule, +) +from virtualship.models import Schedule +from virtualship.models.checkpoint import Checkpoint from virtualship.utils import ( CHECKPOINT, _get_expedition, + expedition_cost, + get_instrument_class, ) -from .checkpoint import Checkpoint -from .expedition_cost import expedition_cost -from .input_data import InputData -from .simulate_measurements import simulate_measurements -from .simulate_schedule import ScheduleProblem, simulate_schedule - # projection used to sail between waypoints projection = pyproj.Geod(ellps="WGS84") -def do_expedition(expedition_dir: str | Path, input_data: Path | None = None) -> None: +# parcels logger (suppress INFO messages to prevent log being flooded) +external_logger = logging.getLogger("parcels.tools.loggers") +external_logger.setLevel(logging.WARNING) + +# copernicusmarine logger (suppress INFO messages to prevent log being flooded) +logging.getLogger("copernicusmarine").setLevel("ERROR") + + +def _run(expedition_dir: str | Path) -> None: """ Perform an expedition, providing terminal feedback and file output. :param expedition_dir: The base directory for the expedition. - :param input_data: Input data folder (override used for testing). """ + # ################################# TEMPORARY TIMER: START ################################# + import time + + start_time = time.time() + print("[TIMER] Expedition started...") + # ################################# TEMPORARY TIMER: START ################################# + print("\n╔═════════════════════════════════════════════════╗") print("║ VIRTUALSHIP EXPEDITION STATUS ║") print("╚═════════════════════════════════════════════════╝") @@ -40,7 +56,7 @@ def do_expedition(expedition_dir: str | Path, input_data: Path | None = None) -> expedition = _get_expedition(expedition_dir) # Verify instruments_config file is consistent with schedule - expedition.instruments_config.verify(expedition.schedule) + expedition.instruments_config.verify(expedition) # load last checkpoint checkpoint = _load_checkpoint(expedition_dir) @@ -50,19 +66,9 @@ def do_expedition(expedition_dir: str | Path, input_data: Path | None = None) -> # verify that schedule and checkpoint match checkpoint.verify(expedition.schedule) - # load fieldsets - loaded_input_data = _load_input_data( - expedition_dir=expedition_dir, - expedition=expedition, - input_data=input_data, - ) - print("\n---- WAYPOINT VERIFICATION ----") - # verify schedule is valid - expedition.schedule.verify( - expedition.ship_config.ship_speed_knots, loaded_input_data - ) + expedition.schedule.verify(expedition.ship_config.ship_speed_knots) # simulate the schedule schedule_results = simulate_schedule( @@ -92,26 +98,38 @@ def do_expedition(expedition_dir: str | Path, input_data: Path | None = None) -> print("\n----- EXPEDITION SUMMARY ------") - # calculate expedition cost in US$ - assert expedition.schedule.waypoints[0].time is not None, ( - "First waypoint has no time. This should not be possible as it should have been verified before." - ) - time_past = schedule_results.time - expedition.schedule.waypoints[0].time - cost = expedition_cost(schedule_results, time_past) - with open(expedition_dir.joinpath("results", "cost.txt"), "w") as file: - file.writelines(f"cost: {cost} US$") - print(f"\nExpedition duration: {time_past}\nExpedition cost: US$ {cost:,.0f}.") + # expedition cost in US$ + _write_expedition_cost(expedition, schedule_results, expedition_dir) print("\n--- MEASUREMENT SIMULATIONS ---") # simulate measurements print("\nSimulating measurements. This may take a while...\n") - simulate_measurements( - expedition_dir, - expedition.instruments_config, - loaded_input_data, - schedule_results.measurements_to_simulate, - ) + + instruments_in_expedition = expedition.get_instruments() + + for itype in instruments_in_expedition: + # get instrument class + instrument_class = get_instrument_class(itype) + if instrument_class is None: + raise RuntimeError(f"No instrument class found for type {itype}.") + + # get measurements to simulate + attr = MeasurementsToSimulate.get_attr_for_instrumenttype(itype) + measurements = getattr(schedule_results.measurements_to_simulate, attr) + + # initialise instrument + instrument = instrument_class( + expedition=expedition, + directory=expedition_dir, + ) + + # execute simulation + instrument.execute( + measurements=measurements, + out_path=expedition_dir.joinpath("results", f"{itype.name.lower()}.zarr"), + ) + print("\nAll measurement simulations are complete.") print("\n----- EXPEDITION RESULTS ------") @@ -121,41 +139,11 @@ def do_expedition(expedition_dir: str | Path, input_data: Path | None = None) -> ) print("\n------------- END -------------\n") - -def _load_input_data( - expedition_dir: Path, - expedition: Expedition, - input_data: Path | None, -) -> InputData: - """ - Load the input data. - - :param expedition_dir: Directory of the expedition. - :param expedition: Expedition object. - :param input_data: Folder containing input data. - :return: InputData object. - """ - if input_data is None: - space_time_region_hash = get_space_time_region_hash( - expedition.schedule.space_time_region - ) - input_data = get_existing_download(expedition_dir, space_time_region_hash) - - assert input_data is not None, ( - "Input data hasn't been found. Have you run the `virtualship fetch` command?" - ) - - return InputData.load( - directory=input_data, - load_adcp=expedition.instruments_config.adcp_config is not None, - load_argo_float=expedition.instruments_config.argo_float_config is not None, - load_ctd=expedition.instruments_config.ctd_config is not None, - load_ctd_bgc=expedition.instruments_config.ctd_bgc_config is not None, - load_drifter=expedition.instruments_config.drifter_config is not None, - load_xbt=expedition.instruments_config.xbt_config is not None, - load_ship_underwater_st=expedition.instruments_config.ship_underwater_st_config - is not None, - ) + ################################# TEMPORARY TIMER: END ################################# + end_time = time.time() + elapsed = end_time - start_time + print(f"[TIMER] Expedition completed in {elapsed:.2f} seconds.") + ################################# TEMPORARY TIMER: END ################################# def _load_checkpoint(expedition_dir: Path) -> Checkpoint | None: @@ -169,3 +157,15 @@ def _load_checkpoint(expedition_dir: Path) -> Checkpoint | None: def _save_checkpoint(checkpoint: Checkpoint, expedition_dir: Path) -> None: file_path = expedition_dir.joinpath(CHECKPOINT) checkpoint.to_yaml(file_path) + + +def _write_expedition_cost(expedition, schedule_results, expedition_dir): + """Calculate the expedition cost, write it to a file, and print summary.""" + assert expedition.schedule.waypoints[0].time is not None, ( + "First waypoint has no time. This should not be possible as it should have been verified before." + ) + time_past = schedule_results.time - expedition.schedule.waypoints[0].time + cost = expedition_cost(schedule_results, time_past) + with open(expedition_dir.joinpath("results", "cost.txt"), "w") as file: + file.writelines(f"cost: {cost} US$") + print(f"\nExpedition duration: {time_past}\nExpedition cost: US$ {cost:,.0f}.") diff --git a/src/virtualship/cli/commands.py b/src/virtualship/cli/commands.py index 3e83be3b..2b0ce8e9 100644 --- a/src/virtualship/cli/commands.py +++ b/src/virtualship/cli/commands.py @@ -3,9 +3,8 @@ import click from virtualship import utils -from virtualship.cli._fetch import _fetch from virtualship.cli._plan import _plan -from virtualship.expedition.do_expedition import do_expedition +from virtualship.cli._run import _run from virtualship.utils import ( EXPEDITION, mfp_to_yaml, @@ -69,48 +68,19 @@ def init(path, from_mfp): ) def plan(path): """ - Launch UI to help build schedule and ship config files. + Launch UI to help build expedition configuration (YAML) file. Should you encounter any issues with using this tool, please report an issue describing the problem to the VirtualShip issue tracker at: https://github.com/OceanParcels/virtualship/issues" """ _plan(Path(path)) -@click.command() -@click.argument( - "path", - type=click.Path(exists=True, file_okay=False, dir_okay=True, readable=True), -) -@click.option( - "--username", - type=str, - default=None, - help="Copernicus Marine username.", -) -@click.option( - "--password", - type=str, - default=None, - help="Copernicus Marine password.", -) -def fetch(path: str | Path, username: str | None, password: str | None) -> None: - """ - Download input data for an expedition. - - Entrypoint for the tool to download data based on space-time region provided in the - schedule file. Data is downloaded from Copernicus Marine, credentials for which can be - obtained via registration: https://data.marine.copernicus.eu/register . Credentials can - be provided on prompt, via command line arguments, or via a YAML config file. Run - `virtualship fetch` on a expedition for more info. - """ - _fetch(path, username, password) - - +# TODO: also add option to 'stream' via link to dir elsewhere, e.g. simlink or path to data stored elsewhere that isn't expedition dir! @click.command() @click.argument( "path", type=click.Path(exists=True, file_okay=False, dir_okay=True, readable=True), ) def run(path): - """Run the expedition.""" - do_expedition(Path(path)) + """Execute the expedition simulations.""" + _run(Path(path)) diff --git a/src/virtualship/cli/main.py b/src/virtualship/cli/main.py index 6ee12eff..a02a5ffb 100644 --- a/src/virtualship/cli/main.py +++ b/src/virtualship/cli/main.py @@ -11,7 +11,6 @@ def cli(): cli.add_command(commands.init) cli.add_command(commands.plan) -cli.add_command(commands.fetch) cli.add_command(commands.run) if __name__ == "__main__": diff --git a/src/virtualship/errors.py b/src/virtualship/errors.py index 6026936e..fad5863d 100644 --- a/src/virtualship/errors.py +++ b/src/virtualship/errors.py @@ -22,7 +22,7 @@ class ScheduleError(RuntimeError): pass -class ConfigError(RuntimeError): +class InstrumentsConfigError(RuntimeError): """An error in the config.""" pass @@ -40,6 +40,12 @@ class UnexpectedError(Exception): pass +class UnderwayConfigsError(Exception): + """Error raised when underway instrument configurations (ADCP or underwater ST) are missing.""" + + pass + + class CopernicusCatalogueError(Exception): """Error raised when a relevant product is not found in the Copernicus Catalogue.""" diff --git a/src/virtualship/expedition/__init__.py b/src/virtualship/expedition/__init__.py index 43d24844..7f072bbf 100644 --- a/src/virtualship/expedition/__init__.py +++ b/src/virtualship/expedition/__init__.py @@ -1,9 +1 @@ -"""Everything for simulating an expedition.""" - -from .do_expedition import do_expedition -from .input_data import InputData - -__all__ = [ - "InputData", - "do_expedition", -] +"""Simulating an expedition.""" diff --git a/src/virtualship/expedition/expedition_cost.py b/src/virtualship/expedition/expedition_cost.py deleted file mode 100644 index cab6ab7d..00000000 --- a/src/virtualship/expedition/expedition_cost.py +++ /dev/null @@ -1,27 +0,0 @@ -"""expedition_cost function.""" - -from datetime import timedelta - -from .simulate_schedule import ScheduleOk - - -def expedition_cost(schedule_results: ScheduleOk, time_past: timedelta) -> float: - """ - Calculate the cost of the expedition in US$. - - :param schedule_results: Results from schedule simulation. - :param time_past: Time the expedition took. - :returns: The calculated cost of the expedition in US$. - """ - SHIP_COST_PER_DAY = 30000 - DRIFTER_DEPLOY_COST = 2500 - ARGO_DEPLOY_COST = 15000 - - ship_cost = SHIP_COST_PER_DAY / 24 * time_past.total_seconds() // 3600 - num_argos = len(schedule_results.measurements_to_simulate.argo_floats) - argo_cost = num_argos * ARGO_DEPLOY_COST - num_drifters = len(schedule_results.measurements_to_simulate.drifters) - drifter_cost = num_drifters * DRIFTER_DEPLOY_COST - - cost = ship_cost + argo_cost + drifter_cost - return cost diff --git a/src/virtualship/expedition/input_data.py b/src/virtualship/expedition/input_data.py deleted file mode 100644 index fa48e0a7..00000000 --- a/src/virtualship/expedition/input_data.py +++ /dev/null @@ -1,255 +0,0 @@ -"""InputData class.""" - -from __future__ import annotations - -from dataclasses import dataclass -from pathlib import Path - -from parcels import Field, FieldSet - - -@dataclass -class InputData: - """A collection of fieldsets that function as input data for simulation.""" - - 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 - - @classmethod - def load( - cls, - directory: str | Path, - 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, - ) -> InputData: - """ - Create an instance of this class from netCDF files. - - For now this function makes a lot of assumption about file location and contents. - - :param directory: Input data directory. - :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. - """ - directory = Path(directory) - if load_drifter: - drifter_fieldset = cls._load_drifter_fieldset(directory) - else: - drifter_fieldset = None - if load_argo_float: - 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: - adcp_fieldset = ship_fieldset - else: - adcp_fieldset = None - if load_ctd: - ctd_fieldset = ship_fieldset - else: - ctd_fieldset = None - if load_ship_underwater_st: - ship_underwater_st_fieldset = ship_fieldset - else: - ship_underwater_st_fieldset = None - if load_xbt: - xbt_fieldset = ship_fieldset - else: - xbt_fieldset = None - - return InputData( - 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, - ) - - @classmethod - def _load_ship_fieldset(cls, directory: Path) -> FieldSet: - filenames = { - "U": directory.joinpath("ship_uv.nc"), - "V": directory.joinpath("ship_uv.nc"), - "S": directory.joinpath("ship_s.nc"), - "T": directory.joinpath("ship_t.nc"), - } - variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} - dimensions = { - "lon": "longitude", - "lat": "latitude", - "time": "time", - "depth": "depth", - } - - # create the fieldset and set interpolation methods - fieldset = FieldSet.from_netcdf( - filenames, variables, dimensions, allow_time_extrapolation=True - ) - fieldset.T.interp_method = "linear_invdist_land_tracer" - fieldset.S.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_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_chl.nc"), - "no3": directory.joinpath("ctd_bgc_no3.nc"), - "po4": directory.joinpath("ctd_bgc_po4.nc"), - "ph": directory.joinpath("ctd_bgc_ph.nc"), - "phyc": directory.joinpath("ctd_bgc_phyc.nc"), - "nppv": directory.joinpath("ctd_bgc_nppv.nc"), - } - variables = { - "U": "uo", - "V": "vo", - "o2": "o2", - "chl": "chl", - "no3": "no3", - "po4": "po4", - "ph": "ph", - "phyc": "phyc", - "nppv": "nppv", - } - 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" - fieldset.no3.interp_method = "linear_invdist_land_tracer" - fieldset.po4.interp_method = "linear_invdist_land_tracer" - fieldset.ph.interp_method = "linear_invdist_land_tracer" - fieldset.phyc.interp_method = "linear_invdist_land_tracer" - fieldset.nppv.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 = { - "U": directory.joinpath("drifter_uv.nc"), - "V": directory.joinpath("drifter_uv.nc"), - "T": directory.joinpath("drifter_t.nc"), - } - variables = {"U": "uo", "V": "vo", "T": "thetao"} - dimensions = { - "lon": "longitude", - "lat": "latitude", - "time": "time", - "depth": "depth", - } - - fieldset = FieldSet.from_netcdf( - filenames, variables, dimensions, allow_time_extrapolation=False - ) - fieldset.T.interp_method = "linear_invdist_land_tracer" - - # make depth negative - for g in fieldset.gridset.grids: - g.negate_depth() - - # read in data already - fieldset.computeTimeChunk(0, 1) - - return fieldset - - @classmethod - def _load_argo_float_fieldset(cls, directory: Path) -> FieldSet: - filenames = { - "U": directory.joinpath("argo_float_uv.nc"), - "V": directory.joinpath("argo_float_uv.nc"), - "S": directory.joinpath("argo_float_s.nc"), - "T": directory.joinpath("argo_float_t.nc"), - } - variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} - dimensions = { - "lon": "longitude", - "lat": "latitude", - "time": "time", - "depth": "depth", - } - - fieldset = FieldSet.from_netcdf( - filenames, variables, dimensions, allow_time_extrapolation=False - ) - fieldset.T.interp_method = "linear_invdist_land_tracer" - fieldset.S.interp_method = "linear_invdist_land_tracer" - - # make depth negative - for g in fieldset.gridset.grids: - if max(g.depth) > 0: - g.negate_depth() - - # read in data already - fieldset.computeTimeChunk(0, 1) - - return fieldset diff --git a/src/virtualship/expedition/simulate_measurements.py b/src/virtualship/expedition/simulate_measurements.py deleted file mode 100644 index 6cb2e488..00000000 --- a/src/virtualship/expedition/simulate_measurements.py +++ /dev/null @@ -1,162 +0,0 @@ -"""simulate_measurements function.""" - -from __future__ import annotations - -import logging -from datetime import timedelta -from pathlib import Path -from typing import TYPE_CHECKING - -from yaspin import yaspin - -from virtualship.instruments.adcp import simulate_adcp -from virtualship.instruments.argo_float import simulate_argo_floats -from virtualship.instruments.ctd import simulate_ctd -from virtualship.instruments.ctd_bgc import simulate_ctd_bgc -from virtualship.instruments.drifter import simulate_drifters -from virtualship.instruments.ship_underwater_st import simulate_ship_underwater_st -from virtualship.instruments.xbt import simulate_xbt -from virtualship.models import InstrumentsConfig -from virtualship.utils import ship_spinner - -from .simulate_schedule import MeasurementsToSimulate - -if TYPE_CHECKING: - from .input_data import InputData - -# parcels logger (suppress INFO messages to prevent log being flooded) -external_logger = logging.getLogger("parcels.tools.loggers") -external_logger.setLevel(logging.WARNING) - - -def simulate_measurements( - expedition_dir: str | Path, - instruments_config: InstrumentsConfig, - input_data: InputData, - measurements: MeasurementsToSimulate, -) -> None: - """ - Simulate measurements using Parcels. - - Saves everything in expedition_dir/results. - - :param expedition_dir: Base directory of the expedition. - :param input_data: Input data for simulation. - :param measurements: The measurements to simulate. - :raises RuntimeError: In case fieldsets of configuration is not provided. Make sure to check this before calling this function. - """ - if isinstance(expedition_dir, str): - expedition_dir = Path(expedition_dir) - - if len(measurements.ship_underwater_sts) > 0: - if instruments_config.ship_underwater_st_config is None: - raise RuntimeError("No configuration for ship underwater ST provided.") - if input_data.ship_underwater_st_fieldset is None: - raise RuntimeError("No fieldset for ship underwater ST provided.") - with yaspin( - text="Simulating onboard temperature and salinity measurements... ", - side="right", - spinner=ship_spinner, - ) as spinner: - simulate_ship_underwater_st( - fieldset=input_data.ship_underwater_st_fieldset, - out_path=expedition_dir.joinpath("results", "ship_underwater_st.zarr"), - depth=-2, - sample_points=measurements.ship_underwater_sts, - ) - spinner.ok("✅") - - if len(measurements.adcps) > 0: - if instruments_config.adcp_config is None: - raise RuntimeError("No configuration for ADCP provided.") - if input_data.adcp_fieldset is None: - raise RuntimeError("No fieldset for ADCP provided.") - with yaspin( - text="Simulating onboard ADCP... ", side="right", spinner=ship_spinner - ) as spinner: - simulate_adcp( - fieldset=input_data.adcp_fieldset, - out_path=expedition_dir.joinpath("results", "adcp.zarr"), - max_depth=instruments_config.adcp_config.max_depth_meter, - min_depth=-5, - num_bins=instruments_config.adcp_config.num_bins, - sample_points=measurements.adcps, - ) - spinner.ok("✅") - - if len(measurements.ctds) > 0: - if instruments_config.ctd_config is None: - raise RuntimeError("No configuration for CTD provided.") - if input_data.ctd_fieldset is None: - raise RuntimeError("No fieldset for CTD provided.") - with yaspin( - text="Simulating CTD casts... ", side="right", spinner=ship_spinner - ) as spinner: - simulate_ctd( - out_path=expedition_dir.joinpath("results", "ctd.zarr"), - fieldset=input_data.ctd_fieldset, - ctds=measurements.ctds, - outputdt=timedelta(seconds=10), - ) - spinner.ok("✅") - - if len(measurements.ctd_bgcs) > 0: - if instruments_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.") - with yaspin( - text="Simulating BGC CTD casts... ", side="right", spinner=ship_spinner - ) as spinner: - 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), - ) - spinner.ok("✅") - - if len(measurements.xbts) > 0: - if instruments_config.xbt_config is None: - raise RuntimeError("No configuration for XBTs provided.") - if input_data.xbt_fieldset is None: - raise RuntimeError("No fieldset for XBTs provided.") - with yaspin( - text="Simulating XBTs... ", side="right", spinner=ship_spinner - ) as spinner: - simulate_xbt( - out_path=expedition_dir.joinpath("results", "xbts.zarr"), - fieldset=input_data.xbt_fieldset, - xbts=measurements.xbts, - outputdt=timedelta(seconds=1), - ) - spinner.ok("✅") - - if len(measurements.drifters) > 0: - print("Simulating drifters... ") - if instruments_config.drifter_config is None: - raise RuntimeError("No configuration for drifters provided.") - if input_data.drifter_fieldset is None: - raise RuntimeError("No fieldset for drifters provided.") - simulate_drifters( - out_path=expedition_dir.joinpath("results", "drifters.zarr"), - fieldset=input_data.drifter_fieldset, - drifters=measurements.drifters, - outputdt=timedelta(hours=5), - dt=timedelta(minutes=5), - endtime=None, - ) - - if len(measurements.argo_floats) > 0: - print("Simulating argo floats... ") - if instruments_config.argo_float_config is None: - raise RuntimeError("No configuration for argo floats provided.") - if input_data.argo_float_fieldset is None: - raise RuntimeError("No fieldset for argo floats provided.") - simulate_argo_floats( - out_path=expedition_dir.joinpath("results", "argo_floats.zarr"), - argo_floats=measurements.argo_floats, - fieldset=input_data.argo_float_fieldset, - outputdt=timedelta(minutes=5), - endtime=None, - ) diff --git a/src/virtualship/expedition/simulate_schedule.py b/src/virtualship/expedition/simulate_schedule.py index 3b78c5c7..3b8fa78a 100644 --- a/src/virtualship/expedition/simulate_schedule.py +++ b/src/virtualship/expedition/simulate_schedule.py @@ -4,6 +4,7 @@ from dataclasses import dataclass, field from datetime import datetime, timedelta +from typing import ClassVar import pyproj @@ -11,10 +12,10 @@ from virtualship.instruments.ctd import CTD from virtualship.instruments.ctd_bgc import CTD_BGC from virtualship.instruments.drifter import Drifter +from virtualship.instruments.types import InstrumentType from virtualship.instruments.xbt import XBT from virtualship.models import ( Expedition, - InstrumentType, Location, Spacetime, Waypoint, @@ -39,7 +40,26 @@ class ScheduleProblem: @dataclass class MeasurementsToSimulate: - """The measurements to simulate, as concluded from schedule simulation.""" + """ + The measurements to simulate, as concluded from schedule simulation. + + Provides a mapping from InstrumentType to the correct attribute name for robust access. + """ + + _instrumenttype_to_attr: ClassVar[dict] = { + InstrumentType.ADCP: "adcps", + InstrumentType.UNDERWATER_ST: "ship_underwater_sts", + InstrumentType.ARGO_FLOAT: "argo_floats", + InstrumentType.DRIFTER: "drifters", + InstrumentType.CTD: "ctds", + InstrumentType.CTD_BGC: "ctd_bgcs", + InstrumentType.XBT: "xbts", + } + + @classmethod + def get_attr_for_instrumenttype(cls, instrument_type): + """Return the attribute name for a given InstrumentType.""" + return cls._instrumenttype_to_attr[instrument_type] adcps: list[Spacetime] = field(default_factory=list, init=False) ship_underwater_sts: list[Spacetime] = field(default_factory=list, init=False) diff --git a/src/virtualship/instruments/__init__.py b/src/virtualship/instruments/__init__.py index 6a6ffbca..b593ed38 100644 --- a/src/virtualship/instruments/__init__.py +++ b/src/virtualship/instruments/__init__.py @@ -1,6 +1,14 @@ -"""Measurement instrument that can be used with Parcels.""" +"""Instruments in VirtualShip.""" -from . import adcp, argo_float, ctd, ctd_bgc, drifter, ship_underwater_st, xbt +from . import ( + adcp, + argo_float, + ctd, + ctd_bgc, + drifter, + ship_underwater_st, + xbt, +) __all__ = [ "adcp", diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py index af2c285e..e524c5a2 100644 --- a/src/virtualship/instruments/adcp.py +++ b/src/virtualship/instruments/adcp.py @@ -1,14 +1,32 @@ -"""ADCP instrument.""" - -from pathlib import Path +from dataclasses import dataclass +from typing import ClassVar import numpy as np -from parcels import FieldSet, ParticleSet, ScipyParticle, Variable +from parcels import ParticleSet, ScipyParticle, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.utils import ( + register_instrument, +) + +# ===================================================== +# SECTION: Dataclass +# ===================================================== + + +@dataclass +class ADCP: + """ADCP configuration.""" + + name: ClassVar[str] = "ADCP" + + +# ===================================================== +# SECTION: Particle Class +# ===================================================== -from virtualship.models import Spacetime -# we specifically use ScipyParticle because we have many small calls to execute -# there is some overhead with JITParticle and this ends up being significantly faster _ADCPParticle = ScipyParticle.add_variables( [ Variable("U", dtype=np.float32, initial=np.nan), @@ -16,6 +34,10 @@ ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + def _sample_velocity(particle, fieldset, time): particle.U, particle.V = fieldset.UV.eval( @@ -23,56 +45,72 @@ def _sample_velocity(particle, fieldset, time): ) -def simulate_adcp( - fieldset: FieldSet, - out_path: str | Path, - max_depth: float, - min_depth: float, - num_bins: int, - sample_points: list[Spacetime], -) -> None: - """ - Use Parcels to simulate an ADCP in a fieldset. - - :param fieldset: The fieldset to simulate the ADCP in. - :param out_path: The path to write the results to. - :param max_depth: Maximum depth the ADCP can measure. - :param min_depth: Minimum depth the ADCP can measure. - :param num_bins: How many samples to take in the complete range between max_depth and min_depth. - :param sample_points: The places and times to sample at. - """ - sample_points.sort(key=lambda p: p.time) - - bins = np.linspace(max_depth, min_depth, num_bins) - num_particles = len(bins) - particleset = ParticleSet.from_list( - fieldset=fieldset, - pclass=_ADCPParticle, - lon=np.full( - num_particles, 0.0 - ), # initial lat/lon are irrelevant and will be overruled later. - lat=np.full(num_particles, 0.0), - depth=bins, - time=0, # same for time - ) +# ===================================================== +# SECTION: Instrument Class +# ===================================================== - # define output file for the simulation - # outputdt set to infinite as we just want to write at the end of every call to 'execute' - out_file = particleset.ParticleFile(name=out_path, outputdt=np.inf) - for point in sample_points: - particleset.lon_nextloop[:] = point.location.lon - particleset.lat_nextloop[:] = point.location.lat - particleset.time_nextloop[:] = fieldset.time_origin.reltime( - np.datetime64(point.time) - ) +@register_instrument(InstrumentType.ADCP) +class ADCPInstrument(Instrument): + """ADCP instrument class.""" + + def __init__(self, expedition, directory): + """Initialize ADCPInstrument.""" + filenames = { + "U": f"{ADCP.name}_uv.nc", + "V": f"{ADCP.name}_uv.nc", + } + variables = {"U": "uo", "V": "vo"} - # perform one step using the particleset - # dt and runtime are set so exactly one step is made. - particleset.execute( - [_sample_velocity], - dt=1, - runtime=1, + super().__init__( + ADCP.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=False, + allow_time_extrapolation=True, verbose_progress=False, - output_file=out_file, + buffer_spec=None, + limit_spec=None, + ) + + def simulate(self, measurements, out_path) -> None: + """Simulate ADCP measurements.""" + MAX_DEPTH = self.expedition.instruments_config.adcp_config.max_depth_meter + MIN_DEPTH = -5.0 + NUM_BINS = self.expedition.instruments_config.adcp_config.num_bins + + measurements.sort(key=lambda p: p.time) + + fieldset = self.load_input_data() + + bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS) + num_particles = len(bins) + particleset = ParticleSet.from_list( + fieldset=fieldset, + pclass=_ADCPParticle, + lon=np.full( + num_particles, 0.0 + ), # initial lat/lon are irrelevant and will be overruled later.s + lat=np.full(num_particles, 0.0), + depth=bins, + time=0, ) + + out_file = particleset.ParticleFile(name=out_path, outputdt=np.inf) + + for point in measurements: + particleset.lon_nextloop[:] = point.location.lon + particleset.lat_nextloop[:] = point.location.lat + particleset.time_nextloop[:] = fieldset.time_origin.reltime( + np.datetime64(point.time) + ) + + particleset.execute( + [_sample_velocity], + dt=1, + runtime=1, + verbose_progress=self.verbose_progress, + output_file=out_file, + ) diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py index d0976367..f064e81a 100644 --- a/src/virtualship/instruments/argo_float.py +++ b/src/virtualship/instruments/argo_float.py @@ -1,27 +1,32 @@ -"""Argo float instrument.""" - import math from dataclasses import dataclass -from datetime import datetime, timedelta -from pathlib import Path +from datetime import timedelta +from typing import ClassVar import numpy as np from parcels import ( AdvectionRK4, - FieldSet, JITParticle, ParticleSet, StatusCode, Variable, ) -from virtualship.models import Spacetime +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.models.spacetime import Spacetime +from virtualship.utils import register_instrument + +# ===================================================== +# SECTION: Dataclass +# ===================================================== @dataclass class ArgoFloat: - """Configuration for a single Argo float.""" + """Argo float configuration.""" + name: ClassVar[str] = "ArgoFloat" spacetime: Spacetime min_depth: float max_depth: float @@ -31,6 +36,10 @@ class ArgoFloat: drift_days: float +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _ArgoParticle = JITParticle.add_variables( [ Variable("cycle_phase", dtype=np.int32, initial=0.0), @@ -47,6 +56,10 @@ class ArgoFloat: ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + def _argo_float_vertical_movement(particle, fieldset, time): if particle.cycle_phase == 0: @@ -115,72 +128,100 @@ def _check_error(particle, fieldset, time): particle.delete() -def simulate_argo_floats( - fieldset: FieldSet, - out_path: str | Path, - argo_floats: list[ArgoFloat], - outputdt: timedelta, - endtime: datetime | None, -) -> None: - """ - Use Parcels to simulate a set of Argo floats in a fieldset. - - :param fieldset: The fieldset to simulate the Argo floats in. - :param out_path: The path to write the results to. - :param argo_floats: A list of Argo floats to simulate. - :param outputdt: Interval which dictates the update frequency of file output during simulation - :param endtime: Stop at this time, or if None, continue until the end of the fieldset. - """ - DT = 10.0 # dt of Argo float simulation integrator - - if len(argo_floats) == 0: - print( - "No Argo floats provided. Parcels currently crashes when providing an empty particle set, so no argo floats simulation will be done and no files will be created." +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.ARGO_FLOAT) +class ArgoFloatInstrument(Instrument): + """ArgoFloat instrument class.""" + + def __init__(self, expedition, directory): + """Initialize ArgoFloatInstrument.""" + filenames = { + "U": f"{ArgoFloat.name}_uv.nc", + "V": f"{ArgoFloat.name}_uv.nc", + "S": f"{ArgoFloat.name}_s.nc", + "T": f"{ArgoFloat.name}_t.nc", + } + variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} + buffer_spec = { + "latlon": 6.0, # [degrees] + "time": 21.0, # [days] + } + + super().__init__( + ArgoFloat.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=False, + allow_time_extrapolation=False, + verbose_progress=True, + buffer_spec=buffer_spec, + limit_spec=None, + ) + + def simulate(self, measurements, out_path) -> None: + """Simulate Argo float measurements.""" + DT = 10.0 # dt of Argo float simulation integrator + OUTPUT_DT = timedelta(minutes=5) + ENDTIME = None + + if len(measurements) == 0: + print( + "No Argo floats provided. Parcels currently crashes when providing an empty particle set, so no argo floats simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + + fieldset = self.load_input_data() + + # define parcel particles + argo_float_particleset = ParticleSet( + fieldset=fieldset, + pclass=_ArgoParticle, + lat=[argo.spacetime.location.lat for argo in measurements], + lon=[argo.spacetime.location.lon for argo in measurements], + depth=[argo.min_depth for argo in measurements], + time=[argo.spacetime.time for argo in measurements], + min_depth=[argo.min_depth for argo in measurements], + max_depth=[argo.max_depth for argo in measurements], + drift_depth=[argo.drift_depth for argo in measurements], + vertical_speed=[argo.vertical_speed for argo in measurements], + cycle_days=[argo.cycle_days for argo in measurements], + drift_days=[argo.drift_days for argo in measurements], + ) + + # define output file for the simulation + out_file = argo_float_particleset.ParticleFile( + name=out_path, + outputdt=OUTPUT_DT, + chunks=[len(argo_float_particleset), 100], + ) + + # get earliest between fieldset end time and provide end time + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + if ENDTIME is None: + actual_endtime = fieldset_endtime + elif ENDTIME > fieldset_endtime: + print("WARN: Requested end time later than fieldset end time.") + actual_endtime = fieldset_endtime + else: + actual_endtime = np.timedelta64(ENDTIME) + + # execute simulation + argo_float_particleset.execute( + [ + _argo_float_vertical_movement, + AdvectionRK4, + _keep_at_surface, + _check_error, + ], + endtime=actual_endtime, + dt=DT, + output_file=out_file, + verbose_progress=self.verbose_progress, ) - # TODO when Parcels supports it this check can be removed. - return - - # define parcel particles - argo_float_particleset = ParticleSet( - fieldset=fieldset, - pclass=_ArgoParticle, - lat=[argo.spacetime.location.lat for argo in argo_floats], - lon=[argo.spacetime.location.lon for argo in argo_floats], - depth=[argo.min_depth for argo in argo_floats], - time=[argo.spacetime.time for argo in argo_floats], - min_depth=[argo.min_depth for argo in argo_floats], - max_depth=[argo.max_depth for argo in argo_floats], - drift_depth=[argo.drift_depth for argo in argo_floats], - vertical_speed=[argo.vertical_speed for argo in argo_floats], - cycle_days=[argo.cycle_days for argo in argo_floats], - drift_days=[argo.drift_days for argo in argo_floats], - ) - - # define output file for the simulation - out_file = argo_float_particleset.ParticleFile( - name=out_path, outputdt=outputdt, chunks=[len(argo_float_particleset), 100] - ) - - # get earliest between fieldset end time and provide end time - fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) - if endtime is None: - actual_endtime = fieldset_endtime - elif endtime > fieldset_endtime: - print("WARN: Requested end time later than fieldset end time.") - actual_endtime = fieldset_endtime - else: - actual_endtime = np.timedelta64(endtime) - - # execute simulation - argo_float_particleset.execute( - [ - _argo_float_vertical_movement, - AdvectionRK4, - _keep_at_surface, - _check_error, - ], - endtime=actual_endtime, - dt=DT, - output_file=out_file, - verbose_progress=True, - ) diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py new file mode 100644 index 00000000..f5a28b2f --- /dev/null +++ b/src/virtualship/instruments/base.py @@ -0,0 +1,178 @@ +import abc +from datetime import timedelta +from pathlib import Path +from typing import TYPE_CHECKING + +import copernicusmarine +import numpy as np +import xarray as xr +from parcels import FieldSet +from yaspin import yaspin + +from virtualship.utils import ( + COPERNICUSMARINE_BGC_VARIABLES, + COPERNICUSMARINE_PHYS_VARIABLES, + _get_bathy_data, + _select_product_id, + ship_spinner, +) + +if TYPE_CHECKING: + from virtualship.models import Expedition + + +class Instrument(abc.ABC): + """Base class for instruments and their simulation.""" + + #! TODO List: + # TODO: how is this handling credentials?! Seems to work already, are these set up from my previous instances of using copernicusmarine? Therefore users will only have to do it once too? + + def __init__( + self, + name: str, + expedition: "Expedition", + directory: Path | str, + filenames: dict, + variables: dict, + add_bathymetry: bool, + allow_time_extrapolation: bool, + verbose_progress: bool, + buffer_spec: dict | None = None, + limit_spec: dict | None = None, + ): + """Initialise instrument.""" + self.name = name + self.expedition = expedition + self.directory = directory + self.filenames = filenames + self.variables = variables + self.dimensions = { + "lon": "longitude", + "lat": "latitude", + "time": "time", + "depth": "depth", + } # same dimensions for all instruments + self.add_bathymetry = add_bathymetry + self.allow_time_extrapolation = allow_time_extrapolation + self.verbose_progress = verbose_progress + self.buffer_spec = buffer_spec + self.limit_spec = limit_spec + + def load_input_data(self) -> FieldSet: + """Load and return the input data as a FieldSet for the instrument.""" + try: + datasets = [] + for var in self.variables.values(): + physical = True if var in COPERNICUSMARINE_PHYS_VARIABLES else False + datasets.append(self._get_copernicus_ds(physical=physical, var=var)) + + # make sure time dims are matched if BGC variables are present (different monthly/daily resolutions can impact fieldset_endtime in simulate) + all_keys = set().union(*(ds.keys() for ds in datasets)) + if all_keys & set(COPERNICUSMARINE_BGC_VARIABLES): + datasets = self._align_temporal(datasets) + + ds_concat = xr.merge(datasets) # TODO: deal with WARNINGS? + + fieldset = FieldSet.from_xarray_dataset( + ds_concat, self.variables, self.dimensions, mesh="spherical" + ) + + except Exception as e: + raise FileNotFoundError( + f"Failed to load input data directly from Copernicus Marine for instrument '{self.name}'. " + f"Please check your credentials, network connection, and variable names. Original error: {e}" + ) from e + + # interpolation methods + for var in (v for v in self.variables if v not in ("U", "V")): + getattr(fieldset, var).interp_method = "linear_invdist_land_tracer" + # depth negative + for g in fieldset.gridset.grids: + g.negate_depth() + + # bathymetry data + if self.add_bathymetry: + bathymetry_field = _get_bathy_data( + self.expedition.schedule.space_time_region, + latlon_buffer=self.buffer_spec.get("latlon") + if self.buffer_spec + else None, + ).bathymetry + bathymetry_field.data = -bathymetry_field.data + fieldset.add_field(bathymetry_field) + + return fieldset + + @abc.abstractmethod + def simulate(self, data_dir: Path, measurements: list, out_path: str | Path): + """Simulate instrument measurements.""" + + def execute(self, measurements: list, out_path: str | Path) -> None: + """Run instrument simulation.""" + if not self.verbose_progress: + with yaspin( + text=f"Simulating {self.name} measurements... ", + side="right", + spinner=ship_spinner, + ) as spinner: + self.simulate(measurements, out_path) + spinner.ok("✅\n") + else: + print(f"Simulating {self.name} measurements... ") + self.simulate(measurements, out_path) + print("\n") + + def _get_copernicus_ds( + self, + physical: bool, + var: str, + ) -> xr.Dataset: + """Get Copernicus Marine dataset for direct ingestion.""" + product_id = _select_product_id( + physical=physical, + schedule_start=self.expedition.schedule.space_time_region.time_range.start_time, + schedule_end=self.expedition.schedule.space_time_region.time_range.end_time, + variable=var if not physical else None, + ) + + latlon_buffer = self.buffer_spec.get("latlon") if self.buffer_spec else 0.0 + time_buffer = self.buffer_spec.get("time") if self.buffer_spec else 0.0 + + depth_min = self.limit_spec.get("depth_min") if self.limit_spec else None + depth_max = self.limit_spec.get("depth_max") if self.limit_spec else None + + return copernicusmarine.open_dataset( + dataset_id=product_id, + minimum_longitude=self.expedition.schedule.space_time_region.spatial_range.minimum_longitude + - latlon_buffer, + maximum_longitude=self.expedition.schedule.space_time_region.spatial_range.maximum_longitude + + latlon_buffer, + minimum_latitude=self.expedition.schedule.space_time_region.spatial_range.minimum_latitude + - latlon_buffer, + maximum_latitude=self.expedition.schedule.space_time_region.spatial_range.maximum_latitude + + latlon_buffer, + variables=[var], + start_datetime=self.expedition.schedule.space_time_region.time_range.start_time, + end_datetime=self.expedition.schedule.space_time_region.time_range.end_time + + timedelta(days=time_buffer), + minimum_depth=depth_min, + maximum_depth=depth_max, + coordinates_selection_method="outside", + ) + + def _align_temporal(self, datasets: list[xr.Dataset]) -> list[xr.Dataset]: + """Align monthly and daily time dims of multiple datasets (by repeating monthly values daily).""" + reference_time = datasets[ + np.argmax(ds.time for ds in datasets) + ].time # daily timeseries + + datasets_aligned = [] + for ds in datasets: + if not np.array_equal(ds.time, reference_time): + # TODO: NEED TO CHOOSE BEST METHOD HERE + # ds = ds.resample(time="1D").ffill().reindex(time=reference_time) + # ds = ds.resample(time="1D").ffill() + ds = ds.reindex({"time": reference_time}, method="nearest") + datasets_aligned.append(ds) + + return datasets_aligned diff --git a/src/virtualship/instruments/ctd.py b/src/virtualship/instruments/ctd.py index 41185007..5ec02e41 100644 --- a/src/virtualship/instruments/ctd.py +++ b/src/virtualship/instruments/ctd.py @@ -1,24 +1,36 @@ -"""CTD instrument.""" - from dataclasses import dataclass from datetime import timedelta -from pathlib import Path +from typing import TYPE_CHECKING, ClassVar import numpy as np -from parcels import FieldSet, JITParticle, ParticleSet, Variable +from parcels import JITParticle, ParticleSet, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType -from virtualship.models import Spacetime +if TYPE_CHECKING: + from virtualship.models.spacetime import Spacetime +from virtualship.utils import add_dummy_UV, register_instrument + +# ===================================================== +# SECTION: Dataclass +# ===================================================== @dataclass class CTD: - """Configuration for a single CTD.""" + """CTD configuration.""" - spacetime: Spacetime + name: ClassVar[str] = "CTD" + spacetime: "Spacetime" min_depth: float max_depth: float +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _CTDParticle = JITParticle.add_variables( [ Variable("salinity", dtype=np.float32, initial=np.nan), @@ -31,6 +43,11 @@ class CTD: ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + + def _sample_temperature(particle, fieldset, time): particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] @@ -53,85 +70,118 @@ def _ctd_cast(particle, fieldset, time): particle.delete() -def simulate_ctd( - fieldset: FieldSet, - out_path: str | Path, - ctds: list[CTD], - outputdt: timedelta, -) -> None: - """ - Use Parcels to simulate a set of CTDs in a fieldset. - - :param fieldset: The fieldset to simulate the CTDs in. - :param out_path: The path to write the results to. - :param ctds: A list of CTDs to simulate. - :param outputdt: Interval which dictates the update frequency of file output during simulation - :raises ValueError: Whenever provided 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(ctds) == 0: - print( - "No CTDs provided. Parcels currently crashes when providing an empty particle set, so no CTD simulation will be done and no files will be created." +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.CTD) +class CTDInstrument(Instrument): + """CTD instrument class.""" + + def __init__(self, expedition, directory): + """Initialize CTDInstrument.""" + filenames = { + "S": f"{CTD.name}_s.nc", + "T": f"{CTD.name}_t.nc", + } + variables = {"S": "so", "T": "thetao"} + + super().__init__( + CTD.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=True, + allow_time_extrapolation=True, + verbose_progress=False, + buffer_spec=None, + limit_spec=None, ) - # 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.spacetime.time) >= fieldset_starttime for ctd in ctds] - ): - raise ValueError("CTD deployed before fieldset starts.") - - # depth the ctd will go to. shallowest between ctd max depth and bathymetry. - max_depths = [ - max( - ctd.max_depth, - fieldset.bathymetry.eval( - z=0, y=ctd.spacetime.location.lat, x=ctd.spacetime.location.lon, time=0 - ), + + def simulate(self, measurements, out_path) -> None: + """Simulate CTD measurements.""" + WINCH_SPEED = 1.0 # sink and rise speed in m/s + DT = 10.0 # dt of CTD simulation integrator + OUTPUT_DT = timedelta(seconds=10) # output dt for CTD simulation + + if len(measurements) == 0: + print( + "No CTDs provided. Parcels currently crashes when providing an empty particle set, so no CTD simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + + fieldset = self.load_input_data() + + # add dummy U + add_dummy_UV(fieldset) # TODO: parcels v3 bodge; remove when parcels v4 is used + + fieldset_starttime = fieldset.T.grid.time_origin.fulltime( + fieldset.T.grid.time_full[0] + ) + fieldset_endtime = fieldset.T.grid.time_origin.fulltime( + fieldset.T.grid.time_full[-1] ) - for ctd in ctds - ] - # 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"CTD max_depth or bathymetry shallower than maximum {-DT * WINCH_SPEED}" + # deploy time for all ctds should be later than fieldset start time + if not all( + [ + np.datetime64(ctd.spacetime.time) >= fieldset_starttime + for ctd in measurements + ] + ): + raise ValueError("CTD deployed before fieldset starts.") + + # depth the ctd will go to. shallowest between ctd max depth and bathymetry. + max_depths = [ + max( + ctd.max_depth, + fieldset.bathymetry.eval( + z=0, + y=ctd.spacetime.location.lat, + x=ctd.spacetime.location.lon, + time=0, + ), + ) + for ctd in measurements + ] + + # 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"CTD max_depth or bathymetry shallower than maximum {-DT * WINCH_SPEED}" + ) + + # define parcel particles + ctd_particleset = ParticleSet( + fieldset=fieldset, + pclass=_CTDParticle, + lon=[ctd.spacetime.location.lon for ctd in measurements], + lat=[ctd.spacetime.location.lat for ctd in measurements], + depth=[ctd.min_depth for ctd in measurements], + time=[ctd.spacetime.time for ctd in measurements], + max_depth=max_depths, + min_depth=[ctd.min_depth for ctd in measurements], + winch_speed=[WINCH_SPEED for _ in measurements], ) - # define parcel particles - ctd_particleset = ParticleSet( - fieldset=fieldset, - pclass=_CTDParticle, - lon=[ctd.spacetime.location.lon for ctd in ctds], - lat=[ctd.spacetime.location.lat for ctd in ctds], - depth=[ctd.min_depth for ctd in ctds], - time=[ctd.spacetime.time for ctd in ctds], - max_depth=max_depths, - min_depth=[ctd.min_depth for ctd in ctds], - winch_speed=[WINCH_SPEED for _ in ctds], - ) - - # define output file for the simulation - out_file = ctd_particleset.ParticleFile(name=out_path, outputdt=outputdt) - - # execute simulation - ctd_particleset.execute( - [_sample_salinity, _sample_temperature, _ctd_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_particleset.particledata) != 0: - raise ValueError( - "Simulation ended before CTD resurfaced. This most likely means the field time dimension did not match the simulation time span." + # define output file for the simulation + out_file = ctd_particleset.ParticleFile(name=out_path, outputdt=OUTPUT_DT) + + # execute simulation + ctd_particleset.execute( + [_sample_salinity, _sample_temperature, _ctd_cast], + endtime=fieldset_endtime, + dt=DT, + verbose_progress=self.verbose_progress, + output_file=out_file, ) + + # there should be no particles left, as they delete themselves when they resurface + if len(ctd_particleset.particledata) != 0: + raise ValueError( + "Simulation ended before CTD resurfaced. This most likely means the field time dimension did not match the simulation time span." + ) diff --git a/src/virtualship/instruments/ctd_bgc.py b/src/virtualship/instruments/ctd_bgc.py index 0a34f61b..33e14aca 100644 --- a/src/virtualship/instruments/ctd_bgc.py +++ b/src/virtualship/instruments/ctd_bgc.py @@ -1,24 +1,34 @@ -"""CTD_BGC instrument.""" - from dataclasses import dataclass from datetime import timedelta -from pathlib import Path +from typing import ClassVar import numpy as np -from parcels import FieldSet, JITParticle, ParticleSet, Variable +from parcels import JITParticle, ParticleSet, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.models.spacetime import Spacetime +from virtualship.utils import add_dummy_UV, register_instrument -from virtualship.models import Spacetime +# ===================================================== +# SECTION: Dataclass +# ===================================================== @dataclass class CTD_BGC: - """Configuration for a single BGC CTD.""" + """CTD_BGC configuration.""" + name: ClassVar[str] = "CTD_BGC" spacetime: Spacetime min_depth: float max_depth: float +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _CTD_BGCParticle = JITParticle.add_variables( [ Variable("o2", dtype=np.float32, initial=np.nan), @@ -35,6 +45,10 @@ class CTD_BGC: ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + def _sample_o2(particle, fieldset, time): particle.o2 = fieldset.o2[time, particle.depth, particle.lat, particle.lon] @@ -78,100 +92,141 @@ def _ctd_bgc_cast(particle, fieldset, time): 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." +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.CTD_BGC) +class CTD_BGCInstrument(Instrument): + """CTD_BGC instrument class.""" + + def __init__(self, expedition, directory): + """Initialize CTD_BGCInstrument.""" + filenames = { + "o2": f"{CTD_BGC.name}_o2.nc", + "chl": f"{CTD_BGC.name}_chl.nc", + "no3": f"{CTD_BGC.name}_no3.nc", + "po4": f"{CTD_BGC.name}_po4.nc", + "ph": f"{CTD_BGC.name}_ph.nc", + "phyc": f"{CTD_BGC.name}_phyc.nc", + "nppv": f"{CTD_BGC.name}_nppv.nc", + } + variables = { + "o2": "o2", + "chl": "chl", + "no3": "no3", + "po4": "po4", + "ph": "ph", + "phyc": "phyc", + "nppv": "nppv", + } + super().__init__( + CTD_BGC.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=True, + allow_time_extrapolation=True, + verbose_progress=False, + buffer_spec=None, + limit_spec=None, ) - # 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]) + def simulate(self, measurements, out_path) -> None: + """Simulate BGC CTD measurements using Parcels.""" + WINCH_SPEED = 1.0 # sink and rise speed in m/s + DT = 10.0 # dt of CTD_BGC simulation integrator + OUTPUT_DT = timedelta(seconds=10) # output dt for CTD_BGC simulation - # 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, - ), + if len(measurements) == 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 = self.load_input_data() + + # add dummy U + add_dummy_UV(fieldset) # TODO: parcels v3 bodge; remove when parcels v4 is used + + fieldset_starttime = fieldset.o2.grid.time_origin.fulltime( + fieldset.o2.grid.time_full[0] + ) + fieldset_endtime = fieldset.o2.grid.time_origin.fulltime( + fieldset.o2.grid.time_full[-1] ) - 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}" + # 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 measurements + ] + ): + 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 measurements + ] + + # 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 measurements], + lat=[ctd_bgc.spacetime.location.lat for ctd_bgc in measurements], + depth=[ctd_bgc.min_depth for ctd_bgc in measurements], + time=[ctd_bgc.spacetime.time for ctd_bgc in measurements], + max_depth=max_depths, + min_depth=[ctd_bgc.min_depth for ctd_bgc in measurements], + winch_speed=[WINCH_SPEED for _ in measurements], ) - # 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, - _sample_nitrate, - _sample_phosphate, - _sample_ph, - _sample_phytoplankton, - _sample_primary_production, - _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." + # define output file for the simulation + out_file = ctd_bgc_particleset.ParticleFile(name=out_path, outputdt=OUTPUT_DT) + + breakpoint() + + # execute simulation + ctd_bgc_particleset.execute( + [ + _sample_o2, + _sample_chlorophyll, + _sample_nitrate, + _sample_phosphate, + _sample_ph, + _sample_phytoplankton, + _sample_primary_production, + _ctd_bgc_cast, + ], + endtime=fieldset_endtime, + dt=DT, + verbose_progress=self.verbose_progress, + 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/instruments/drifter.py b/src/virtualship/instruments/drifter.py index 5aef240f..035a81b0 100644 --- a/src/virtualship/instruments/drifter.py +++ b/src/virtualship/instruments/drifter.py @@ -1,24 +1,34 @@ -"""Drifter instrument.""" - from dataclasses import dataclass -from datetime import datetime, timedelta -from pathlib import Path +from datetime import timedelta +from typing import ClassVar import numpy as np -from parcels import AdvectionRK4, FieldSet, JITParticle, ParticleSet, Variable +from parcels import AdvectionRK4, JITParticle, ParticleSet, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.models.spacetime import Spacetime +from virtualship.utils import register_instrument -from virtualship.models import Spacetime +# ===================================================== +# SECTION: Dataclass +# ===================================================== @dataclass class Drifter: - """Configuration for a single Drifter.""" + """Drifter configuration.""" + name: ClassVar[str] = "Drifter" spacetime: Spacetime depth: float # depth at which it floats and samples lifetime: timedelta | None # if none, lifetime is infinite +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _DrifterParticle = JITParticle.add_variables( [ Variable("temperature", dtype=np.float32, initial=np.nan), @@ -28,6 +38,10 @@ class Drifter: ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + def _sample_temperature(particle, fieldset, time): particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] @@ -40,74 +54,107 @@ def _check_lifetime(particle, fieldset, time): particle.delete() -def simulate_drifters( - fieldset: FieldSet, - out_path: str | Path, - drifters: list[Drifter], - outputdt: timedelta, - dt: timedelta, - endtime: datetime | None = None, -) -> None: - """ - Use Parcels to simulate a set of drifters in a fieldset. - - :param fieldset: The fieldset to simulate the Drifters in. - :param out_path: The path to write the results to. - :param drifters: A list of drifters to simulate. - :param outputdt: Interval which dictates the update frequency of file output during simulation. - :param dt: Dt for integration. - :param endtime: Stop at this time, or if None, continue until the end of the fieldset or until all drifters ended. If this is earlier than the last drifter ended or later than the end of the fieldset, a warning will be printed. - """ - if len(drifters) == 0: - print( - "No drifters provided. Parcels currently crashes when providing an empty particle set, so no drifter simulation will be done and no files will be created." +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.DRIFTER) +class DrifterInstrument(Instrument): + """Drifter instrument class.""" + + def __init__(self, expedition, directory): + """Initialize DrifterInstrument.""" + filenames = { + "U": f"{Drifter.name}_uv.nc", + "V": f"{Drifter.name}_uv.nc", + "T": f"{Drifter.name}_t.nc", + } + variables = {"U": "uo", "V": "vo", "T": "thetao"} + buffer_spec = { + "latlon": 6.0, # [degrees] + "time": 21.0, # [days] + } + limit_spec = { + "depth_min": 1.0, # [meters] + "depth_max": 1.0, # [meters] + } + + super().__init__( + Drifter.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=False, + allow_time_extrapolation=False, + verbose_progress=True, + buffer_spec=buffer_spec, + limit_spec=limit_spec, + ) + + def simulate(self, measurements, out_path) -> None: + """Simulate Drifter measurements.""" + OUTPUT_DT = timedelta(hours=5) + DT = timedelta(minutes=5) + ENDTIME = None + + if len(measurements) == 0: + print( + "No drifters provided. Parcels currently crashes when providing an empty particle set, so no drifter simulation will be done and no files will be created." + ) + # TODO when Parcels supports it this check can be removed. + return + + fieldset = self.load_input_data() + + # define parcel particles + drifter_particleset = ParticleSet( + fieldset=fieldset, + pclass=_DrifterParticle, + lat=[drifter.spacetime.location.lat for drifter in measurements], + lon=[drifter.spacetime.location.lon for drifter in measurements], + depth=[drifter.depth for drifter in measurements], + time=[drifter.spacetime.time for drifter in measurements], + has_lifetime=[ + 1 if drifter.lifetime is not None else 0 for drifter in measurements + ], + lifetime=[ + 0 if drifter.lifetime is None else drifter.lifetime.total_seconds() + for drifter in measurements + ], ) - # TODO when Parcels supports it this check can be removed. - return - - # define parcel particles - drifter_particleset = ParticleSet( - fieldset=fieldset, - pclass=_DrifterParticle, - lat=[drifter.spacetime.location.lat for drifter in drifters], - lon=[drifter.spacetime.location.lon for drifter in drifters], - depth=[drifter.depth for drifter in drifters], - time=[drifter.spacetime.time for drifter in drifters], - has_lifetime=[1 if drifter.lifetime is not None else 0 for drifter in drifters], - lifetime=[ - 0 if drifter.lifetime is None else drifter.lifetime.total_seconds() - for drifter in drifters - ], - ) - - # define output file for the simulation - out_file = drifter_particleset.ParticleFile( - name=out_path, outputdt=outputdt, chunks=[len(drifter_particleset), 100] - ) - - # get earliest between fieldset end time and provide end time - fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) - if endtime is None: - actual_endtime = fieldset_endtime - elif endtime > fieldset_endtime: - print("WARN: Requested end time later than fieldset end time.") - actual_endtime = fieldset_endtime - else: - actual_endtime = np.timedelta64(endtime) - - # execute simulation - drifter_particleset.execute( - [AdvectionRK4, _sample_temperature, _check_lifetime], - endtime=actual_endtime, - dt=dt, - output_file=out_file, - verbose_progress=True, - ) - - # if there are more particles left than the number of drifters with an indefinite endtime, warn the user - if len(drifter_particleset.particledata) > len( - [d for d in drifters if d.lifetime is None] - ): - print( - "WARN: Some drifters had a life time beyond the end time of the fieldset or the requested end time." + + # define output file for the simulation + out_file = drifter_particleset.ParticleFile( + name=out_path, + outputdt=OUTPUT_DT, + chunks=[len(drifter_particleset), 100], ) + + # get earliest between fieldset end time and prescribed end time + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) + if ENDTIME is None: + actual_endtime = fieldset_endtime + elif ENDTIME > fieldset_endtime: + print("WARN: Requested end time later than fieldset end time.") + actual_endtime = fieldset_endtime + else: + actual_endtime = np.timedelta64(ENDTIME) + + # execute simulation + drifter_particleset.execute( + [AdvectionRK4, _sample_temperature, _check_lifetime], + endtime=actual_endtime, + dt=DT, + output_file=out_file, + verbose_progress=self.verbose_progress, + ) + + # if there are more particles left than the number of drifters with an indefinite endtime, warn the user + if len(drifter_particleset.particledata) > len( + [d for d in measurements if d.lifetime is None] + ): + print( + "WARN: Some drifters had a life time beyond the end time of the fieldset or the requested end time." + ) diff --git a/src/virtualship/instruments/ship_underwater_st.py b/src/virtualship/instruments/ship_underwater_st.py index 7b08ad4b..1067c03e 100644 --- a/src/virtualship/instruments/ship_underwater_st.py +++ b/src/virtualship/instruments/ship_underwater_st.py @@ -1,14 +1,29 @@ -"""Ship salinity and temperature.""" - -from pathlib import Path +from dataclasses import dataclass +from typing import ClassVar import numpy as np -from parcels import FieldSet, ParticleSet, ScipyParticle, Variable +from parcels import ParticleSet, ScipyParticle, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.utils import add_dummy_UV, register_instrument + +# ===================================================== +# SECTION: Dataclass +# ===================================================== + + +@dataclass +class Underwater_ST: + """Underwater_ST configuration.""" + + name: ClassVar[str] = "Underwater_ST" -from virtualship.models import Spacetime -# we specifically use ScipyParticle because we have many small calls to execute -# there is some overhead with JITParticle and this ends up being significantly faster +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _ShipSTParticle = ScipyParticle.add_variables( [ Variable("S", dtype=np.float32, initial=np.nan), @@ -16,6 +31,10 @@ ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + # define function sampling Salinity def _sample_salinity(particle, fieldset, time): @@ -27,50 +46,69 @@ def _sample_temperature(particle, fieldset, time): particle.T = fieldset.T[time, particle.depth, particle.lat, particle.lon] -def simulate_ship_underwater_st( - fieldset: FieldSet, - out_path: str | Path, - depth: float, - sample_points: list[Spacetime], -) -> None: - """ - Use Parcels to simulate underway data, measuring salinity and temperature at the given depth along the ship track in a fieldset. - - :param fieldset: The fieldset to simulate the sampling in. - :param out_path: The path to write the results to. - :param depth: The depth at which to measure. 0 is water surface, negative is into the water. - :param sample_points: The places and times to sample at. - """ - sample_points.sort(key=lambda p: p.time) - - particleset = ParticleSet.from_list( - fieldset=fieldset, - pclass=_ShipSTParticle, - lon=0.0, # initial lat/lon are irrelevant and will be overruled later - lat=0.0, - depth=depth, - time=0, # same for time - ) - - # define output file for the simulation - # outputdt set to infinie as we want to just want to write at the end of every call to 'execute' - out_file = particleset.ParticleFile(name=out_path, outputdt=np.inf) - - # iterate over each point, manually set lat lon time, then - # execute the particle set for one step, performing one set of measurement - for point in sample_points: - particleset.lon_nextloop[:] = point.location.lon - particleset.lat_nextloop[:] = point.location.lat - particleset.time_nextloop[:] = fieldset.time_origin.reltime( - np.datetime64(point.time) - ) +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.UNDERWATER_ST) +class Underwater_STInstrument(Instrument): + """Underwater_ST instrument class.""" - # perform one step using the particleset - # dt and runtime are set so exactly one step is made. - particleset.execute( - [_sample_salinity, _sample_temperature], - dt=1, - runtime=1, + def __init__(self, expedition, directory): + """Initialize Underwater_STInstrument.""" + filenames = { + "S": f"{Underwater_ST.name}_s.nc", + "T": f"{Underwater_ST.name}_t.nc", + } + variables = {"S": "so", "T": "thetao"} + + super().__init__( + Underwater_ST.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=False, + allow_time_extrapolation=True, verbose_progress=False, - output_file=out_file, + buffer_spec=None, + limit_spec=None, + ) + + def simulate(self, measurements, out_path) -> None: + """Simulate underway salinity and temperature measurements.""" + DEPTH = -2.0 + + measurements.sort(key=lambda p: p.time) + + fieldset = self.load_input_data() + + # add dummy U + add_dummy_UV(fieldset) # TODO: parcels v3 bodge; remove when parcels v4 is used + + particleset = ParticleSet.from_list( + fieldset=fieldset, + pclass=_ShipSTParticle, + lon=0.0, + lat=0.0, + depth=DEPTH, + time=0, ) + + out_file = particleset.ParticleFile(name=out_path, outputdt=np.inf) + + for point in measurements: + particleset.lon_nextloop[:] = point.location.lon + particleset.lat_nextloop[:] = point.location.lat + particleset.time_nextloop[:] = fieldset.time_origin.reltime( + np.datetime64(point.time) + ) + + particleset.execute( + [_sample_salinity, _sample_temperature], + dt=1, + runtime=1, + verbose_progress=self.verbose_progress, + output_file=out_file, + ) diff --git a/src/virtualship/instruments/types.py b/src/virtualship/instruments/types.py new file mode 100644 index 00000000..9ae221e9 --- /dev/null +++ b/src/virtualship/instruments/types.py @@ -0,0 +1,18 @@ +from enum import Enum + + +class InstrumentType(Enum): + """Types of the instruments.""" + + CTD = "CTD" + CTD_BGC = "CTD_BGC" + DRIFTER = "DRIFTER" + ARGO_FLOAT = "ARGO_FLOAT" + XBT = "XBT" + ADCP = "ADCP" + UNDERWATER_ST = "UNDERWATER_ST" + + @property + def is_underway(self) -> bool: + """Return True if instrument is an underway instrument (ADCP, UNDERWATER_ST).""" + return self in {InstrumentType.ADCP, InstrumentType.UNDERWATER_ST} diff --git a/src/virtualship/instruments/xbt.py b/src/virtualship/instruments/xbt.py index 6d75be8c..9ec76527 100644 --- a/src/virtualship/instruments/xbt.py +++ b/src/virtualship/instruments/xbt.py @@ -1,19 +1,25 @@ -"""XBT instrument.""" - from dataclasses import dataclass from datetime import timedelta -from pathlib import Path +from typing import ClassVar import numpy as np -from parcels import FieldSet, JITParticle, ParticleSet, Variable +from parcels import JITParticle, ParticleSet, Variable + +from virtualship.instruments.base import Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.models.spacetime import Spacetime +from virtualship.utils import register_instrument -from virtualship.models import Spacetime +# ===================================================== +# SECTION: Dataclass +# ===================================================== @dataclass class XBT: - """Configuration for a single XBT.""" + """XBT configuration.""" + name: ClassVar[str] = "XBT" spacetime: Spacetime min_depth: float max_depth: float @@ -21,6 +27,10 @@ class XBT: deceleration_coefficient: float +# ===================================================== +# SECTION: Particle Class +# ===================================================== + _XBTParticle = JITParticle.add_variables( [ Variable("temperature", dtype=np.float32, initial=np.nan), @@ -31,6 +41,10 @@ class XBT: ] ) +# ===================================================== +# SECTION: Kernels +# ===================================================== + def _sample_temperature(particle, fieldset, time): particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] @@ -54,88 +68,112 @@ def _xbt_cast(particle, fieldset, time): particle_ddepth = particle.max_depth - particle.depth -def simulate_xbt( - fieldset: FieldSet, - out_path: str | Path, - xbts: list[XBT], - outputdt: timedelta, -) -> None: - """ - Use Parcels to simulate a set of XBTs in a fieldset. - - :param fieldset: The fieldset to simulate the XBTs in. - :param out_path: The path to write the results to. - :param xbts: A list of XBTs to simulate. - :param outputdt: Interval which dictates the update frequency of file output during simulation - :raises ValueError: Whenever provided XBTs, fieldset, are not compatible with this function. - """ - DT = 10.0 # dt of XBT simulation integrator - - if len(xbts) == 0: - print( - "No XBTs provided. Parcels currently crashes when providing an empty particle set, so no XBT simulation will be done and no files will be created." +# ===================================================== +# SECTION: Instrument Class +# ===================================================== + + +@register_instrument(InstrumentType.XBT) +class XBTInstrument(Instrument): + """XBT instrument class.""" + + def __init__(self, expedition, directory): + """Initialize XBTInstrument.""" + filenames = { + "U": f"{XBT.name}_uv.nc", + "V": f"{XBT.name}_uv.nc", + "S": f"{XBT.name}_s.nc", + "T": f"{XBT.name}_t.nc", + } + variables = {"U": "uo", "V": "vo", "S": "so", "T": "thetao"} + super().__init__( + XBT.name, + expedition, + directory, + filenames, + variables, + add_bathymetry=True, + allow_time_extrapolation=True, + verbose_progress=False, + buffer_spec=None, + limit_spec=None, ) - # 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 xbts should be later than fieldset start time - if not all( - [np.datetime64(xbt.spacetime.time) >= fieldset_starttime for xbt in xbts] - ): - raise ValueError("XBT deployed before fieldset starts.") - - # depth the xbt will go to. shallowest between xbt max depth and bathymetry. - max_depths = [ - max( - xbt.max_depth, - fieldset.bathymetry.eval( - z=0, y=xbt.spacetime.location.lat, x=xbt.spacetime.location.lon, time=0 - ), - ) - for xbt in xbts - ] - # initial fall speeds - initial_fall_speeds = [xbt.fall_speed for xbt in xbts] + def simulate(self, measurements, out_path) -> None: + """Simulate XBT measurements.""" + DT = 10.0 # dt of XBT simulation integrator + OUTPUT_DT = timedelta(seconds=1) - # XBT depth can not be too shallow, because kernel would break. - # This shallow is not useful anyway, no need to support. - for max_depth, fall_speed in zip(max_depths, initial_fall_speeds, strict=False): - if not max_depth <= -DT * fall_speed: - raise ValueError( - f"XBT max_depth or bathymetry shallower than maximum {-DT * fall_speed}" + if len(measurements) == 0: + print( + "No XBTs provided. Parcels currently crashes when providing an empty particle set, so no XBT simulation will be done and no files will be created." ) + # TODO when Parcels supports it this check can be removed. + return + + fieldset = self.load_input_data() + + 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 xbts should be later than fieldset start time + if not all( + [ + np.datetime64(xbt.spacetime.time) >= fieldset_starttime + for xbt in measurements + ] + ): + raise ValueError("XBT deployed before fieldset starts.") + + # depth the xbt will go to. shallowest between xbt max depth and bathymetry. + max_depths = [ + max( + xbt.max_depth, + fieldset.bathymetry.eval( + z=0, + y=xbt.spacetime.location.lat, + x=xbt.spacetime.location.lon, + time=0, + ), + ) + for xbt in measurements + ] + + # initial fall speeds + initial_fall_speeds = [xbt.fall_speed for xbt in measurements] + + # XBT depth can not be too shallow, because kernel would break. + for max_depth, fall_speed in zip(max_depths, initial_fall_speeds, strict=False): + if not max_depth <= -DT * fall_speed: + raise ValueError( + f"XBT max_depth or bathymetry shallower than minimum {-DT * fall_speed}. It is likely the XBT cannot be deployed in this area, which is too shallow." + ) + + # define xbt particles + xbt_particleset = ParticleSet( + fieldset=fieldset, + pclass=_XBTParticle, + lon=[xbt.spacetime.location.lon for xbt in measurements], + lat=[xbt.spacetime.location.lat for xbt in measurements], + depth=[xbt.min_depth for xbt in measurements], + time=[xbt.spacetime.time for xbt in measurements], + max_depth=max_depths, + min_depth=[xbt.min_depth for xbt in measurements], + fall_speed=[xbt.fall_speed for xbt in measurements], + ) - # define xbt particles - xbt_particleset = ParticleSet( - fieldset=fieldset, - pclass=_XBTParticle, - lon=[xbt.spacetime.location.lon for xbt in xbts], - lat=[xbt.spacetime.location.lat for xbt in xbts], - depth=[xbt.min_depth for xbt in xbts], - time=[xbt.spacetime.time for xbt in xbts], - max_depth=max_depths, - min_depth=[xbt.min_depth for xbt in xbts], - fall_speed=[xbt.fall_speed for xbt in xbts], - ) - - # define output file for the simulation - out_file = xbt_particleset.ParticleFile(name=out_path, outputdt=outputdt) - - # execute simulation - xbt_particleset.execute( - [_sample_temperature, _xbt_cast], - endtime=fieldset_endtime, - dt=DT, - verbose_progress=False, - output_file=out_file, - ) + out_file = xbt_particleset.ParticleFile(name=out_path, outputdt=OUTPUT_DT) - # there should be no particles left, as they delete themselves when they finish profiling - if len(xbt_particleset.particledata) != 0: - raise ValueError( - "Simulation ended before XBT finished profiling. This most likely means the field time dimension did not match the simulation time span." + xbt_particleset.execute( + [_sample_temperature, _xbt_cast], + endtime=fieldset_endtime, + dt=DT, + verbose_progress=self.verbose_progress, + output_file=out_file, ) + + # there should be no particles left, as they delete themselves when they finish profiling + if len(xbt_particleset.particledata) != 0: + raise ValueError( + "Simulation ended before XBT finished profiling. This most likely means the field time dimension did not match the simulation time span." + ) diff --git a/src/virtualship/models/__init__.py b/src/virtualship/models/__init__.py index a2f1546c..5eaabb85 100644 --- a/src/virtualship/models/__init__.py +++ b/src/virtualship/models/__init__.py @@ -8,7 +8,6 @@ DrifterConfig, Expedition, InstrumentsConfig, - InstrumentType, Schedule, ShipConfig, ShipUnderwaterSTConfig, @@ -30,7 +29,6 @@ "Schedule", "ShipConfig", "Waypoint", - "InstrumentType", "ArgoFloatConfig", "ADCPConfig", "CTDConfig", diff --git a/src/virtualship/expedition/checkpoint.py b/src/virtualship/models/checkpoint.py similarity index 95% rename from src/virtualship/expedition/checkpoint.py rename to src/virtualship/models/checkpoint.py index 6daf1a9b..98fe1ae0 100644 --- a/src/virtualship/expedition/checkpoint.py +++ b/src/virtualship/models/checkpoint.py @@ -8,7 +8,8 @@ import yaml from virtualship.errors import CheckpointError -from virtualship.models import InstrumentType, Schedule +from virtualship.instruments.types import InstrumentType +from virtualship.models import Schedule class _YamlDumper(yaml.SafeDumper): diff --git a/src/virtualship/models/expedition.py b/src/virtualship/models/expedition.py index 2e073b84..0cb13955 100644 --- a/src/virtualship/models/expedition.py +++ b/src/virtualship/models/expedition.py @@ -2,23 +2,21 @@ import itertools from datetime import datetime, timedelta -from enum import Enum from typing import TYPE_CHECKING import pydantic import pyproj import yaml -from virtualship.errors import ConfigError, ScheduleError -from virtualship.utils import _validate_numeric_mins_to_timedelta +from virtualship.errors import InstrumentsConfigError, ScheduleError +from virtualship.instruments.types import InstrumentType +from virtualship.utils import _get_bathy_data, _validate_numeric_mins_to_timedelta from .location import Location from .space_time_region import SpaceTimeRegion if TYPE_CHECKING: - from parcels import FieldSet - - from virtualship.expedition.input_data import InputData + pass projection: pyproj.Geod = pyproj.Geod(ellps="WGS84") @@ -45,6 +43,28 @@ def from_yaml(cls, file_path: str) -> Expedition: data = yaml.safe_load(file) return Expedition(**data) + def get_instruments(self) -> set[InstrumentType]: + """Return a set of unique InstrumentType enums used in the expedition.""" + instruments_in_expedition = [] + # from waypoints + for waypoint in self.schedule.waypoints: + if waypoint.instrument: + for instrument in waypoint.instrument: + if instrument: + instruments_in_expedition.append(instrument) + + # check for underway instruments and add if present in expeditions + try: + if self.instruments_config.adcp_config is not None: + instruments_in_expedition.append(InstrumentType.ADCP) + if self.instruments_config.ship_underwater_st_config is not None: + instruments_in_expedition.append(InstrumentType.UNDERWATER_ST) + return sorted(set(instruments_in_expedition), key=lambda x: x.name) + except Exception as e: + raise InstrumentsConfigError( + "Underway instrument config attribute(s) are missing from YAML. Must be Config object or None." + ) from e + class ShipConfig(pydantic.BaseModel): """Configuration of the ship.""" @@ -64,23 +84,12 @@ class Schedule(pydantic.BaseModel): model_config = pydantic.ConfigDict(extra="forbid") - def get_instruments(self) -> set[InstrumentType]: - """Return a set of unique InstrumentType enums used in the schedule.""" - instruments_in_schedule = [] - for waypoint in self.waypoints: - if waypoint.instrument: - for instrument in waypoint.instrument: - if instrument: - instruments_in_schedule.append(instrument) - return set(instruments_in_schedule) - def verify( self, ship_speed: float, - input_data: InputData | None, + ignore_missing_bathymetry: bool = False, *, check_space_time_region: bool = False, - ignore_missing_fieldsets: bool = False, ) -> None: """ Verify the feasibility and correctness of the schedule's waypoints. @@ -91,20 +100,12 @@ def verify( 3. Waypoint times are in ascending order. 4. All waypoints are in water (not on land). 5. The ship can arrive on time at each waypoint given its speed. - - :param ship_speed: The ship's speed in knots. - :param input_data: An InputData object containing fieldsets used to check if waypoints are on water. - :param check_space_time_region: whether to check for missing space_time_region. - :param ignore_missing_fieldsets: whether to ignore warning for missing field sets. - :raises PlanningError: If any of the verification checks fail, indicating infeasible or incorrect waypoints. - :raises NotImplementedError: If an instrument in the schedule is not implemented. - :return: None. The method doesn't return a value but raises exceptions if verification fails. """ print("\nVerifying route... ") if check_space_time_region and self.space_time_region is None: raise ScheduleError( - "space_time_region not found in schedule, please define it to fetch the data." + "space_time_region not found in schedule, please define it to proceed." ) if len(self.waypoints) == 0: @@ -125,45 +126,33 @@ def verify( f"Waypoint(s) {', '.join(f'#{i + 1}' for i in invalid_i)}: each waypoint should be timed after all previous waypoints", ) - # check if all waypoints are in water - # this is done by picking an arbitrary provided fieldset and checking if UV is not zero - - # get all available fieldsets - available_fieldsets = [] - if input_data is not None: - fieldsets = [ - input_data.adcp_fieldset, - input_data.argo_float_fieldset, - input_data.ctd_fieldset, - input_data.drifter_fieldset, - input_data.ship_underwater_st_fieldset, - ] - for fs in fieldsets: - if fs is not None: - available_fieldsets.append(fs) - - # check if there are any fieldsets, else it's an error - if len(available_fieldsets) == 0: - if not ignore_missing_fieldsets: - print( - "Cannot verify because no fieldsets have been loaded. This is probably " - "because you are not using any instruments in your schedule. This is not a problem, " - "but carefully check your waypoint locations manually." - ) + # check if all waypoints are in water using bathymetry data + # TODO: write test that checks that will flag when waypoint is on land!! [add to existing suite of fail .verify() tests in test_expedition.py] + land_waypoints = [] + if not ignore_missing_bathymetry: + try: + bathymetry_field = _get_bathy_data( + self.space_time_region, latlon_buffer=None + ).bathymetry # via copernicusmarine + except Exception as e: + raise ScheduleError( + f"Problem loading bathymetry data (used to verify waypoints are in water) directly via copernicusmarine. \n\n original message: {e}" + ) from e + + for wp_i, wp in enumerate(self.waypoints): + try: + bathymetry_field.eval( + 0, # time + 0, # depth (surface) + wp.location.lat, + wp.location.lon, + ) + except Exception: + land_waypoints.append((wp_i, wp)) - else: - # pick any - fieldset = available_fieldsets[0] - # get waypoints with 0 UV - land_waypoints = [ - (wp_i, wp) - for wp_i, wp in enumerate(self.waypoints) - if _is_on_land_zero_uv(fieldset, wp) - ] - # raise an error if there are any if len(land_waypoints) > 0: raise ScheduleError( - f"The following waypoints are on land: {['#' + str(wp_i) + ' ' + str(wp) for (wp_i, wp) in land_waypoints]}" + f"The following waypoint(s) throw(s) error(s): {['#' + str(wp_i + 1) + ' ' + str(wp) for (wp_i, wp) in land_waypoints]}\n\nINFO: They are likely on land (bathymetry data cannot be interpolated to their location(s)).\n" ) # check that ship will arrive on time at each waypoint (in case no unexpected event happen) @@ -213,16 +202,6 @@ def serialize_instrument(self, instrument): return instrument.value if instrument else None -class InstrumentType(Enum): - """Types of the instruments.""" - - CTD = "CTD" - CTD_BGC = "CTD_BGC" - DRIFTER = "DRIFTER" - ARGO_FLOAT = "ARGO_FLOAT" - XBT = "XBT" - - class ArgoFloatConfig(pydantic.BaseModel): """Configuration for argos floats.""" @@ -404,53 +383,38 @@ class InstrumentsConfig(pydantic.BaseModel): model_config = pydantic.ConfigDict(extra="forbid") - def verify(self, schedule: Schedule) -> None: + def verify(self, expedition: Expedition) -> None: """ Verify instrument configurations against the schedule. Removes instrument configs not present in the schedule and checks that all scheduled instruments are configured. Raises ConfigError if any scheduled instrument is missing a config. """ - instruments_in_schedule = schedule.get_instruments() + instruments_in_expedition = expedition.get_instruments() instrument_config_map = { InstrumentType.ARGO_FLOAT: "argo_float_config", InstrumentType.DRIFTER: "drifter_config", InstrumentType.XBT: "xbt_config", InstrumentType.CTD: "ctd_config", InstrumentType.CTD_BGC: "ctd_bgc_config", + InstrumentType.ADCP: "adcp_config", + InstrumentType.UNDERWATER_ST: "ship_underwater_st_config", } # Remove configs for unused instruments for inst_type, config_attr in instrument_config_map.items(): - if hasattr(self, config_attr) and inst_type not in instruments_in_schedule: - print( - f"{inst_type.value} configuration provided but not in schedule. Removing config." - ) + if ( + hasattr(self, config_attr) + and inst_type not in instruments_in_expedition + ): setattr(self, config_attr, None) # Check all scheduled instruments are configured - for inst_type in instruments_in_schedule: + for inst_type in instruments_in_expedition: config_attr = instrument_config_map.get(inst_type) if ( not config_attr or not hasattr(self, config_attr) or getattr(self, config_attr) is None ): - raise ConfigError( - f"Schedule includes instrument '{inst_type.value}', but instruments_config does not provide configuration for it." + raise InstrumentsConfigError( + f"Expedition includes instrument '{inst_type.value}', but instruments_config does not provide configuration for it." ) - - -def _is_on_land_zero_uv(fieldset: FieldSet, waypoint: Waypoint) -> bool: - """ - Check if waypoint is on land by assuming zero velocity means land. - - :param fieldset: The fieldset to sample the velocity from. - :param waypoint: The waypoint to check. - :returns: If the waypoint is on land. - """ - return fieldset.UV.eval( - 0, - fieldset.gridset.grids[0].depth[0], - waypoint.location.lat, - waypoint.location.lon, - applyConversion=False, - ) == (0.0, 0.0) diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py index 0a39d035..22e8f498 100644 --- a/src/virtualship/utils.py +++ b/src/virtualship/utils.py @@ -8,9 +8,17 @@ from pathlib import Path from typing import TYPE_CHECKING, TextIO +import copernicusmarine +import numpy as np +from parcels import FieldSet + +from virtualship.errors import CopernicusCatalogueError + if TYPE_CHECKING: + from virtualship.expedition.simulate_schedule import ScheduleOk from virtualship.models import Expedition + import pandas as pd import yaml from pydantic import BaseModel @@ -239,3 +247,212 @@ def _get_expedition(expedition_dir: Path) -> Expedition: "🚢 ", ], ) + + +# InstrumentType -> Instrument registry and registration utilities. +INSTRUMENT_CLASS_MAP = {} + + +def register_instrument(instrument_type): + def decorator(cls): + INSTRUMENT_CLASS_MAP[instrument_type] = cls + return cls + + return decorator + + +def get_instrument_class(instrument_type): + return INSTRUMENT_CLASS_MAP.get(instrument_type) + + +def add_dummy_UV(fieldset: FieldSet): + """Add a dummy U and V field to a FieldSet to satisfy parcels FieldSet completeness checks.""" + if "U" not in fieldset.__dict__.keys(): + for uv_var in ["U", "V"]: + dummy_field = getattr( + FieldSet.from_data( + {"U": 0, "V": 0}, {"lon": 0, "lat": 0}, mesh="spherical" + ), + uv_var, + ) + fieldset.add_field(dummy_field) + try: + fieldset.time_origin = ( + fieldset.T.grid.time_origin + if "T" in fieldset.__dict__.keys() + else fieldset.o2.grid.time_origin + ) + except Exception: + raise ValueError( + "Cannot determine time_origin for dummy UV fields. Assert T or o2 exists in fieldset." + ) from None + + +# Copernicus Marine product IDs + +PRODUCT_IDS = { + "phys": { + "reanalysis": "cmems_mod_glo_phy_my_0.083deg_P1D-m", + "reanalysis_interim": "cmems_mod_glo_phy_myint_0.083deg_P1D-m", + "analysis": "cmems_mod_glo_phy_anfc_0.083deg_P1D-m", + }, + "bgc": { + "reanalysis": "cmems_mod_glo_bgc_my_0.25deg_P1D-m", + "reanalysis_interim": "cmems_mod_glo_bgc_myint_0.25deg_P1D-m", + "analysis": None, # will be set per variable + }, +} + +BGC_ANALYSIS_IDS = { + "o2": "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1D-m", + "chl": "cmems_mod_glo_bgc-pft_anfc_0.25deg_P1D-m", + "no3": "cmems_mod_glo_bgc-nut_anfc_0.25deg_P1D-m", + "po4": "cmems_mod_glo_bgc-nut_anfc_0.25deg_P1D-m", + "ph": "cmems_mod_glo_bgc-car_anfc_0.25deg_P1D-m", + "phyc": "cmems_mod_glo_bgc-pft_anfc_0.25deg_P1D-m", + "nppv": "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1D-m", +} + +MONTHLY_BGC_REANALYSIS_IDS = { + "ph": "cmems_mod_glo_bgc_my_0.25deg_P1M-m", + "phyc": "cmems_mod_glo_bgc_my_0.25deg_P1M-m", +} +MONTHLY_BGC_REANALYSIS_INTERIM_IDS = { + "ph": "cmems_mod_glo_bgc_myint_0.25deg_P1M-m", + "phyc": "cmems_mod_glo_bgc_myint_0.25deg_P1M-m", +} + +# variables used in VirtualShip which are physical or biogeochemical variables, respectively +COPERNICUSMARINE_PHYS_VARIABLES = ["uo", "vo", "so", "thetao"] +COPERNICUSMARINE_BGC_VARIABLES = ["o2", "chl", "no3", "po4", "ph", "phyc", "nppv"] + + +def _select_product_id( + physical: bool, + schedule_start, + schedule_end, + username: str | None = None, + password: str | None = None, + variable: str | None = None, +) -> str: + """Determine which copernicus product id should be selected (reanalysis, reanalysis-interim, analysis & forecast), for prescribed schedule and physical vs. BGC.""" + key = "phys" if physical else "bgc" + selected_id = None + + for period, pid in PRODUCT_IDS[key].items(): + # for BGC analysis, set pid per variable + if key == "bgc" and period == "analysis": + if variable is None or variable not in BGC_ANALYSIS_IDS: + continue + pid = BGC_ANALYSIS_IDS[variable] + # for BGC reanalysis, check if requires monthly product + if ( + key == "bgc" + and period == "reanalysis" + and variable in MONTHLY_BGC_REANALYSIS_IDS + ): + monthly_pid = MONTHLY_BGC_REANALYSIS_IDS[variable] + ds_monthly = copernicusmarine.open_dataset( + monthly_pid, + username=username, + password=password, + ) + time_end_monthly = ds_monthly["time"][-1].values + if np.datetime64(schedule_end) <= time_end_monthly: + pid = monthly_pid + # for BGC reanalysis_interim, check if requires monthly product + if ( + key == "bgc" + and period == "reanalysis_interim" + and variable in MONTHLY_BGC_REANALYSIS_INTERIM_IDS + ): + monthly_pid = MONTHLY_BGC_REANALYSIS_INTERIM_IDS[variable] + ds_monthly = copernicusmarine.open_dataset( + monthly_pid, username=username, password=password + ) + time_end_monthly = ds_monthly["time"][-1].values + if np.datetime64(schedule_end) <= time_end_monthly: + pid = monthly_pid + if pid is None: + continue + ds = copernicusmarine.open_dataset(pid, username=username, password=password) + time_end = ds["time"][-1].values + if np.datetime64(schedule_end) <= time_end: + selected_id = pid + break + + if selected_id is None: + raise CopernicusCatalogueError( + "No suitable product found in the Copernicus Marine Catalogue for the scheduled time and variable." + ) + + if _start_end_in_product_timerange( + selected_id, schedule_start, schedule_end, username, password + ): + return selected_id + else: + return ( + PRODUCT_IDS["phys"]["analysis"] if physical else BGC_ANALYSIS_IDS[variable] + ) + + +def _start_end_in_product_timerange( + selected_id, schedule_start, schedule_end, username, password +): + ds_selected = copernicusmarine.open_dataset( + selected_id, username=username, password=password + ) + time_values = ds_selected["time"].values + import numpy as np + + time_min, time_max = np.min(time_values), np.max(time_values) + return ( + np.datetime64(schedule_start) >= time_min + and np.datetime64(schedule_end) <= time_max + ) + + +def _get_bathy_data(space_time_region, latlon_buffer: float | None = None) -> FieldSet: + """Bathymetry data 'streamed' directly from Copernicus Marine.""" + ds_bathymetry = copernicusmarine.open_dataset( + dataset_id="cmems_mod_glo_phy_my_0.083deg_static", + minimum_longitude=space_time_region.spatial_range.minimum_longitude + - (latlon_buffer if latlon_buffer is not None else 0), + maximum_longitude=space_time_region.spatial_range.maximum_longitude + + (latlon_buffer if latlon_buffer is not None else 0), + minimum_latitude=space_time_region.spatial_range.minimum_latitude + - (latlon_buffer if latlon_buffer is not None else 0), + maximum_latitude=space_time_region.spatial_range.maximum_latitude + + (latlon_buffer if latlon_buffer is not None else 0), + variables=["deptho"], + start_datetime=space_time_region.time_range.start_time, + end_datetime=space_time_region.time_range.end_time, + coordinates_selection_method="outside", + ) + bathymetry_variables = {"bathymetry": "deptho"} + bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"} + return FieldSet.from_xarray_dataset( + ds_bathymetry, bathymetry_variables, bathymetry_dimensions + ) + + +def expedition_cost(schedule_results: ScheduleOk, time_past: timedelta) -> float: + """ + Calculate the cost of the expedition in US$. + + :param schedule_results: Results from schedule simulation. + :param time_past: Time the expedition took. + :returns: The calculated cost of the expedition in US$. + """ + SHIP_COST_PER_DAY = 30000 + DRIFTER_DEPLOY_COST = 2500 + ARGO_DEPLOY_COST = 15000 + + ship_cost = SHIP_COST_PER_DAY / 24 * time_past.total_seconds() // 3600 + num_argos = len(schedule_results.measurements_to_simulate.argo_floats) + argo_cost = num_argos * ARGO_DEPLOY_COST + num_drifters = len(schedule_results.measurements_to_simulate.drifters) + drifter_cost = num_drifters * DRIFTER_DEPLOY_COST + + cost = ship_cost + argo_cost + drifter_cost + return cost diff --git a/tests/cli/test_fetch.py b/tests/cli/test_fetch.py index 9725601a..4c041dbe 100644 --- a/tests/cli/test_fetch.py +++ b/tests/cli/test_fetch.py @@ -17,8 +17,6 @@ get_existing_download, hash_model, hash_to_filename, - select_product_id, - start_end_in_product_timerange, ) from virtualship.models import Expedition from virtualship.utils import EXPEDITION, get_example_expedition @@ -98,31 +96,6 @@ def test_complete_download(tmp_path): assert_complete_download(tmp_path) -@pytest.mark.usefixtures("copernicus_no_download") -def test_select_product_id(expedition): - """Should return the physical reanalysis product id via the timings prescribed in the static schedule.yaml file.""" - result = select_product_id( - physical=True, - schedule_start=expedition.schedule.space_time_region.time_range.start_time, - schedule_end=expedition.schedule.space_time_region.time_range.end_time, - username="test", - password="test", - ) - assert result == "cmems_mod_glo_phy_my_0.083deg_P1D-m" - - -@pytest.mark.usefixtures("copernicus_no_download") -def test_start_end_in_product_timerange(expedition): - """Should return True for valid range ass determined by the static schedule.yaml file.""" - assert start_end_in_product_timerange( - selected_id="cmems_mod_glo_phy_my_0.083deg_P1D-m", - schedule_start=expedition.schedule.space_time_region.time_range.start_time, - schedule_end=expedition.schedule.space_time_region.time_range.end_time, - username="test", - password="test", - ) - - def test_assert_complete_download_complete(tmp_path): # Setup DownloadMetadata(download_complete=True).to_yaml(tmp_path / DOWNLOAD_METADATA) diff --git a/tests/expedition/test_expedition.py b/tests/expedition/test_expedition.py index a4643e03..4bed35e7 100644 --- a/tests/expedition/test_expedition.py +++ b/tests/expedition/test_expedition.py @@ -4,9 +4,14 @@ import pyproj import pytest -from virtualship.errors import ConfigError, ScheduleError +from virtualship.errors import InstrumentsConfigError, ScheduleError from virtualship.expedition.do_expedition import _load_input_data -from virtualship.models import Expedition, Location, Schedule, Waypoint +from virtualship.models import ( + Expedition, + Location, + Schedule, + Waypoint, +) from virtualship.utils import EXPEDITION, _get_expedition, get_example_expedition projection = pyproj.Geod(ellps="WGS84") @@ -56,6 +61,7 @@ def test_verify_schedule() -> None: def test_get_instruments() -> None: + get_expedition = _get_expedition(expedition_dir) schedule = Schedule( waypoints=[ Waypoint(location=Location(0, 0), instrument=["CTD"]), @@ -63,12 +69,21 @@ def test_get_instruments() -> None: Waypoint(location=Location(1, 0), instrument=["CTD"]), ] ) - - assert set(instrument.name for instrument in schedule.get_instruments()) == { - "CTD", - "XBT", - "ARGO_FLOAT", - } + expedition = Expedition( + schedule=schedule, + instruments_config=get_expedition.instruments_config, + ship_config=get_expedition.ship_config, + ) + assert ( + set(instrument.name for instrument in expedition.get_instruments()) + == { + "CTD", + "UNDERWATER_ST", # not added above but underway instruments are auto present from instruments_config in expedition_dir/expedition.yaml + "ADCP", # as above + "ARGO_FLOAT", + "XBT", + } + ) @pytest.mark.parametrize( @@ -165,15 +180,15 @@ def test_verify_schedule_errors( @pytest.fixture -def schedule(tmp_file): +def expedition(tmp_file): with open(tmp_file, "w") as file: file.write(get_example_expedition()) - return Expedition.from_yaml(tmp_file).schedule + return Expedition.from_yaml(tmp_file) @pytest.fixture -def schedule_no_xbt(schedule): - for waypoint in schedule.waypoints: +def expedition_no_xbt(expedition): + for waypoint in expedition.schedule.waypoints: if waypoint.instrument and any( instrument.name == "XBT" for instrument in waypoint.instrument ): @@ -183,54 +198,57 @@ def schedule_no_xbt(schedule): if instrument.name != "XBT" ] - return schedule + return expedition @pytest.fixture -def instruments_config(tmp_file): - with open(tmp_file, "w") as file: - file.write(get_example_expedition()) - return Expedition.from_yaml(tmp_file).instruments_config +def instruments_config_no_xbt(expedition): + delattr(expedition.instruments_config, "xbt_config") + return expedition.instruments_config @pytest.fixture -def instruments_config_no_xbt(instruments_config): - delattr(instruments_config, "xbt_config") - return instruments_config +def instruments_config_no_ctd(expedition): + delattr(expedition.instruments_config, "ctd_config") + return expedition.instruments_config @pytest.fixture -def instruments_config_no_ctd(instruments_config): - delattr(instruments_config, "ctd_config") - return instruments_config +def instruments_config_no_ctd_bgc(expedition): + delattr(expedition.instruments_config, "ctd_bgc_config") + return expedition.instruments_config @pytest.fixture -def instruments_config_no_ctd_bgc(instruments_config): - delattr(instruments_config, "ctd_bgc_config") - return instruments_config +def instruments_config_no_argo_float(expedition): + delattr(expedition.instruments_config, "argo_float_config") + return expedition.instruments_config @pytest.fixture -def instruments_config_no_argo_float(instruments_config): - delattr(instruments_config, "argo_float_config") - return instruments_config +def instruments_config_no_drifter(expedition): + delattr(expedition.instruments_config, "drifter_config") + return expedition.instruments_config @pytest.fixture -def instruments_config_no_drifter(instruments_config): - delattr(instruments_config, "drifter_config") - return instruments_config +def instruments_config_no_adcp(expedition): + delattr(expedition.instruments_config, "adcp_config") + return expedition.instruments_config -def test_verify_instruments_config(instruments_config, schedule) -> None: - instruments_config.verify(schedule) +@pytest.fixture +def instruments_config_no_underwater_st(expedition): + delattr(expedition.instruments_config, "ship_underwater_st_config") + return expedition.instruments_config -def test_verify_instruments_config_no_instrument( - instruments_config, schedule_no_xbt -) -> None: - instruments_config.verify(schedule_no_xbt) +def test_verify_instruments_config(expedition) -> None: + expedition.instruments_config.verify(expedition) + + +def test_verify_instruments_config_no_instrument(expedition, expedition_no_xbt) -> None: + expedition.instruments_config.verify(expedition_no_xbt) @pytest.mark.parametrize( @@ -238,40 +256,52 @@ def test_verify_instruments_config_no_instrument( [ pytest.param( "instruments_config_no_xbt", - ConfigError, - "Schedule includes instrument 'XBT', but instruments_config does not provide configuration for it.", - id="ShipConfigNoXBT", + InstrumentsConfigError, + "Expedition includes instrument 'XBT', but instruments_config does not provide configuration for it.", + id="InstrumentsConfigNoXBT", ), pytest.param( "instruments_config_no_ctd", - ConfigError, - "Schedule includes instrument 'CTD', but instruments_config does not provide configuration for it.", - id="ShipConfigNoCTD", + InstrumentsConfigError, + "Expedition includes instrument 'CTD', but instruments_config does not provide configuration for it.", + id="InstrumentsConfigNoCTD", ), pytest.param( "instruments_config_no_ctd_bgc", - ConfigError, - "Schedule includes instrument 'CTD_BGC', but instruments_config does not provide configuration for it.", - id="ShipConfigNoCTD_BGC", + InstrumentsConfigError, + "Expedition includes instrument 'CTD_BGC', but instruments_config does not provide configuration for it.", + id="InstrumentsConfigNoCTD_BGC", ), pytest.param( "instruments_config_no_argo_float", - ConfigError, - "Schedule includes instrument 'ARGO_FLOAT', but instruments_config does not provide configuration for it.", - id="ShipConfigNoARGO_FLOAT", + InstrumentsConfigError, + "Expedition includes instrument 'ARGO_FLOAT', but instruments_config does not provide configuration for it.", + id="InstrumentsConfigNoARGO_FLOAT", ), pytest.param( "instruments_config_no_drifter", - ConfigError, - "Schedule includes instrument 'DRIFTER', but instruments_config does not provide configuration for it.", - id="ShipConfigNoDRIFTER", + InstrumentsConfigError, + "Expedition includes instrument 'DRIFTER', but instruments_config does not provide configuration for it.", + id="InstrumentsConfigNoDRIFTER", + ), + pytest.param( + "instruments_config_no_adcp", + InstrumentsConfigError, + r"Underway instrument config attribute\(s\) are missing from YAML\. Must be Config object or None\.", + id="InstrumentsConfigNoADCP", + ), + pytest.param( + "instruments_config_no_underwater_st", + InstrumentsConfigError, + r"Underway instrument config attribute\(s\) are missing from YAML\. Must be Config object or None\.", + id="InstrumentsConfigNoUNDERWATER_ST", ), ], ) def test_verify_instruments_config_errors( - request, schedule, instruments_config_fixture, error, match + request, expedition, instruments_config_fixture, error, match ) -> None: instruments_config = request.getfixturevalue(instruments_config_fixture) with pytest.raises(error, match=match): - instruments_config.verify(schedule) + instruments_config.verify(expedition) diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py new file mode 100644 index 00000000..0a58efcd --- /dev/null +++ b/tests/instruments/test_base.py @@ -0,0 +1,150 @@ +import datetime +from pathlib import Path +from unittest.mock import MagicMock, patch + +import numpy as np +import pytest +import xarray as xr + +from virtualship.instruments.base import InputDataset, Instrument +from virtualship.instruments.types import InstrumentType +from virtualship.models.space_time_region import ( + SpaceTimeRegion, + SpatialRange, + TimeRange, +) +from virtualship.utils import get_input_dataset_class + +# test dataclass, particle class, kernels, etc. are defined for each instrument + + +def test_all_instruments_have_input_class(): + for instrument in InstrumentType: + input_class = get_input_dataset_class(instrument) + assert input_class is not None, f"No input_class for {instrument}" + + +# test InputDataset class + + +class DummyInputDataset(InputDataset): + """A minimal InputDataset subclass for testing purposes.""" + + def get_datasets_dict(self): + """Return a dummy datasets dict for testing.""" + return { + "dummy": { + "physical": True, + "variables": ["var1"], + "output_filename": "dummy.nc", + } + } + + +@pytest.fixture +def dummy_space_time_region(): + spatial_range = SpatialRange( + minimum_longitude=0, + maximum_longitude=1, + minimum_latitude=0, + maximum_latitude=1, + minimum_depth=0, + maximum_depth=10, + ) + base_time = datetime.datetime.strptime("1950-01-01", "%Y-%m-%d") + time_range = TimeRange( + start_time=base_time, + end_time=base_time + datetime.timedelta(hours=1), + ) + return SpaceTimeRegion( + spatial_range=spatial_range, + time_range=time_range, + ) + + +def test_dummyinputdataset_initialization(dummy_space_time_region): + ds = DummyInputDataset( + name="test", + latlon_buffer=0.5, + datetime_buffer=1, + min_depth=0, + max_depth=10, + data_dir=".", + credentials={"username": "u", "password": "p"}, + space_time_region=dummy_space_time_region, + ) + assert ds.name == "test" + assert ds.latlon_buffer == 0.5 + assert ds.datetime_buffer == 1 + assert ds.min_depth == 0 + assert ds.max_depth == 10 + assert ds.data_dir == "." + assert ds.credentials["username"] == "u" + + +@patch("virtualship.instruments.base.copernicusmarine.open_dataset") +@patch("virtualship.instruments.base.copernicusmarine.subset") +def test_download_data_calls_subset( + mock_subset, mock_open_dataset, dummy_space_time_region +): + """Test that download_data calls the subset function correctly, will also test Copernicus Marine product id search logic.""" + mock_open_dataset.return_value = xr.Dataset( + { + "time": ( + "time", + [ + np.datetime64("1993-01-01T00:00:00"), + np.datetime64("2023-01-01T01:00:00"), + ], + ) + } + ) + ds = DummyInputDataset( + name="test", + latlon_buffer=0.5, + datetime_buffer=1, + min_depth=0, + max_depth=10, + data_dir=".", + credentials={"username": "u", "password": "p"}, + space_time_region=dummy_space_time_region, + ) + ds.download_data() + assert mock_subset.called + + +# test Instrument class + + +class DummyInstrument(Instrument): + """Minimal concrete Instrument for testing.""" + + def simulate(self, data_dir, measurements, out_path): + """Dummy simulate implementation for test.""" + self.simulate_called = True + + +@patch("virtualship.instruments.base.FieldSet") +@patch("virtualship.instruments.base.get_existing_download") +@patch("virtualship.instruments.base.get_space_time_region_hash") +def test_load_input_data_calls(mock_hash, mock_get_download, mock_FieldSet): + """Test Instrument.load_input_data with mocks.""" + mock_hash.return_value = "hash" + mock_get_download.return_value = Path("/tmp/data") + mock_fieldset = MagicMock() + mock_FieldSet.from_netcdf.return_value = mock_fieldset + mock_fieldset.gridset.grids = [MagicMock(negate_depth=MagicMock())] + mock_fieldset.__getitem__.side_effect = lambda k: MagicMock() + dummy = DummyInstrument( + name="test", + expedition=MagicMock(schedule=MagicMock(space_time_region=MagicMock())), + directory="/tmp", + filenames={"A": "a.nc"}, + variables={"A": "a"}, + add_bathymetry=False, + allow_time_extrapolation=False, + verbose_progress=False, + ) + fieldset = dummy.load_input_data() + assert mock_FieldSet.from_netcdf.called + assert fieldset == mock_fieldset diff --git a/tests/test_utils.py b/tests/test_utils.py index 0dcebd79..bb8208f6 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,4 +1,4 @@ -from virtualship.models import Expedition +from virtualship.models.expedition import Expedition from virtualship.utils import get_example_expedition @@ -12,3 +12,19 @@ def test_valid_example_expedition(tmp_path): file.write(get_example_expedition()) Expedition.from_yaml(path) + + +def test_instrument_registry_updates(): + from virtualship import utils + + class DummyInputDataset: + pass + + class DummyInstrument: + pass + + utils.register_input_dataset("DUMMY_TYPE")(DummyInputDataset) + utils.register_instrument("DUMMY_TYPE")(DummyInstrument) + + assert utils.INPUT_DATASET_MAP["DUMMY_TYPE"] is DummyInputDataset + assert utils.INSTRUMENT_CLASS_MAP["DUMMY_TYPE"] is DummyInstrument