From 4a9fb76a159ea9b358d8107cb93a95b32db8b581 Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Tue, 27 Sep 2022 00:29:17 -0700 Subject: [PATCH 1/7] LongTermAverageSource initial implementation --- .../data_sources/OceanCurrentField.py | 4 +- .../OceanCurrentSource/OceanCurrentSource.py | 38 +++++++++++++++++++ 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentField.py b/ocean_navigation_simulator/data_sources/OceanCurrentField.py index 2c5e1c5d..4e94d5d3 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentField.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentField.py @@ -3,7 +3,7 @@ from ocean_navigation_simulator.data_sources.DataField import DataField from ocean_navigation_simulator.data_sources.OceanCurrentSource.OceanCurrentSource import OceanCurrentSource from ocean_navigation_simulator.data_sources.OceanCurrentSource.OceanCurrentSource import \ - HindcastFileSource, HindcastOpendapSource, ForecastFileSource + HindcastFileSource, HindcastOpendapSource, ForecastFileSource, LongTermAverageSource import ocean_navigation_simulator.data_sources.OceanCurrentSource.AnalyticalOceanCurrents as AnalyticalSources @@ -47,6 +47,8 @@ def instantiate_source_from_dict(source_dict: Dict) -> OceanCurrentSource: return HindcastFileSource(source_dict) elif source_dict['source'] == 'forecast_files': return ForecastFileSource(source_dict) + elif source_dict['source'] == 'longterm_average': + return LongTermAverageSource(source_dict) elif source_dict['source'] == 'analytical': specific_analytical_current = getattr(AnalyticalSources, source_dict['source_settings']['name']) return specific_analytical_current(source_dict) diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py index 85094a14..da558dc7 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py @@ -2,6 +2,7 @@ import os from typing import List, AnyStr, Optional, Union import logging +from calendar import month import casadi as ca import dask.array.core @@ -309,6 +310,43 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean return OceanCurrentVector(u=self.u_curr_func(spatio_temporal_point.to_spatio_temporal_casadi_input()), v=self.v_curr_func(spatio_temporal_point.to_spatio_temporal_casadi_input())) +class LongTermAverageSource(OceanCurrentSource): # figure out inheritance + def __init__(self, source_config_dict: dict): + print("++++++") + print(source_config_dict) + self.forecast_data_source = ForecastFileSource(source_config_dict['forecast_dict']) + self.monthly_avg_data_source = HindcastFileSource(source_config_dict['monthly_avg_dict']) # defaults currents to normal + self.source_config_dict = source_config_dict + # self.t_0 = source_config_dict['t0'] # not sure what to do here + + def get_data_over_area(self, x_interval: List[float], y_interval: List[float], + t_interval: List[Union[datetime.datetime, int]], + spatial_resolution: Optional[float] = None, + temporal_resolution: Optional[float] = None) -> xr: + # Query as much forecast data as is possible + print("T_INTERVAL: " + str(t_interval)) + forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, spatial_resolution, temporal_resolution) + # last_time = forecast_dataframe["time"].resample(how="last") + # print("301") + # print(last_time) + remaining_t_interval = [get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]), t_interval[1]] # may not work + print("REMAINING T INTERVAL: " + str(remaining_t_interval)) + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, remaining_t_interval, spatial_resolution, temporal_resolution) + + return xr.concat([forecast_dataframe, monthly_average_dataframe], dim="time") + + def check_for_most_recent_fmrc_dataframe(self, time: datetime.datetime) -> int: + """Helper function to check update the self.OceanCurrent if a new forecast is available at + the specified input time. + Args: + time: datetime object + """ + return self.forecast_data_source.check_for_most_recent_fmrc_dataframe(time) + + # Not sure if I can just all this + def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> OceanCurrentVector: + """We overwrite it because we don't want that Forecast needs caching...""" + return self.forecast_data_source.get_data_at_point(spatio_temporal_point==spatio_temporal_point) # Helper functions across the OceanCurrentSource objects def get_file_dicts(folder: AnyStr, currents='normal') -> List[dict]: From 9d453579d8d3dac3e638478be83c52fe86d62adf Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Tue, 27 Sep 2022 00:42:11 -0700 Subject: [PATCH 2/7] testing implementation of LongTermAverageSource --- .../gulf_of_mexico_LongTermAverageSource.yaml | 40 +++++++++++++++++++ scripts/tutorial/controller/hj_planner.py | 2 +- 2 files changed, 41 insertions(+), 1 deletion(-) create mode 100644 config/arena/gulf_of_mexico_LongTermAverageSource.yaml diff --git a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml new file mode 100644 index 00000000..7dde0aaa --- /dev/null +++ b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml @@ -0,0 +1,40 @@ +casadi_cache_dict: + deg_around_x_t: 0.5 + time_around_x_t: 172800. # 3600 * 24 * 5 + +platform_dict: + battery_cap_in_wh: 400.0 + u_max_in_mps: 0.1 + motor_efficiency: 1.0 + solar_panel_size: 0.5 + solar_efficiency: 0.2 + drag_factor: 675.0 + dt_in_s: 600.0 + +use_geographic_coordinate_system: True + +spatial_boundary: null + +ocean_dict: + hindcast: + field: 'OceanCurrents' + source: 'hindcast_files' + source_settings: + folder: "data/hindcast_test" + source: "HYCOM" + type: "hindcast" + forecast: + field: 'OceanCurrents' + source: 'longterm_average' + source_settings: + folder: "data/forecast_test" + source: "HYCOM" + type: "hindcast" + +solar_dict: + hindcast: null + forecast: null + +seaweed_dict: + hindcast: null + forecast: null \ No newline at end of file diff --git a/scripts/tutorial/controller/hj_planner.py b/scripts/tutorial/controller/hj_planner.py index 241315c0..5007ff2e 100644 --- a/scripts/tutorial/controller/hj_planner.py +++ b/scripts/tutorial/controller/hj_planner.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt # Initialize the Arena (holds all data sources and the platform, everything except controller) -arena = ArenaFactory.create(scenario_name='gulf_of_mexico_HYCOM_hindcast_local') +arena = ArenaFactory.create(scenario_name='gulf_of_mexico_LongTermAverageSource') # we can also download the respective files directly to a temp folder, then t_interval needs to be set # % Specify Navigation Problem x_0 = PlatformState(lon=units.Distance(deg=-82.5), lat=units.Distance(deg=23.7), From 55650fd2afb1a9fd1a51565aecf4924816e8f8e3 Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Mon, 3 Oct 2022 00:53:52 -0700 Subject: [PATCH 3/7] updating format of yaml file and debugging --- config/arena/gulf_of_mexico_LongTermAverageSource.yaml | 9 ++++++++- .../OceanCurrentSource/OceanCurrentSource.py | 8 +++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml index 7dde0aaa..34897b8d 100644 --- a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml +++ b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml @@ -25,10 +25,17 @@ ocean_dict: type: "hindcast" forecast: field: 'OceanCurrents' - source: 'longterm_average' + source: 'forecast_files' source_settings: folder: "data/forecast_test" source: "HYCOM" + type: "forecast" + monthly_avg: + field: 'OceanCurrents' + source: 'longterm_average' + source_settings: + folder: "data/monthly_average" + source: "HYCOM" type: "hindcast" solar_dict: diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py index da558dc7..172bccc9 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py @@ -312,10 +312,8 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean class LongTermAverageSource(OceanCurrentSource): # figure out inheritance def __init__(self, source_config_dict: dict): - print("++++++") - print(source_config_dict) - self.forecast_data_source = ForecastFileSource(source_config_dict['forecast_dict']) - self.monthly_avg_data_source = HindcastFileSource(source_config_dict['monthly_avg_dict']) # defaults currents to normal + self.forecast_data_source = ForecastFileSource(source_config_dict['forecast']) + self.monthly_avg_data_source = HindcastFileSource(source_config_dict['monthly_avg']) # defaults currents to normal self.source_config_dict = source_config_dict # self.t_0 = source_config_dict['t0'] # not sure what to do here @@ -355,7 +353,7 @@ def get_file_dicts(folder: AnyStr, currents='normal') -> List[dict]: {'t_range': [, T], 'file': ,'y_range': [min_lat, max_lat], 'x_range': [min_lon, max_lon]} """ # get a list of files from the folder - files_list = [folder + f for f in os.listdir(folder) if + files_list = [folder + '/' + f for f in os.listdir(folder) if (os.path.isfile(os.path.join(folder, f)) and f != '.DS_Store')] # iterate over all files to extract the grids and put them in an ordered list of dicts From a5ee2e9a593d14340485cf6e5d30e3945301bb94 Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Mon, 3 Oct 2022 17:32:06 -0700 Subject: [PATCH 4/7] debugged LongTermAverage further --- .../gulf_of_mexico_LongTermAverageSource.yaml | 25 ++-- .../OceanCurrentSource/OceanCurrentSource.py | 25 ++-- scripts/tutorial/controller/hj_planner.py | 114 +++++++++++------- 3 files changed, 99 insertions(+), 65 deletions(-) diff --git a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml index 34897b8d..40dfcbe6 100644 --- a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml +++ b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml @@ -23,20 +23,25 @@ ocean_dict: folder: "data/hindcast_test" source: "HYCOM" type: "hindcast" + forecast: - field: 'OceanCurrents' - source: 'forecast_files' - source_settings: - folder: "data/forecast_test" - source: "HYCOM" - type: "forecast" - monthly_avg: field: 'OceanCurrents' source: 'longterm_average' source_settings: - folder: "data/monthly_average" - source: "HYCOM" - type: "hindcast" + forecast: + field: 'OceanCurrents' + source: 'forecast_files' + source_settings: + folder: "data/forecast_test" + source: "HYCOM" + type: "forecast" + average: + field: 'OceanCurrents' + source: 'longterm_average' + source_settings: + folder: "data/monthly_average" + source: "HYCOM" + type: "average" solar_dict: hindcast: null diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py index 172bccc9..8934d888 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py @@ -1,4 +1,5 @@ import datetime +from multiprocessing.sharedctypes import Value import os from typing import List, AnyStr, Optional, Union import logging @@ -45,7 +46,6 @@ def initialize_casadi_functions(self, grid: List[List[float]], array: xr) -> Non grid: list of the 3 grids [time, y_grid, x_grid] for the xr data array: xarray object containing the sub-setted data for the next cached round """ - self.u_curr_func = ca.interpolant('u_curr', 'linear', grid, array['water_u'].values.ravel(order='F')) self.v_curr_func = ca.interpolant('v_curr', 'linear', grid, array['water_v'].values.ravel(order='F')) @@ -312,8 +312,9 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean class LongTermAverageSource(OceanCurrentSource): # figure out inheritance def __init__(self, source_config_dict: dict): - self.forecast_data_source = ForecastFileSource(source_config_dict['forecast']) - self.monthly_avg_data_source = HindcastFileSource(source_config_dict['monthly_avg']) # defaults currents to normal + self.u_curr_func, self.v_curr_func = [None] * 2 + self.forecast_data_source = ForecastFileSource(source_config_dict['source_settings']['forecast']) + self.monthly_avg_data_source = HindcastFileSource(source_config_dict['source_settings']['average']) # defaults currents to normal self.source_config_dict = source_config_dict # self.t_0 = source_config_dict['t0'] # not sure what to do here @@ -323,14 +324,20 @@ def get_data_over_area(self, x_interval: List[float], y_interval: List[float], temporal_resolution: Optional[float] = None) -> xr: # Query as much forecast data as is possible print("T_INTERVAL: " + str(t_interval)) - forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, spatial_resolution, temporal_resolution) - # last_time = forecast_dataframe["time"].resample(how="last") - # print("301") - # print(last_time) - remaining_t_interval = [get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]), t_interval[1]] # may not work + try: + forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, spatial_resolution, temporal_resolution) + end_forecast_time = get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]) + except ValueError: + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, t_interval, spatial_resolution, temporal_resolution) + return monthly_average_dataframe + + + if end_forecast_time >= t_interval[1]: + return forecast_dataframe + + remaining_t_interval = [end_forecast_time, t_interval[1]] # may not work print("REMAINING T INTERVAL: " + str(remaining_t_interval)) monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, remaining_t_interval, spatial_resolution, temporal_resolution) - return xr.concat([forecast_dataframe, monthly_average_dataframe], dim="time") def check_for_most_recent_fmrc_dataframe(self, time: datetime.datetime) -> int: diff --git a/scripts/tutorial/controller/hj_planner.py b/scripts/tutorial/controller/hj_planner.py index 5007ff2e..9ff4a4e5 100644 --- a/scripts/tutorial/controller/hj_planner.py +++ b/scripts/tutorial/controller/hj_planner.py @@ -31,57 +31,79 @@ x_T=x_T, deg_around_x0_xT_box=1, temp_horizon_in_s=3600) +# x_0.date_time = 2021-11-24 12:00:00+00:00 ax = arena.ocean_field.hindcast_data_source.plot_data_at_time_over_area( time=x_0.date_time, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) problem.plot(ax=ax) plt.show() +# %% Try to do the same as above for the longterm average +print("Longterm Average") +t_interval, lat_bnds, lon_bnds = arena.ocean_field.forecast_data_source.convert_to_x_y_time_bounds( + x_0=x_0.to_spatio_temporal_point(), + x_T=x_T, deg_around_x0_xT_box=1, + temp_horizon_in_s=3600) + +# x_0.date_time = 2021-11-24 12:00:00+00:00 -> should just plot the forecast stuff +time_forecast = datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc) +ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( + time=time_forecast, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) +problem.plot(ax=ax) +plt.show() + +time_average = datetime.datetime(2021, 12, 1, 12, 0, tzinfo=datetime.timezone.utc) +ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( + time=time_average, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) +problem.plot(ax=ax) +plt.show() + +d = arena.ocean_field.forecast_data_source.get_data_over_area(x_interval=lon_bnds, y_interval=lat_bnds, t_interval=[time_forecast, time_average]) # %% Instantiate the HJ Planner -specific_settings = { - 'replan_on_new_fmrc': True, - 'replan_every_X_seconds': False, - 'direction': 'multi-time-reach-back', - 'n_time_vector': 200, - # Note that this is the number of time-intervals, the vector is +1 longer because of init_time - 'deg_around_xt_xT_box': 1., # area over which to run HJ_reachability - 'accuracy': 'high', - 'artificial_dissipation_scheme': 'local_local', - 'T_goal_in_seconds': 3600 * 24 * 5, - 'use_geographic_coordinate_system': True, - 'progress_bar': True, - 'initial_set_radii': [0.1, 0.1], # this is in deg lat, lon. Note: for Backwards-Reachability this should be bigger. - # Note: grid_res should always be bigger than initial_set_radii, otherwise reachability behaves weirdly. - 'grid_res': 0.02, # Note: this is in deg lat, lon (HYCOM Global is 0.083 and Mexico 0.04) - 'd_max': 0.0, - # 'EVM_threshold': 0.3 # in m/s error when floating in forecasted vs sensed currents - # 'fwd_back_buffer_in_seconds': 0.5, # this is the time added to the earliest_to_reach as buffer for forward-backward - 'platform_dict': arena.platform.platform_dict -} -planner = HJReach2DPlanner(problem=problem, specific_settings=specific_settings) +# specific_settings = { +# 'replan_on_new_fmrc': True, +# 'replan_every_X_seconds': False, +# 'direction': 'multi-time-reach-back', +# 'n_time_vector': 200, +# # Note that this is the number of time-intervals, the vector is +1 longer because of init_time +# 'deg_around_xt_xT_box': 1., # area over which to run HJ_reachability +# 'accuracy': 'high', +# 'artificial_dissipation_scheme': 'local_local', +# 'T_goal_in_seconds': 3600 * 24 * 5, +# 'use_geographic_coordinate_system': True, +# 'progress_bar': True, +# 'initial_set_radii': [0.1, 0.1], # this is in deg lat, lon. Note: for Backwards-Reachability this should be bigger. +# # Note: grid_res should always be bigger than initial_set_radii, otherwise reachability behaves weirdly. +# 'grid_res': 0.02, # Note: this is in deg lat, lon (HYCOM Global is 0.083 and Mexico 0.04) +# 'd_max': 0.0, +# # 'EVM_threshold': 0.3 # in m/s error when floating in forecasted vs sensed currents +# # 'fwd_back_buffer_in_seconds': 0.5, # this is the time added to the earliest_to_reach as buffer for forward-backward +# 'platform_dict': arena.platform.platform_dict +# } +# planner = HJReach2DPlanner(problem=problem, specific_settings=specific_settings) -# % Run reachability planner -observation = arena.reset(platform_state=x_0) -action = planner.get_action(observation=observation) -# %% Various plotting of the reachability computations -planner.plot_reachability_snapshot(rel_time_in_seconds=0, granularity_in_h=5, alpha_color=1, time_to_reach=True, fig_size_inches=(12, 12), plot_in_h=True) -# planner.plot_reachability_snapshot_over_currents(rel_time_in_seconds=0, granularity_in_h=5, time_to_reach=False) -# planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, filename="test_reach_animation.mp4") -# planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, with_opt_ctrl=True, -# filename="test_reach_animation_w_ctrl.mp4", forward_time=True) -#%% save planner state and reload it -# Save it to a folder -# planner.save_planner_state("saved_planner/") -# Load it from the folder -# loaded_planner = HJReach2DPlanner.from_saved_planner_state(folder="saved_planner/", problem=problem) -# loaded_planner.last_data_source = arena.ocean_field.hindcast_data_source +# # % Run reachability planner # observation = arena.reset(platform_state=x_0) -# loaded_planner._update_current_data(observation=observation) -# planner = loaded_planner -# %% Let the controller run closed-loop within the arena (the simulation loop) -for i in tqdm(range(int(3600 * 24 * 3 / 600))): # 3 days - action = planner.get_action(observation=observation) - observation = arena.step(action) +# action = planner.get_action(observation=observation) +# # %% Various plotting of the reachability computations +# planner.plot_reachability_snapshot(rel_time_in_seconds=0, granularity_in_h=5, alpha_color=1, time_to_reach=True, fig_size_inches=(12, 12), plot_in_h=True) +# # planner.plot_reachability_snapshot_over_currents(rel_time_in_seconds=0, granularity_in_h=5, time_to_reach=False) +# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, filename="test_reach_animation.mp4") +# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, with_opt_ctrl=True, +# # filename="test_reach_animation_w_ctrl.mp4", forward_time=True) +# #%% save planner state and reload it +# # Save it to a folder +# # planner.save_planner_state("saved_planner/") +# # Load it from the folder +# # loaded_planner = HJReach2DPlanner.from_saved_planner_state(folder="saved_planner/", problem=problem) +# # loaded_planner.last_data_source = arena.ocean_field.hindcast_data_source +# # observation = arena.reset(platform_state=x_0) +# # loaded_planner._update_current_data(observation=observation) +# # planner = loaded_planner +# # %% Let the controller run closed-loop within the arena (the simulation loop) +# for i in tqdm(range(int(3600 * 24 * 3 / 600))): # 3 days +# action = planner.get_action(observation=observation) +# observation = arena.step(action) -#%% Plot the arena trajectory on the map -arena.plot_all_on_map(problem=problem) -#%% Animate the trajectory -arena.animate_trajectory(problem=problem, temporal_resolution=7200) +# #%% Plot the arena trajectory on the map +# arena.plot_all_on_map(problem=problem) +# #%% Animate the trajectory +# arena.animate_trajectory(problem=problem, temporal_resolution=7200) From 1d75c7761b9f73e2ec939b372dd08ba7b0d1e288 Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Sat, 8 Oct 2022 21:21:03 -0700 Subject: [PATCH 5/7] remove print statements --- .../data_sources/OceanCurrentSource/OceanCurrentSource.py | 2 -- scripts/tutorial/controller/hj_planner.py | 3 ++- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py index 8934d888..857a87f1 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py @@ -323,7 +323,6 @@ def get_data_over_area(self, x_interval: List[float], y_interval: List[float], spatial_resolution: Optional[float] = None, temporal_resolution: Optional[float] = None) -> xr: # Query as much forecast data as is possible - print("T_INTERVAL: " + str(t_interval)) try: forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, spatial_resolution, temporal_resolution) end_forecast_time = get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]) @@ -336,7 +335,6 @@ def get_data_over_area(self, x_interval: List[float], y_interval: List[float], return forecast_dataframe remaining_t_interval = [end_forecast_time, t_interval[1]] # may not work - print("REMAINING T INTERVAL: " + str(remaining_t_interval)) monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, remaining_t_interval, spatial_resolution, temporal_resolution) return xr.concat([forecast_dataframe, monthly_average_dataframe], dim="time") diff --git a/scripts/tutorial/controller/hj_planner.py b/scripts/tutorial/controller/hj_planner.py index 9ff4a4e5..826389a8 100644 --- a/scripts/tutorial/controller/hj_planner.py +++ b/scripts/tutorial/controller/hj_planner.py @@ -43,7 +43,7 @@ x_T=x_T, deg_around_x0_xT_box=1, temp_horizon_in_s=3600) -# x_0.date_time = 2021-11-24 12:00:00+00:00 -> should just plot the forecast stuff +# # x_0.date_time = 2021-11-24 12:00:00+00:00 -> should just plot the forecast stuff time_forecast = datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc) ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( time=time_forecast, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) @@ -57,6 +57,7 @@ plt.show() d = arena.ocean_field.forecast_data_source.get_data_over_area(x_interval=lon_bnds, y_interval=lat_bnds, t_interval=[time_forecast, time_average]) +print(d) # %% Instantiate the HJ Planner # specific_settings = { # 'replan_on_new_fmrc': True, From e0326bc973a5714aec7883f9e39d3deb1cca0f90 Mon Sep 17 00:00:00 2001 From: nishap1225 Date: Sat, 8 Oct 2022 21:23:15 -0700 Subject: [PATCH 6/7] move testing script to scripts/nisha folder and reset hj_planner.py --- scripts/nisha/avg_data_source.py | 110 ++++++++++++++++++++ scripts/tutorial/controller/hj_planner.py | 117 +++++++++------------- 2 files changed, 157 insertions(+), 70 deletions(-) create mode 100644 scripts/nisha/avg_data_source.py diff --git a/scripts/nisha/avg_data_source.py b/scripts/nisha/avg_data_source.py new file mode 100644 index 00000000..826389a8 --- /dev/null +++ b/scripts/nisha/avg_data_source.py @@ -0,0 +1,110 @@ +import datetime +import numpy as np +from tqdm import tqdm + +from ocean_navigation_simulator.environment.ArenaFactory import ArenaFactory +from ocean_navigation_simulator.environment.Platform import PlatformState +from ocean_navigation_simulator.environment.PlatformState import SpatialPoint +from ocean_navigation_simulator.environment.NavigationProblem import NavigationProblem +from ocean_navigation_simulator.controllers.hj_planners.HJReach2DPlanner import HJReach2DPlanner +from ocean_navigation_simulator.utils import units +import matplotlib.pyplot as plt + +# Initialize the Arena (holds all data sources and the platform, everything except controller) +arena = ArenaFactory.create(scenario_name='gulf_of_mexico_LongTermAverageSource') +# we can also download the respective files directly to a temp folder, then t_interval needs to be set +# % Specify Navigation Problem +x_0 = PlatformState(lon=units.Distance(deg=-82.5), lat=units.Distance(deg=23.7), + date_time=datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc)) +x_T = SpatialPoint(lon=units.Distance(deg=-80.3), lat=units.Distance(deg=24.6)) + +problem = NavigationProblem( + start_state=x_0, + end_region=x_T, + target_radius=0.1, + timeout=datetime.timedelta(days=2), + platform_dict=arena.platform.platform_dict) + +# %% Plot the problem on the map +t_interval, lat_bnds, lon_bnds = arena.ocean_field.hindcast_data_source.convert_to_x_y_time_bounds( + x_0=x_0.to_spatio_temporal_point(), + x_T=x_T, deg_around_x0_xT_box=1, + temp_horizon_in_s=3600) + +# x_0.date_time = 2021-11-24 12:00:00+00:00 +ax = arena.ocean_field.hindcast_data_source.plot_data_at_time_over_area( + time=x_0.date_time, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) +problem.plot(ax=ax) +plt.show() +# %% Try to do the same as above for the longterm average +print("Longterm Average") +t_interval, lat_bnds, lon_bnds = arena.ocean_field.forecast_data_source.convert_to_x_y_time_bounds( + x_0=x_0.to_spatio_temporal_point(), + x_T=x_T, deg_around_x0_xT_box=1, + temp_horizon_in_s=3600) + +# # x_0.date_time = 2021-11-24 12:00:00+00:00 -> should just plot the forecast stuff +time_forecast = datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc) +ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( + time=time_forecast, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) +problem.plot(ax=ax) +plt.show() + +time_average = datetime.datetime(2021, 12, 1, 12, 0, tzinfo=datetime.timezone.utc) +ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( + time=time_average, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) +problem.plot(ax=ax) +plt.show() + +d = arena.ocean_field.forecast_data_source.get_data_over_area(x_interval=lon_bnds, y_interval=lat_bnds, t_interval=[time_forecast, time_average]) +print(d) +# %% Instantiate the HJ Planner +# specific_settings = { +# 'replan_on_new_fmrc': True, +# 'replan_every_X_seconds': False, +# 'direction': 'multi-time-reach-back', +# 'n_time_vector': 200, +# # Note that this is the number of time-intervals, the vector is +1 longer because of init_time +# 'deg_around_xt_xT_box': 1., # area over which to run HJ_reachability +# 'accuracy': 'high', +# 'artificial_dissipation_scheme': 'local_local', +# 'T_goal_in_seconds': 3600 * 24 * 5, +# 'use_geographic_coordinate_system': True, +# 'progress_bar': True, +# 'initial_set_radii': [0.1, 0.1], # this is in deg lat, lon. Note: for Backwards-Reachability this should be bigger. +# # Note: grid_res should always be bigger than initial_set_radii, otherwise reachability behaves weirdly. +# 'grid_res': 0.02, # Note: this is in deg lat, lon (HYCOM Global is 0.083 and Mexico 0.04) +# 'd_max': 0.0, +# # 'EVM_threshold': 0.3 # in m/s error when floating in forecasted vs sensed currents +# # 'fwd_back_buffer_in_seconds': 0.5, # this is the time added to the earliest_to_reach as buffer for forward-backward +# 'platform_dict': arena.platform.platform_dict +# } +# planner = HJReach2DPlanner(problem=problem, specific_settings=specific_settings) + +# # % Run reachability planner +# observation = arena.reset(platform_state=x_0) +# action = planner.get_action(observation=observation) +# # %% Various plotting of the reachability computations +# planner.plot_reachability_snapshot(rel_time_in_seconds=0, granularity_in_h=5, alpha_color=1, time_to_reach=True, fig_size_inches=(12, 12), plot_in_h=True) +# # planner.plot_reachability_snapshot_over_currents(rel_time_in_seconds=0, granularity_in_h=5, time_to_reach=False) +# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, filename="test_reach_animation.mp4") +# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, with_opt_ctrl=True, +# # filename="test_reach_animation_w_ctrl.mp4", forward_time=True) +# #%% save planner state and reload it +# # Save it to a folder +# # planner.save_planner_state("saved_planner/") +# # Load it from the folder +# # loaded_planner = HJReach2DPlanner.from_saved_planner_state(folder="saved_planner/", problem=problem) +# # loaded_planner.last_data_source = arena.ocean_field.hindcast_data_source +# # observation = arena.reset(platform_state=x_0) +# # loaded_planner._update_current_data(observation=observation) +# # planner = loaded_planner +# # %% Let the controller run closed-loop within the arena (the simulation loop) +# for i in tqdm(range(int(3600 * 24 * 3 / 600))): # 3 days +# action = planner.get_action(observation=observation) +# observation = arena.step(action) + +# #%% Plot the arena trajectory on the map +# arena.plot_all_on_map(problem=problem) +# #%% Animate the trajectory +# arena.animate_trajectory(problem=problem, temporal_resolution=7200) diff --git a/scripts/tutorial/controller/hj_planner.py b/scripts/tutorial/controller/hj_planner.py index 826389a8..e0cc8b28 100644 --- a/scripts/tutorial/controller/hj_planner.py +++ b/scripts/tutorial/controller/hj_planner.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt # Initialize the Arena (holds all data sources and the platform, everything except controller) -arena = ArenaFactory.create(scenario_name='gulf_of_mexico_LongTermAverageSource') +arena = ArenaFactory.create(scenario_name='gulf_of_mexico_HYCOM_hindcast_local') # we can also download the respective files directly to a temp folder, then t_interval needs to be set # % Specify Navigation Problem x_0 = PlatformState(lon=units.Distance(deg=-82.5), lat=units.Distance(deg=23.7), @@ -31,80 +31,57 @@ x_T=x_T, deg_around_x0_xT_box=1, temp_horizon_in_s=3600) -# x_0.date_time = 2021-11-24 12:00:00+00:00 ax = arena.ocean_field.hindcast_data_source.plot_data_at_time_over_area( time=x_0.date_time, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) problem.plot(ax=ax) plt.show() -# %% Try to do the same as above for the longterm average -print("Longterm Average") -t_interval, lat_bnds, lon_bnds = arena.ocean_field.forecast_data_source.convert_to_x_y_time_bounds( - x_0=x_0.to_spatio_temporal_point(), - x_T=x_T, deg_around_x0_xT_box=1, - temp_horizon_in_s=3600) - -# # x_0.date_time = 2021-11-24 12:00:00+00:00 -> should just plot the forecast stuff -time_forecast = datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc) -ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( - time=time_forecast, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) -problem.plot(ax=ax) -plt.show() - -time_average = datetime.datetime(2021, 12, 1, 12, 0, tzinfo=datetime.timezone.utc) -ax = arena.ocean_field.forecast_data_source.plot_data_at_time_over_area( - time=time_average, x_interval=lon_bnds, y_interval=lat_bnds, return_ax=True) -problem.plot(ax=ax) -plt.show() - -d = arena.ocean_field.forecast_data_source.get_data_over_area(x_interval=lon_bnds, y_interval=lat_bnds, t_interval=[time_forecast, time_average]) -print(d) # %% Instantiate the HJ Planner -# specific_settings = { -# 'replan_on_new_fmrc': True, -# 'replan_every_X_seconds': False, -# 'direction': 'multi-time-reach-back', -# 'n_time_vector': 200, -# # Note that this is the number of time-intervals, the vector is +1 longer because of init_time -# 'deg_around_xt_xT_box': 1., # area over which to run HJ_reachability -# 'accuracy': 'high', -# 'artificial_dissipation_scheme': 'local_local', -# 'T_goal_in_seconds': 3600 * 24 * 5, -# 'use_geographic_coordinate_system': True, -# 'progress_bar': True, -# 'initial_set_radii': [0.1, 0.1], # this is in deg lat, lon. Note: for Backwards-Reachability this should be bigger. -# # Note: grid_res should always be bigger than initial_set_radii, otherwise reachability behaves weirdly. -# 'grid_res': 0.02, # Note: this is in deg lat, lon (HYCOM Global is 0.083 and Mexico 0.04) -# 'd_max': 0.0, -# # 'EVM_threshold': 0.3 # in m/s error when floating in forecasted vs sensed currents -# # 'fwd_back_buffer_in_seconds': 0.5, # this is the time added to the earliest_to_reach as buffer for forward-backward -# 'platform_dict': arena.platform.platform_dict -# } -# planner = HJReach2DPlanner(problem=problem, specific_settings=specific_settings) +specific_settings = { + 'replan_on_new_fmrc': True, + 'replan_every_X_seconds': False, + 'direction': 'multi-time-reach-back', + 'n_time_vector': 200, + # Note that this is the number of time-intervals, the vector is +1 longer because of init_time + 'deg_around_xt_xT_box': 1., # area over which to run HJ_reachability + 'accuracy': 'high', + 'artificial_dissipation_scheme': 'local_local', + 'T_goal_in_seconds': 3600 * 24 * 5, + 'use_geographic_coordinate_system': True, + 'progress_bar': True, + 'initial_set_radii': [0.1, 0.1], # this is in deg lat, lon. Note: for Backwards-Reachability this should be bigger. + # Note: grid_res should always be bigger than initial_set_radii, otherwise reachability behaves weirdly. + 'grid_res': 0.02, # Note: this is in deg lat, lon (HYCOM Global is 0.083 and Mexico 0.04) + 'd_max': 0.0, + # 'EVM_threshold': 0.3 # in m/s error when floating in forecasted vs sensed currents + # 'fwd_back_buffer_in_seconds': 0.5, # this is the time added to the earliest_to_reach as buffer for forward-backward + 'platform_dict': arena.platform.platform_dict +} +planner = HJReach2DPlanner(problem=problem, specific_settings=specific_settings) -# # % Run reachability planner +# % Run reachability planner +observation = arena.reset(platform_state=x_0) +action = planner.get_action(observation=observation) +# %% Various plotting of the reachability computations +planner.plot_reachability_snapshot(rel_time_in_seconds=0, granularity_in_h=5, alpha_color=1, time_to_reach=True, fig_size_inches=(12, 12), plot_in_h=True) +# planner.plot_reachability_snapshot_over_currents(rel_time_in_seconds=0, granularity_in_h=5, time_to_reach=False) +# planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, filename="test_reach_animation.mp4") +# planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, with_opt_ctrl=True, +# filename="test_reach_animation_w_ctrl.mp4", forward_time=True) +#%% save planner state and reload it +# Save it to a folder +# planner.save_planner_state("saved_planner/") +# Load it from the folder +# loaded_planner = HJReach2DPlanner.from_saved_planner_state(folder="saved_planner/", problem=problem) +# loaded_planner.last_data_source = arena.ocean_field.hindcast_data_source # observation = arena.reset(platform_state=x_0) -# action = planner.get_action(observation=observation) -# # %% Various plotting of the reachability computations -# planner.plot_reachability_snapshot(rel_time_in_seconds=0, granularity_in_h=5, alpha_color=1, time_to_reach=True, fig_size_inches=(12, 12), plot_in_h=True) -# # planner.plot_reachability_snapshot_over_currents(rel_time_in_seconds=0, granularity_in_h=5, time_to_reach=False) -# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, filename="test_reach_animation.mp4") -# # planner.plot_reachability_animation(time_to_reach=True, granularity_in_h=5, with_opt_ctrl=True, -# # filename="test_reach_animation_w_ctrl.mp4", forward_time=True) -# #%% save planner state and reload it -# # Save it to a folder -# # planner.save_planner_state("saved_planner/") -# # Load it from the folder -# # loaded_planner = HJReach2DPlanner.from_saved_planner_state(folder="saved_planner/", problem=problem) -# # loaded_planner.last_data_source = arena.ocean_field.hindcast_data_source -# # observation = arena.reset(platform_state=x_0) -# # loaded_planner._update_current_data(observation=observation) -# # planner = loaded_planner -# # %% Let the controller run closed-loop within the arena (the simulation loop) -# for i in tqdm(range(int(3600 * 24 * 3 / 600))): # 3 days -# action = planner.get_action(observation=observation) -# observation = arena.step(action) +# loaded_planner._update_current_data(observation=observation) +# planner = loaded_planner +# %% Let the controller run closed-loop within the arena (the simulation loop) +for i in tqdm(range(int(3600 * 24 * 3 / 600))): # 3 days + action = planner.get_action(observation=observation) + observation = arena.step(action) -# #%% Plot the arena trajectory on the map -# arena.plot_all_on_map(problem=problem) -# #%% Animate the trajectory -# arena.animate_trajectory(problem=problem, temporal_resolution=7200) +#%% Plot the arena trajectory on the map +arena.plot_all_on_map(problem=problem) +#%% Animate the trajectory +arena.animate_trajectory(problem=problem, temporal_resolution=7200) \ No newline at end of file From 864b2f6c455a38dc04d25767fcadc1e033002550 Mon Sep 17 00:00:00 2001 From: MariusWiggert Date: Fri, 14 Oct 2022 14:38:13 -0700 Subject: [PATCH 7/7] small iterations on Nishas long-term average code --- .../gulf_of_mexico_LongTermAverageSource.yaml | 6 +- .../data_sources/DataSource.py | 8 +- .../OceanCurrentSource/OceanCurrentSource.py | 21 ++- .../data_sources/ocean_current_sources.py | 161 ++++++++++++++++++ 4 files changed, 183 insertions(+), 13 deletions(-) diff --git a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml index 40dfcbe6..14029d71 100644 --- a/config/arena/gulf_of_mexico_LongTermAverageSource.yaml +++ b/config/arena/gulf_of_mexico_LongTermAverageSource.yaml @@ -33,15 +33,11 @@ ocean_dict: source: 'forecast_files' source_settings: folder: "data/forecast_test" - source: "HYCOM" - type: "forecast" average: field: 'OceanCurrents' - source: 'longterm_average' + source: 'hindcast_files' source_settings: folder: "data/monthly_average" - source: "HYCOM" - type: "average" solar_dict: hindcast: null diff --git a/ocean_navigation_simulator/data_sources/DataSource.py b/ocean_navigation_simulator/data_sources/DataSource.py index 0c783c11..3c2fdbc6 100644 --- a/ocean_navigation_simulator/data_sources/DataSource.py +++ b/ocean_navigation_simulator/data_sources/DataSource.py @@ -170,14 +170,12 @@ def array_subsetting_sanity_check(array: xr, x_interval: List[float], y_interval # Step 2: Data partially not in the array check if array.coords['lat'].data[0] > y_interval[0] or array.coords['lat'].data[-1] < y_interval[1]: logger.warning( - f"Part of the y requested area is outside of file(file: [{array.coords['lat'].data[0]}, {array.coords['lat'].data[-1]}], requested: [{y_interval[0]}, {y_interval[1]}]).", - RuntimeWarning) + f"Part of the y requested area is outside of file(file: [{array.coords['lat'].data[0]}, {array.coords['lat'].data[-1]}], requested: [{y_interval[0]}, {y_interval[1]}]).") if array.coords['lon'].data[0] > x_interval[0] or array.coords['lon'].data[-1] < x_interval[1]: logger.warning( - f"Part of the x requested area is outside of file (file: [{array.coords['lon'].data[0]}, {array.coords['lon'].data[-1]}], requested: [{x_interval[0]}, {x_interval[1]}]).", - RuntimeWarning) + f"Part of the x requested area is outside of file (file: [{array.coords['lon'].data[0]}, {array.coords['lon'].data[-1]}], requested: [{x_interval[0]}, {x_interval[1]}]).") if units.get_datetime_from_np64(array.coords['time'].data[-1]) < t_interval[1]: - logger.warning("The final time is not part of the subset.".format(t_interval[1]), RuntimeWarning) + logger.warning("The final time {} is not part of the subset.".format(t_interval[1])) def plot_data_at_time_over_area(self, time: Union[datetime.datetime, float], x_interval: List[float], y_interval: List[float], diff --git a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py index 857a87f1..7aa31f2c 100644 --- a/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py +++ b/ocean_navigation_simulator/data_sources/OceanCurrentSource/OceanCurrentSource.py @@ -186,7 +186,7 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean class ForecastFileSource(OceanCurrentSourceXarray): # TODO: Make it work with multiple files for one forecast (a bit of extra logic, but possible) - """Data Source Object that accesses and manages multiple daily HYCOM files as source.""" + """Data Source Object that accesses and manages multiple daily files as source.""" def __init__(self, source_config_dict: dict): super().__init__(source_config_dict) @@ -210,6 +210,18 @@ def get_data_over_area(self, x_interval: List[float], y_interval: List[float], spatial_resolution: Optional[float] = None, temporal_resolution: Optional[float] = None, most_recent_fmrc_at_time: Optional[datetime.datetime] = None) -> xr: + """Function to get the the raw current data over an x, y, and t interval. + Args: + x_interval: List of the lower and upper x area in the respective coordinate units [x_lower, x_upper] + y_interval: List of the lower and upper y area in the respective coordinate units [y_lower, y_upper] + t_interval: List of the lower and upper datetime requested [t_0, t_T] in datetime + spatial_resolution: spatial resolution in the same units as x and y interval + temporal_resolution: temporal resolution in seconds + most_recent_fmrc_at_time: if specified this is the idx of a specific forecast file to get data from it + otherwise the most recent fmrc available at t_interval[0] is used. + Returns: + data_array in xarray format that contains the grid and the values (land is NaN) + """ # format to datetime object if not isinstance(t_interval[0], datetime.datetime): t_interval = [datetime.datetime.fromtimestamp(time, tz=datetime.timezone.utc) for time in t_interval] @@ -267,7 +279,7 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean class HindcastFileSource(OceanCurrentSourceXarray): - """Data Source Object that accesses and manages one or many HYCOM files as source.""" + """Data Source Object that accesses and manages one or many daily files as source.""" def __init__(self, source_config_dict: dict): super().__init__(source_config_dict) @@ -290,6 +302,7 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean class HindcastOpendapSource(OceanCurrentSourceXarray): + """Data Source Object that accesses the data via the opendap framework directly from the server.""" def __init__(self, source_config_dict: dict): super().__init__(source_config_dict) # Step 1: establish the opendap connection with the settings in config_dict @@ -310,7 +323,9 @@ def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> Ocean return OceanCurrentVector(u=self.u_curr_func(spatio_temporal_point.to_spatio_temporal_casadi_input()), v=self.v_curr_func(spatio_temporal_point.to_spatio_temporal_casadi_input())) -class LongTermAverageSource(OceanCurrentSource): # figure out inheritance + +class LongTermAverageSource(OceanCurrentSource): + """""" def __init__(self, source_config_dict: dict): self.u_curr_func, self.v_curr_func = [None] * 2 self.forecast_data_source = ForecastFileSource(source_config_dict['source_settings']['forecast']) diff --git a/scripts/tutorial/data_sources/ocean_current_sources.py b/scripts/tutorial/data_sources/ocean_current_sources.py index b09033cb..57122161 100644 --- a/scripts/tutorial/data_sources/ocean_current_sources.py +++ b/scripts/tutorial/data_sources/ocean_current_sources.py @@ -7,6 +7,167 @@ # For fast interpolation of currents we cache part of the spatio-temporal data around x_t in a casadi function casadi_cache_dict = {'deg_around_x_t': 1, 'time_around_x_t': 3600 * 24 * 1} +import yaml +with open(f'config/arena/gulf_of_mexico_LongTermAverageSource.yaml') as f: + full_config = yaml.load(f, Loader=yaml.FullLoader) + +source_config_dict = full_config['ocean_dict']['forecast'] +forecast_dict = source_config_dict['source_settings']['forecast'] +average_dict = source_config_dict['source_settings']['average'] +for source_dict in [forecast_dict, average_dict]: + source_dict['casadi_cache_settings'] = casadi_cache_dict + source_dict['use_geographic_coordinate_system'] = True +#%% +from ocean_navigation_simulator.data_sources.OceanCurrentSource.OceanCurrentSource import OceanCurrentSource, ForecastFileSource, HindcastFileSource, get_datetime_from_np64 + +# Step 1: Initialize both data_sources +forecast_data_source = ForecastFileSource(forecast_dict) +monthly_avg_data_source = HindcastFileSource(average_dict) +#%% plot to see if it worked +# forecast_data_source.plot_data_at_time_over_area(datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc), +# [-84, -82], [20,25]) +# monthly_avg_data_source.plot_data_at_time_over_area(datetime.datetime(2021, 11, 30, 12, 0, tzinfo=datetime.timezone.utc), +# [-84, -82], [20,25]) +#%% +x_interval= [-84, -82] +y_interval= [20,25] +t_interval= [datetime.datetime(2021, 11, 24, 12, 0, tzinfo=datetime.timezone.utc), + datetime.datetime(2021, 11, 30, 12, 0, tzinfo=datetime.timezone.utc)] +spatial_resolution = 0.1 +temporal_resolution= 3600*6 +#%% +# todo: use keyword inputs for clarity! +forecast_dataframe = forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval,spatial_resolution, temporal_resolution) +# 11-24 to 11-28 +# Now get end_forecast time +end_forecast_time = get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]) +if end_forecast_time >= t_interval[1]: + print("forecast is enough") + +# easiest fix: run it with temp and spat resolution of the forecast... +remaining_t_interval = [end_forecast_time, t_interval[1]] +monthly_average_dataframe = monthly_avg_data_source.get_data_over_area(x_interval, y_interval, + remaining_t_interval, + spatial_resolution, + temporal_resolution) +#%% Now cut out the relevant times: +subset = monthly_average_dataframe.sel( + time=slice(remaining_t_interval[0], remaining_t_interval[1])) +#%% now concat +import xarray as xr +full_frame = xr.concat([forecast_dataframe, subset], dim="time") +# Works but right now it's multiple resolutions! (hourly vs monthly data...) +#%% render it +import datetime +import matplotlib.animation as animation +import numpy as np +import xarray as xr +from functools import partial + +# Step 1: get the data_subset for animation +xarray = full_frame + +# Calculate min and max over the full tempo-spatial array +# get rounded up vmax across the whole dataset (with ` decimals) +xarray = xarray.assign(magnitude=lambda x: (x.water_u ** 2 + x.water_v ** 2) ** 0.5) +vmax = round(xarray['magnitude'].max().item() + 0.049, 1) +vmin = 0 + + +# create global figure object where the animation happens +import matplotlib.pyplot as plt +fig = plt.figure(figsize=(12, 12)) + +render_frame = partial(forecast_data_source.plot_xarray_for_animation, xarray=xarray, vmin=vmin, vmax=vmax, + reset_plot=True) +# set time direction of the animation +frames_vector = np.where(True, np.arange(xarray['time'].size), np.flip(np.arange(xarray['time'].size))) +# create animation function object (it's not yet executed) +ani = animation.FuncAnimation(fig, func=render_frame, frames=frames_vector, repeat=False) + +# render the animation with the keyword arguments +forecast_data_source.render_animation(animation_object=ani, output="average.mp4", fps=10) +#%% +full_frame.sel(time=slice(remaining_t_interval[0] - datetime.timedelta(hours=10), + remaining_t_interval[0] + datetime.timedelta(hours=10))) +#%% + # Query as much forecast data as is possible + try: + forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, + spatial_resolution, temporal_resolution) + end_forecast_time = get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]) + except ValueError: + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, + t_interval, spatial_resolution, + temporal_resolution) + return monthly_average_dataframe + + if end_forecast_time >= t_interval[1]: + return forecast_dataframe + + remaining_t_interval = [end_forecast_time, t_interval[1]] # may not work + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, + remaining_t_interval, + spatial_resolution, + temporal_resolution) + return xr.concat([forecast_dataframe, monthly_average_dataframe], dim="time") + + +#%% +class LongTermAverageSource(OceanCurrentSource): + """""" + + def __init__(self, source_config_dict: dict): + self.u_curr_func, self.v_curr_func = [None] * 2 + self.forecast_data_source = ForecastFileSource(source_config_dict['source_settings']['forecast']) + self.monthly_avg_data_source = HindcastFileSource( + source_config_dict['source_settings']['average']) # defaults currents to normal + self.source_config_dict = source_config_dict + # self.t_0 = source_config_dict['t0'] # not sure what to do here + + def get_data_over_area(self, x_interval: List[float], y_interval: List[float], + t_interval: List[Union[datetime.datetime, int]], + spatial_resolution: Optional[float] = None, + temporal_resolution: Optional[float] = None) -> xr: + # Query as much forecast data as is possible + try: + forecast_dataframe = self.forecast_data_source.get_data_over_area(x_interval, y_interval, t_interval, + spatial_resolution, temporal_resolution) + end_forecast_time = get_datetime_from_np64(forecast_dataframe["time"].to_numpy()[-1]) + except ValueError: + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, + t_interval, spatial_resolution, + temporal_resolution) + return monthly_average_dataframe + + if end_forecast_time >= t_interval[1]: + return forecast_dataframe + + remaining_t_interval = [end_forecast_time, t_interval[1]] # may not work + monthly_average_dataframe = self.monthly_avg_data_source.get_data_over_area(x_interval, y_interval, + remaining_t_interval, + spatial_resolution, + temporal_resolution) + return xr.concat([forecast_dataframe, monthly_average_dataframe], dim="time") + + def check_for_most_recent_fmrc_dataframe(self, time: datetime.datetime) -> int: + """Helper function to check update the self.OceanCurrent if a new forecast is available at + the specified input time. + Args: + time: datetime object + """ + return self.forecast_data_source.check_for_most_recent_fmrc_dataframe(time) + + # Not sure if I can just all this + def get_data_at_point(self, spatio_temporal_point: SpatioTemporalPoint) -> OceanCurrentVector: + """We overwrite it because we don't want that Forecast needs caching...""" + return self.forecast_data_source.get_data_at_point(spatio_temporal_point == spatio_temporal_point) + + + + + + #%% Create the source dict for the ocean current #%% Option 1: Accessing data in the Copernicus (or HYCOM) server directly via opendap -> data loaded when needed