diff --git a/tests/sailship_config.json b/tests/sailship_config.json deleted file mode 100644 index 1f3805dc..00000000 --- a/tests/sailship_config.json +++ /dev/null @@ -1,40 +0,0 @@ -{ - "region_of_interest": { - "North": 64, - "East": -23, - "South": 59, - "West": -43 - }, - "requested_ship_time": { - "start": "2022-01-01T00:00:00", - "end": "2022-01-2T00:00:00" - }, - "route_coordinates": [ - [-23.071289, 63.743631], [-23.081289, 63.743631], [-23.091289, 63.743631] - ], - "underway_data": true, - "ADCP_data": false, - "ADCP_settings": { - "max_depth": -1000, - "bin_size_m": 24 - }, - "CTD_locations": [ - [-23.081289, 63.743631] - ], - "CTD_settings": { - "max_depth": "max" - }, - "drifter_deploylocations": [ - [-23.081289, 63.743631] - ], - "argo_deploylocations": [ - [-23.081289, 63.743631] - ], - "argo_characteristics": { - "driftdepth": -1000, - "maxdepth": -2000, - "vertical_speed": -0.10, - "cycle_days" : 10, - "drift_days": 9 - } - } diff --git a/tests/test_sailship.py b/tests/test_sailship.py index 6cbe6b0d..2f3dd974 100644 --- a/tests/test_sailship.py +++ b/tests/test_sailship.py @@ -5,8 +5,13 @@ import numpy as np from parcels import Field, FieldSet +from virtual_ship import Location from virtual_ship.sailship import sailship -from virtual_ship.virtual_ship_configuration import VirtualShipConfiguration +from virtual_ship.virtual_ship_configuration import ( + ADCPConfig, + ArgoFloatConfig, + VirtualShipConfig, +) def _make_ctd_fieldset(base_time: datetime) -> FieldSet: @@ -77,13 +82,37 @@ def test_sailship() -> None: }, ) - config = VirtualShipConfiguration( - "sailship_config.json", + argo_float_config = ArgoFloatConfig( + fieldset=argo_float_fieldset, + max_depth=-2000, + drift_depth=-1000, + vertical_speed=-0.10, + cycle_days=10, + drift_days=9, + ) + + adcp_config = ADCPConfig(max_depth=-1000, bin_size_m=24) + + config = VirtualShipConfig( + start_time=datetime.datetime.strptime( + "2022-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S" + ), + route_coordinates=[ + Location(latitude=-23.071289, longitude=63.743631), + Location(latitude=-23.081289, longitude=63.743631), + Location(latitude=-23.191289, longitude=63.743631), + ], adcp_fieldset=adcp_fieldset, ship_underwater_st_fieldset=ship_underwater_st_fieldset, ctd_fieldset=ctd_fieldset, drifter_fieldset=drifter_fieldset, - argo_float_fieldset=argo_float_fieldset, + argo_float_deploy_locations=[ + Location(latitude=-23.081289, longitude=63.743631) + ], + drifter_deploy_locations=[Location(latitude=-23.081289, longitude=63.743631)], + ctd_deploy_locations=[Location(latitude=-23.081289, longitude=63.743631)], + argo_float_config=argo_float_config, + adcp_config=adcp_config, ) sailship(config) diff --git a/virtual_ship/costs.py b/virtual_ship/costs.py index 9ba1d970..459dba67 100644 --- a/virtual_ship/costs.py +++ b/virtual_ship/costs.py @@ -2,10 +2,10 @@ from datetime import timedelta -from .virtual_ship_configuration import VirtualShipConfiguration +from .virtual_ship_configuration import VirtualShipConfig -def costs(config: VirtualShipConfiguration, total_time: timedelta): +def costs(config: VirtualShipConfig, total_time: timedelta): """ Calculate the cost of the virtual ship (in US$). @@ -18,8 +18,8 @@ def costs(config: VirtualShipConfiguration, total_time: timedelta): argo_deploy_cost = 15000 ship_cost = ship_cost_per_day / 24 * total_time.total_seconds() // 3600 - argo_cost = len(config.argo_deploylocations) * argo_deploy_cost - drifter_cost = len(config.drifter_deploylocations) * drifter_deploy_cost + argo_cost = len(config.argo_float_deploy_locations) * argo_deploy_cost + drifter_cost = len(config.drifter_deploy_locations) * drifter_deploy_cost cost = ship_cost + argo_cost + drifter_cost return cost diff --git a/virtual_ship/sailship.py b/virtual_ship/sailship.py index b0eb3d6c..9a02c828 100644 --- a/virtual_ship/sailship.py +++ b/virtual_ship/sailship.py @@ -1,12 +1,10 @@ """sailship function.""" -import datetime import os from datetime import timedelta import numpy as np import pyproj -from shapely.geometry import Point, Polygon from .costs import costs from .instruments.adcp import simulate_adcp @@ -17,49 +15,34 @@ from .location import Location from .postprocess import postprocess from .spacetime import Spacetime -from .virtual_ship_configuration import VirtualShipConfiguration +from .virtual_ship_configuration import VirtualShipConfig -def sailship(config: VirtualShipConfiguration): +def sailship(config: VirtualShipConfig): """ Use parcels to simulate the ship, take ctd_instruments and measure ADCP and underwaydata. :param config: The cruise configuration. """ - # TODO this will be in the config later, but for now we don't change the config structure - # from here ----- - argo_locations_list = [ - Location(latitude=argo[1], longitude=argo[0]) - for argo in config.argo_deploylocations - ] - argo_locations = set(argo_locations_list) - if len(argo_locations) != len(argo_locations_list): + config.verify() + + # combine identical instrument deploy location + argo_locations = set(config.argo_float_deploy_locations) + if len(argo_locations) != len(config.argo_float_deploy_locations): print( "WARN: Some argo float deployment locations are identical and have been combined." ) - drifter_locations_list = [ - Location(latitude=drifter[1], longitude=drifter[0]) - for drifter in config.drifter_deploylocations - ] - drifter_locations = set(drifter_locations_list) - if len(drifter_locations) != len(drifter_locations_list): + drifter_locations = set(config.drifter_deploy_locations) + if len(drifter_locations) != len(config.drifter_deploy_locations): print( "WARN: Some drifter deployment locations are identical and have been combined." ) - ctd_locations_list = [ - Location(latitude=ctd[1], longitude=ctd[0]) for ctd in config.CTD_locations - ] - ctd_locations = set(ctd_locations_list) - if len(ctd_locations) != len(ctd_locations_list): + ctd_locations = set(config.ctd_deploy_locations) + if len(drifter_locations) != len(config.ctd_deploy_locations): print("WARN: Some CTD locations are identical and have been combined.") - start_time = datetime.datetime.strptime( - config.requested_ship_time["start"], "%Y-%m-%dT%H:%M:%S" - ) - # until here ---- - # get discrete points along the ships route were sampling and deployments will be performed route_dt = timedelta(minutes=5) route_points = shiproute(config=config, dt=route_dt) @@ -139,12 +122,12 @@ def sailship(config: VirtualShipConfiguration): spacetime=Spacetime( location=route_point, time=time_past.total_seconds() ), - min_depth=-config.argo_float_fieldset.U.depth[0], - max_depth=config.argo_characteristics["maxdepth"], - drift_depth=config.argo_characteristics["driftdepth"], - vertical_speed=config.argo_characteristics["vertical_speed"], - cycle_days=config.argo_characteristics["cycle_days"], - drift_days=config.argo_characteristics["drift_days"], + min_depth=-config.argo_float_config.fieldset.U.depth[0], + max_depth=config.argo_float_config.max_depth, + drift_depth=config.argo_float_config.drift_depth, + vertical_speed=config.argo_float_config.vertical_speed, + cycle_days=config.argo_float_config.cycle_days, + drift_days=config.argo_float_config.drift_days, ) ) argo_locations_visited = argo_locations_visited.union(argos_here) @@ -168,7 +151,7 @@ def sailship(config: VirtualShipConfiguration): CTD( spacetime=Spacetime( location=route_point, - time=start_time + time_past, + time=config.start_time + time_past, ), min_depth=config.ctd_fieldset.U.depth[0], max_depth=config.ctd_fieldset.U.depth[-1], @@ -212,10 +195,9 @@ def sailship(config: VirtualShipConfiguration): simulate_adcp( fieldset=config.adcp_fieldset, out_path=os.path.join("results", "adcp.zarr"), - max_depth=config.ADCP_settings["max_depth"], + max_depth=config.adcp_config.max_depth, min_depth=-5, - num_bins=(-5 - config.ADCP_settings["max_depth"]) - // config.ADCP_settings["bin_size_m"], + num_bins=(-5 - config.adcp_config.max_depth) // config.adcp_config.bin_size_m, sample_points=adcps, ) @@ -241,7 +223,7 @@ def sailship(config: VirtualShipConfiguration): simulate_argo_floats( out_path=os.path.join("results", "argo_floats.zarr"), argo_floats=argo_floats, - fieldset=config.argo_float_fieldset, + fieldset=config.argo_float_config.fieldset, outputdt=timedelta(minutes=5), endtime=None, ) @@ -256,7 +238,7 @@ def sailship(config: VirtualShipConfiguration): print(f"This cruise took {time_past} and would have cost {cost:,.0f} euros.") -def shiproute(config: VirtualShipConfiguration, dt: timedelta) -> list[Location]: +def shiproute(config: VirtualShipConfig, dt: timedelta) -> list[Location]: """ Take in route coordinates and return lat and lon points within region of interest to sample. @@ -264,57 +246,44 @@ def shiproute(config: VirtualShipConfiguration, dt: timedelta) -> list[Location] :param dt: Sailing time between each discrete route point. :returns: lat and lon points within region of interest to sample. """ - # Initialize lists to store intermediate points - lons = [] - lats = [] - - # Loop over station coordinates and calculate intermediate points along great circle path - for i in range(len(config.route_coordinates) - 1): - startlong = config.route_coordinates[i][0] - startlat = config.route_coordinates[i][1] - endlong = config.route_coordinates[i + 1][0] - endlat = config.route_coordinates[i + 1][1] - - # calculate line string along path with segments every 5 min for ADCP measurements - # current cruise speed is 10knots = 5.14 m/s * 60*5 = 1542 m every 5 min - # Realistic time between measurements is 2 min on Pelagia according to Floran - cruise_speed = 5.14 - geod = pyproj.Geod(ellps="WGS84") - azimuth1, azimuth2, distance = geod.inv(startlong, startlat, endlong, endlat) - if distance > (cruise_speed * dt.total_seconds()): - r = geod.inv_intermediate( - startlong, - startlat, - endlong, - endlat, - del_s=1545, - initial_idx=0, - return_back_azimuth=False, - ) - lons = np.append(lons, r.lons) # stored as a list of arrays - lats = np.append(lats, r.lats) - else: - lons = np.append(lons, endlong) - lats = np.append(lats, endlat) - - # initial_idx will add begin point to each list (but not end point to avoid doubling) so add final endpoint manually - lons = np.append(np.hstack(lons), endlong) - lats = np.append(np.hstack(lats), endlat) - - # check if input sample locations are within data availability area, only save if so - north = config.region_of_interest["North"] - east = config.region_of_interest["East"] - south = config.region_of_interest["South"] - west = config.region_of_interest["West"] - poly = Polygon([(west, north), (west, south), (east, south), (east, north)]) - sample_lons = [] - sample_lats = [] - for i in range(len(lons)): - if poly.contains(Point(lons[i], lats[i])): - sample_lons.append(lons[i]) - sample_lats.append(lats[i]) - points = [ - Location(latitude=lat, longitude=lon) - for lat, lon in zip(sample_lats, sample_lons, strict=True) - ] - return points + CRUISE_SPEED = 5.14 + + # discrete points the ship will pass + sample_points: list[Location] = [] + + # projection used to get discrete locations + geod = pyproj.Geod(ellps="WGS84") + + # loop over station coordinates and calculate intermediate points along great circle path + for startloc, endloc in zip(config.route_coordinates, config.route_coordinates[1:]): + # iterate over each coordinate and the next coordinate + # last coordinate has no next coordinate and is skipped + + # get locations between start and end, seperate by 5 minutes of cruising + # excludes final point, but this is added explicitly after this loop + int_points = geod.inv_intermediate( + startloc.lon, + startloc.lat, + endloc.lon, + endloc.lat, + del_s=CRUISE_SPEED * dt.total_seconds(), + initial_idx=0, + return_back_azimuth=False, + ) + + sample_points.extend( + [ + Location(latitude=lat, longitude=lon) + for lat, lon in zip(int_points.lats, int_points.lons, strict=True) + ] + ) + + # explitly include final point which is not added by the previous loop + sample_points.append( + Location( + latitude=config.route_coordinates[-1].lat, + longitude=config.route_coordinates[-1].lon, + ) + ) + + return sample_points diff --git a/virtual_ship/virtual_ship_configuration.py b/virtual_ship/virtual_ship_configuration.py index ca2b8547..8d48d4c0 100644 --- a/virtual_ship/virtual_ship_configuration.py +++ b/virtual_ship/virtual_ship_configuration.py @@ -1,173 +1,156 @@ -"""VirtualShipConfiguration class.""" +"""VirtualShipConfig class.""" -import datetime -import json -from datetime import timedelta +from dataclasses import dataclass +from datetime import datetime +import numpy as np from parcels import FieldSet -from shapely.geometry import Point, Polygon +from .location import Location -class VirtualShipConfiguration: - """Configuration of the virtual ship, initialized from a json file.""" + +@dataclass +class ArgoFloatConfig: + """Configuration for argos floats.""" + + fieldset: FieldSet + max_depth: float + drift_depth: float + vertical_speed: float + cycle_days: float + drift_days: float + + +@dataclass +class ADCPConfig: + """Configuration for ADCP instrument.""" + + max_depth: float + bin_size_m: int + + +@dataclass +class VirtualShipConfig: + """Configuration of the virtual ship.""" + + start_time: datetime + route_coordinates: list[Location] adcp_fieldset: FieldSet ship_underwater_st_fieldset: FieldSet - argo_float_fieldset: FieldSet - drifter_fieldset: FieldSet ctd_fieldset: FieldSet drifter_fieldset: FieldSet - argo_float_fieldset: FieldSet - - def __init__( - self, - json_file, - adcp_fieldset: FieldSet, - ship_underwater_st_fieldset: FieldSet, - ctd_fieldset: FieldSet, - drifter_fieldset: FieldSet, - argo_float_fieldset: FieldSet, - ): + + argo_float_deploy_locations: list[Location] + drifter_deploy_locations: list[Location] + ctd_deploy_locations: list[Location] + + argo_float_config: ArgoFloatConfig + adcp_config: ADCPConfig + + def verify(self) -> None: """ - Initialize this object. - - :param json_file: Path to the JSON file to init from. - :param adcp_fieldset: Fieldset for ADCP measurements. - :param ship_underwater_st_fieldset: Fieldset for ship salinity temperature measurements. - :param ctd_fieldset: Fieldset for CTD measurements. - :param drifter_fieldset: Fieldset for CTD measurements. - :param argo_float_fieldset: Fieldset for argo float measurements. - :raises ValueError: If JSON file not valid. + Verify this configuration is valid. + + :raises ValueError: If not valid. """ - self.adcp_fieldset = adcp_fieldset - self.ship_underwater_st_fieldset = ship_underwater_st_fieldset - self.ctd_fieldset = ctd_fieldset - self.drifter_fieldset = drifter_fieldset - self.argo_float_fieldset = argo_float_fieldset - - with open(json_file, "r") as file: - json_input = json.loads(file.read()) - for key in json_input: - setattr(self, key, json_input[key]) - - # Create a polygon from the region of interest to check coordinates - north = self.region_of_interest["North"] - east = self.region_of_interest["East"] - south = self.region_of_interest["South"] - west = self.region_of_interest["West"] - poly = Polygon([(west, north), (west, south), (east, south), (east, north)]) - # validate input, raise ValueErrors if invalid - if ( - self.region_of_interest["North"] > 90 - or self.region_of_interest["South"] < -90 - or self.region_of_interest["East"] > 180 - or self.region_of_interest["West"] < -180 + if len(self.route_coordinates) < 2: + raise ValueError("Route needs to consist of at least locations.") + + if not all( + [self._is_valid_location(coord) for coord in self.route_coordinates] ): - raise ValueError("Invalid coordinates in region of interest") - if self.requested_ship_time["start"] > self.requested_ship_time["end"]: - raise ValueError("Start time should be before end time") - if datetime.datetime.strptime( - self.requested_ship_time["end"], "%Y-%m-%dT%H:%M:%S" - ) > datetime.datetime.now() + timedelta(days=2): - raise ValueError("End time cannot be more then 2 days into the future") - if datetime.datetime.strptime( - self.requested_ship_time["end"], "%Y-%m-%dT%H:%M:%S" - ) - datetime.datetime.strptime( - self.requested_ship_time["start"], "%Y-%m-%dT%H:%M:%S" - ) > timedelta( - days=21 + raise ValueError("Invalid coordinates in route.") + + if not all( + [ + self._is_valid_location(coord) + for coord in self.argo_float_deploy_locations + ] ): - raise ValueError("The time period is too long, maximum is 21 days") - if len(self.route_coordinates) < 2: - raise ValueError( - "Route needs to consist of at least 2 longitude-latitude coordinate sets" - ) - for coord in self.route_coordinates: - if coord[1] > 90 or coord[1] < -90 or coord[0] > 180 or coord[0] < -180: - raise ValueError("Invalid coordinates in route") - # if not poly.contains(Point(coord)): - # raise ValueError("Route coordinates need to be within the region of interest") - if not len(self.CTD_locations) == 0: - for coord in self.CTD_locations: - if coord not in self.route_coordinates: - raise ValueError("CTD coordinates should be on the route") - if coord[1] > 90 or coord[1] < -90 or coord[0] > 180 or coord[0] < -180: - raise ValueError("Invalid coordinates in route") - if not poly.contains(Point(coord)): - raise ValueError( - "CTD coordinates need to be within the region of interest" - ) - if type(self.CTD_settings["max_depth"]) is not int: - if self.CTD_settings["max_depth"] != "max": - raise ValueError( - 'Specify "max" for maximum depth or a negative integer for max_depth in CTD_settings' + raise ValueError("Argo float deploy locations are not valid coordinates.") + + if not all( + [ + any( + [ + np.isclose(deploy.lat, coord.lat) + and np.isclose(deploy.lon, coord.lon) + for coord in self.route_coordinates + ] ) - if type(self.CTD_settings["max_depth"]) is int: - if self.CTD_settings["max_depth"] > 0: - raise ValueError("Invalid depth for CTD") - if len(self.drifter_deploylocations) > 30: - raise ValueError("Too many drifter deployment locations, maximum is 30") - if not len(self.drifter_deploylocations) == 0: - for coord in self.drifter_deploylocations: - if coord not in self.route_coordinates: - raise ValueError("Drifter coordinates should be on the route") - if coord[1] > 90 or coord[1] < -90 or coord[0] > 180 or coord[0] < -180: - raise ValueError("Invalid coordinates in route") - if not poly.contains(Point(coord)): - raise ValueError( - "Drifter coordinates need to be within the region of interest" - ) - if len(self.argo_deploylocations) > 30: - raise ValueError("Too many argo deployment locations, maximum is 30") - if not len(self.argo_deploylocations) == 0: - for coord in self.argo_deploylocations: - if coord not in self.route_coordinates: - raise ValueError("argo coordinates should be on the route") - if coord[1] > 90 or coord[1] < -90 or coord[0] > 180 or coord[0] < -180: - raise ValueError("Invalid coordinates in route") - if not poly.contains(Point(coord)): - raise ValueError( - "argo coordinates need to be within the region of interest" - ) - if not isinstance(self.underway_data, bool): - raise ValueError("Underway data needs to be true or false") - if not isinstance(self.ADCP_data, bool): - raise ValueError("ADCP data needs to be true or false") - if ( - self.ADCP_settings["bin_size_m"] < 0 - or self.ADCP_settings["bin_size_m"] > 24 - ): - raise ValueError("Invalid bin size for ADCP") - if self.ADCP_settings["max_depth"] > 0: - raise ValueError("Invalid depth for ADCP") - if ( - self.argo_characteristics["driftdepth"] > 0 - or self.argo_characteristics["driftdepth"] < -5727 + for deploy in self.argo_float_deploy_locations + ] ): raise ValueError( - "Specify negative depth. Max drift depth for argo is -5727 m due to data availability" + "Argo float deploy locations are not exactly on route coordinates." ) - if ( - self.argo_characteristics["maxdepth"] > 0 - or self.argo_characteristics["maxdepth"] < -5727 + + if not all( + [self._is_valid_location(coord) for coord in self.drifter_deploy_locations] + ): + raise ValueError("Drifter deploy locations are not valid coordinates.") + + if not all( + [ + any( + [ + np.isclose(deploy.lat, coord.lat) + and np.isclose(deploy.lon, coord.lon) + for coord in self.route_coordinates + ] + ) + for deploy in self.drifter_deploy_locations + ] ): raise ValueError( - "Specify negative depth. Max depth for argo is -5727 m due to data availability" + "Drifter deploy locations are not exactly on route coordinates." ) - if type(self.argo_characteristics["vertical_speed"]) is not float: - raise ValueError("Specify vertical speed for argo with decimals in m/s") - if self.argo_characteristics["vertical_speed"] > 0: - self.argo_characteristics["vertical_speed"] = -self.argo_characteristics[ - "vertical_speed" + + if not all( + [self._is_valid_location(coord) for coord in self.ctd_deploy_locations] + ): + raise ValueError("CTD deploy locations are not valid coordinates.") + + if not all( + [ + any( + [ + np.isclose(deploy.lat, coord.lat) + and np.isclose(deploy.lon, coord.lon) + for coord in self.route_coordinates + ] + ) + for deploy in self.ctd_deploy_locations ] - if ( - abs(self.argo_characteristics["vertical_speed"]) < 0.06 - or abs(self.argo_characteristics["vertical_speed"]) > 0.12 ): raise ValueError( - "Specify a realistic speed for argo, i.e. between -0.06 and -0.12 m/s" + "CTD deploy locations are not exactly on route coordinates." ) - if self.argo_characteristics["cycle_days"] < 0: - raise ValueError("Specify a postitive number of cycle days for argo") - if self.argo_characteristics["drift_days"] < 0: - raise ValueError("Specify a postitive number of drift days for argo") + + if self.argo_float_config.max_depth > 0: + raise ValueError("Argo float max depth must be negative or zero.") + + if self.argo_float_config.drift_depth > 0: + raise ValueError("Argo float drift depth must be negative or zero.") + + if self.argo_float_config.vertical_speed >= 0: + raise ValueError("Argo float vertical speed must be negative.") + + if self.argo_float_config.cycle_days <= 0: + raise ValueError("Argo float cycle days must be larger than zero.") + + if self.argo_float_config.drift_days <= 0: + raise ValueError("Argo drift cycle days must be larger than zero.") + + if self.adcp_config.max_depth > 0: + raise ValueError("ADCP max depth must be negative.") + + @staticmethod + def _is_valid_location(location: Location) -> bool: + return ( + location.lat >= -90 + and location.lat <= 90 + and location.lon >= -180 + and location.lon <= 360 + )