Skip to content

Commit e839285

Browse files
authored
Merge pull request #1773 from knutfrode/dev
[run-ex] It is now possible to provide element-specific environment c…
2 parents 33b871f + 2cc8ebb commit e839285

5 files changed

Lines changed: 105 additions & 17 deletions

File tree

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
#!/usr/bin/env python
2+
"""
3+
Element dependent environment
4+
=============================
5+
"""
6+
7+
from datetime import datetime
8+
import matplotlib.pyplot as plt
9+
import numpy as np
10+
from opendrift.models.oceandrift import OceanDrift
11+
import trajan as ta
12+
13+
#%%
14+
# OpenDrift elements have properties such as lon, lat, z, size etc.
15+
# These elements are moved and changed based on environment properties such as current, wind, temperature etc.
16+
# In principle the element properties do not affect the environment, however,
17+
# for sensitivity studies it can be of interest to let different elements of
18+
# the same simulation be exposed to different environment.
19+
# This is made possible by specifying the (constant) environment values
20+
# for all elements of a seed call as illustrated below
21+
22+
23+
#%%
24+
# First an example with two different seedings with two different environments
25+
o = OceanDrift(loglevel=20)
26+
# First seeding 200 elements that will be exposed to eastward current and a horizontal diffusivity of 10 m2/s
27+
number = 200
28+
o.seed_elements(lon=-60, lat=40, time=datetime(2022,1,1), number=number, radius=10,
29+
environment={'horizontal_diffusivity': 10,
30+
'x_sea_water_velocity': 1,
31+
'y_sea_water_velocity': 0})
32+
# Then seeding 100 elements that will be exposed to northward current and less diffusivity (1 m2/s)
33+
number = 100
34+
o.seed_elements(lon=-60, lat=40, time=datetime(2022,1,1), number=number, radius=10,
35+
environment={'horizontal_diffusivity': 1,
36+
'x_sea_water_velocity': 0,
37+
'y_sea_water_velocity': .5})
38+
o.run(steps=10)
39+
o.plot()
40+
41+
#%%
42+
# Second example with a single seeding where each element will be exposed to different diffusivity values
43+
o = OceanDrift(loglevel=20)
44+
# Seeding 1000 elements that will be exposed to north-eastward current with diffusivities ranging from 0 to 50 m2/s
45+
number = 1000
46+
diffusivity_values = [0, 1, 5, 10, 50]
47+
# Repeating values so that 200 elements get each diffusivity
48+
diffusivities = np.repeat(diffusivity_values, number/len(diffusivity_values))
49+
50+
o.seed_elements(lon=-60, lat=40, time=datetime(2022,1,1), number=number, radius=10,
51+
environment={'horizontal_diffusivity': diffusivities,
52+
'x_sea_water_velocity': .2,
53+
'y_sea_water_velocity': .2})
54+
ds = o.run(steps=10)
55+
ds.traj.plot(land=None, margin=0)
56+
# Plotting the convex hull around end positions separately for each diffusivity value
57+
ds = ds.isel(time=-1)
58+
colors = plt.cm.jet(np.linspace(0, 1, 5))
59+
for d, color in zip(diffusivity_values, colors):
60+
ds.where(ds.horizontal_diffusivity==d).traj.plot.convex_hull(label=f'Diffusivity {d} m2/s', color=color)
61+
plt.legend()
62+
plt.show()
63+
64+
#%%
65+
# The above could also be achieved by performing separate simulations for each value of diffusivity,
66+
# but with more computational overhead/time and more complexity

opendrift/models/basemodel/__init__.py

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -904,6 +904,8 @@ def schedule_elements(self, elements, time):
904904
self.start_time = min_time
905905
logger.debug('Setting simulation start time to %s' % str(min_time))
906906

907+
return elements.ID
908+
907909
def release_elements(self):
908910
"""Activate elements which are scheduled within following timestep."""
909911

@@ -1215,10 +1217,23 @@ def seed_elements(self,
12151217
if prop in self.ElementType.variables:
12161218
kwargs[prop] = seed_config[f'seed:{prop}']['value']
12171219

1220+
environment = kwargs.pop('environment', None)
1221+
12181222
# Creating and scheduling elements
12191223
elements = self.ElementType(lon=lon, lat=lat, **kwargs)
12201224
time_array = np.array(time)
1221-
self.schedule_elements(elements, time)
1225+
new_element_IDs = self.schedule_elements(elements, time)
1226+
1227+
if environment is not None:
1228+
logger.debug(f'enviroment is provided to seed method, adding constant reader for {environment}')
1229+
from opendrift.readers.reader_constant import Reader as ConstantReader
1230+
environment['element_ID'] = new_element_IDs
1231+
cr = ConstantReader(environment)
1232+
# Must change mode tomporarily to be allowed to add a new reader
1233+
tmp_mode = self.mode
1234+
self.mode = opendrift.models.basemodel.Mode.Config
1235+
self.add_reader(cr)
1236+
self.mode = tmp_mode
12221237

12231238
@require_mode(mode=Mode.Ready)
12241239
def seed_cone(self, lon, lat, time, radius=0, number=None, **kwargs):

opendrift/models/basemodel/environment.py

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -621,7 +621,7 @@ def get_environment(self,
621621
# Fetch given variables at given positions from current reader
622622
if reader._element_ID is not None:
623623
logger.debug(f'Setting _element_ID for reader {reader_name}')
624-
reader._element_ID = element_ID
624+
reader._element_ID = element_ID[missing_indices]
625625
try:
626626
logger.debug('Data needed for %i elements' %
627627
len(missing_indices))
@@ -638,8 +638,7 @@ def get_environment(self,
638638
variable_group, profiles = profiles_from_reader,
639639
profiles_depth = profiles_depth, time = time,
640640
lon=lon[missing_indices], lat=lat[missing_indices],
641-
z=z[missing_indices], rotate_to_proj=self.proj_latlon,
642-
element_ID=element_ID)
641+
z=z[missing_indices], rotate_to_proj=self.proj_latlon)
643642

644643
except NotCoveredError as e:
645644
logger.info(e)

opendrift/readers/basereader/variables.py

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -687,7 +687,7 @@ def __check_env_arrays__(self, env):
687687

688688
@abstractmethod
689689
def _get_variables_interpolated_(self, variables, profiles, profiles_depth,
690-
time, reader_x, reader_y, z, element_ID=None):
690+
time, reader_x, reader_y, z):
691691
"""
692692
This method _must_ be implemented by every reader. Usually by
693693
subclassing one of the reader types (e.g.
@@ -710,8 +710,7 @@ def get_variables_interpolated_xy(self,
710710
x=None,
711711
y=None,
712712
z=None,
713-
rotate_to_proj=None,
714-
element_ID=None):
713+
rotate_to_proj=None):
715714
"""
716715
Get variables in native projection of reader.
717716
@@ -862,8 +861,7 @@ def get_variables_interpolated(self,
862861
lon=None,
863862
lat=None,
864863
z=None,
865-
rotate_to_proj=None,
866-
element_ID=None):
864+
rotate_to_proj=None):
867865
"""
868866
`get_variables_interpolated` is the main interface to
869867
:class:`opendrift.basemodel.OpenDriftSimulation`, and is responsible

opendrift/readers/reader_constant.py

Lines changed: 18 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -25,16 +25,18 @@ def __init__(self, parameter_value_map):
2525
"""init with a map {'variable_name': value, ...}
2626
2727
value can also be an array, and in this case the map/dictionary
28-
should also include `element_ID` which corresponds to the elements that
28+
must also include `element_ID` which corresponds to the elements that
2929
shall receive the actual value:
3030
self.environment.<variable_name> --> value[element_ID = self.elements.ID] (pseudo code)
31+
This is however more simply achived by specifying environment when seeding, see:
32+
https://opendrift.github.io/gallery/example_element_dependent_environment.html
3133
3234
"""
3335

3436
for key, var in parameter_value_map.items():
3537
parameter_value_map[key] = np.atleast_1d(var)
3638
self._parameter_value_map = parameter_value_map
37-
self.variables = list(parameter_value_map.keys())
39+
self.variables = list([v for v in parameter_value_map.keys() if v !='element_ID'])
3840
self.proj4 = '+proj=latlong'
3941
self.xmin = -180
4042
self.xmax = 180
@@ -55,14 +57,22 @@ def get_variables(self, requestedVariables, time=None,
5557
x=None, y=None, z=None):
5658

5759
variables = {'time': time, 'x': x, 'y': y, 'z': z}
58-
#variables.update(self._parameter_value_map)
60+
5961
for var in requestedVariables:
62+
variables[var] = np.nan*np.ones(x.shape) # Initialize with NaN
6063
value = self._parameter_value_map[var]
61-
if self._element_ID is None or len(self._parameter_value_map[var]==1): # Same scalar value for all elements
62-
variables[var] = self._parameter_value_map[var]*np.ones(x.shape)
63-
else: # Individual mapping
64-
indices = np.where(np.isin(self._parameter_value_map['element_ID'], self._element_ID))[0]
65-
variables[var] = self._parameter_value_map[var][indices]
64+
65+
if self._element_ID is None: # Same constant value for all elements
66+
variables[var] = value*np.ones(x.shape)
67+
continue
68+
69+
# Individual mapping
70+
indices = np.where(np.isin(self._parameter_value_map['element_ID'], self._element_ID))[0]
71+
ind_opp = np.where(np.isin(self._element_ID, self._parameter_value_map['element_ID']))[0]
72+
if len(value)==1: # Same value for all relevant elements
73+
variables[var][ind_opp] = value
74+
else:
75+
variables[var][ind_opp] = value[indices]
6676

6777
return variables
6878

0 commit comments

Comments
 (0)