Skip to content

Commit

Permalink
Merge pull request #1338 from Ugomartinez/dev
Browse files Browse the repository at this point in the history 
New features for plot function and improove of  wind measurement example
  • Loading branch information
@knutfrode
knutfrode authored Jun 27, 2024
2 parents 59bbd8f + 07b943b commit 38f858e
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 29 deletions.
63 changes: 41 additions & 22 deletions examples/example_wind_measurements.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,17 @@
import matplotlib.pyplot as plt
import matplotlib as mpl

# Preparation of figsize, dpi, fontsize and font style
# Could also be done using mpl style sheet
# https://matplotlib.org/stable/users/explain/customizing.html
mpl.rcParams['figure.figsize'] = (5, 5)
mpl.rcParams['savefig.dpi'] = 300
mpl.rcParams['font.family'] = "sans-serif"
mpl.rcParams['font.sans-serif'] = "DejaVu Sans"
mpl.rcParams['font.size'] = 22
mpl.rcParams['mathtext.fontset'] = "dejavusans"



#%%
# This example shows how to use wind or currents measurments
Expand All @@ -26,17 +35,19 @@


#%%
# We import two times the same model, because one
# Fisrt, we import two times the same model, because one
# will be used for the usual simulation, and one will
# be used for the measurement combined simulation.

reader_norkyst = reader_netCDF_CF_generic.Reader('https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')
reader_norkyst2 = reader_netCDF_CF_generic.Reader('https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')

#Verification
print(reader_norkyst)


#%%
# We import wind measurement data,
# We then import wind measurement data,
# that will be stored in a timeserie reader.
# Don't forget to give the lon and lat of the
# measurement station.
Expand All @@ -46,19 +57,17 @@
wind_speed_angle = np.array([22, 25, 24, 19, 17, 30, 16])


####################################
#File stored wind measurement data
#Get this type of file with printobs
#https://github.com/bohlinger/printobs

#data_wind_bridge = np.genfromtxt("Simulation_data/data_wind_21_mai_2024.txt", dtype = '<U10', skip_header=1)
#temporal_time = np.array([data_wind_bridge[i, 0] + "T" + data_wind_bridge[i, 1] for i in range(len(data_wind_bridge))])
#time_list = [datetime.fromisoformat(temporal_time[i]) for i in range(len(temporal_time))]
#wind_speed_array = data_wind_bridge[:, 3]
#wind_speed_angle = data_wind_bridge[:, 5]
#wind_speed_array = wind_speed_array.astype(float)
#wind_speed_angle = wind_speed_angle.astype(float)
##################################
# Alternative if data are in a csv file

# data_wind_bridge = np.genfromtxt("Simulation_data/data_wind_21_mai_2024.txt", dtype = '<U10', skip_header=1)
# temporal_time = np.array([data_wind_bridge[i, 0] + "T" + data_wind_bridge[i, 1] for i in range(len(data_wind_bridge))])
# time_list = [datetime.fromisoformat(temporal_time[i]) for i in range(len(temporal_time))]
# wind_speed_array = data_wind_bridge[:, 3]
# wind_speed_angle = data_wind_bridge[:, 5]
# wind_speed_array = wind_speed_array.astype(float)
# wind_speed_angle = wind_speed_angle.astype(float)


#Calculation of wind_x and wind_y
wind_speed_angle = np.deg2rad(wind_speed_angle)
Expand All @@ -69,6 +78,11 @@
lon, lat = 5.17, 60.37
parameter_value_map = {'time': time_list, 'wind_speed': wind_speed_array, 'x_wind': wind_x, 'y_wind': wind_y, 'lon':lon, 'lat':lat}
wind_obs_reader = reader_timeseries.Reader(parameter_value_map)

#Preparing future plots
text = [{'s': '* Station', 'x': lon, 'y': lat, 'fontsize': 20, 'color': 'k', 'backgroundcolor': 'white', 'bbox': dict(facecolor='none', edgecolor='none', alpha=0.4), 'zorder': 1000}]

#Verification
print(wind_obs_reader)

#%%
Expand All @@ -89,7 +103,7 @@
time_step = 150
time_step_output = 150

#Seeding configuration
# Seeding configuration
lon = 5.187
lat = 60.385
radius = 0
Expand All @@ -99,17 +113,17 @@


#%%
#Usual simulation
# Usual simulation

o = OceanDrift(loglevel=20)
o.add_reader([reader_norkyst])

#o.set_config('drift:current_uncertainty', .1)
#o.set_config('drift:wind_uncertainty', 1)
# o.set_config('drift:current_uncertainty', .1)
# o.set_config('drift:wind_uncertainty', 1)

o.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)
o.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)
o.plot(fast=False, background=['x_sea_water_velocity', 'y_sea_water_velocity'], buffer = .02)
o.plot(fast = False, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 150, linewidth = 3, title = "", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r"Wind speed $\mathrm{(m.s^{-1})}$", cpad = 0.08, text=text)

#%%
# Measurement modified model simulation
Expand All @@ -124,8 +138,7 @@
o2.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)
o2.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)

o.plot(fast=False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst onlyt', 'Gaussian measurment'], buffer = .02, markersize = 70, linewidth = 1)

o.plot(fast = False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 150, linewidth = 3, title = "", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r"Wind speed $\mathrm{(m.s^{-1})}$", cpad = 0.08, text=text)

#%%
# Here, we generate more particles and look
Expand All @@ -146,7 +159,7 @@
o2.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)
o2.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)

o.plot(fast=False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst onlyt', 'Gaussian measurment'], buffer = .02, markersize = 70, linewidth = 1)
o.plot(fast = False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 70, linewidth = 1, title = "", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r"Wind speed $\mathrm{(m.s^{-1})}$", cpad = 0.08, text=text)

lon_drifters_1, lat_drifters_1 = o.get_lonlats()
lon_drifters_2, lat_drifters_2 = o2.get_lonlats()
Expand All @@ -158,3 +171,9 @@
plt.xlabel("Skillscore")
plt.ylabel("Histogram (percentages)")
plt.show()

plt.figure()
plt.hist(skillscores, bins = 100, range = (0.9, 0.95), facecolor = 'none', edgecolor = 'C0')
plt.xlabel("Skillscore")
plt.ylabel("Histogram (percentages)")
plt.show()
22 changes: 15 additions & 7 deletions opendrift/models/basemodel/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2300,6 +2300,8 @@ def set_up_map(self,
lscale=None,
fast=False,
hide_landmask=False,
xlocs = None,
ylocs = None,
**kwargs):
"""
Generate Figure instance on which trajectories are plotted.
Expand Down Expand Up @@ -2455,7 +2457,7 @@ def set_up_map(self,
latmax, fast, ocean_color,
land_color, lscale, globe)

gl = ax.gridlines(ccrs.PlateCarree(globe=globe), draw_labels=True)
gl = ax.gridlines(ccrs.PlateCarree(globe=globe), draw_labels=True, xlocs = xlocs, ylocs = ylocs)
gl.top_labels = None

fig.canvas.draw()
Expand Down Expand Up @@ -2539,6 +2541,8 @@ def animation(self,
fast=False,
blit=False,
frames=None,
xlocs = None,
ylocs = None,
**kwargs):
"""Animate last run."""

Expand Down Expand Up @@ -2600,7 +2604,7 @@ def animation(self,
# Find map coordinates and plot points with empty data
fig, ax, crs, x, y, index_of_first, index_of_last = \
self.set_up_map(buffer=buffer, corners=corners, lscale=lscale,
fast=fast, hide_landmask=hide_landmask, **kwargs)
fast=fast, hide_landmask=hide_landmask, xlocs = xlocs, ylocs = ylocs, **kwargs)

gcrs = ccrs.PlateCarree(globe=crs.globe)

Expand Down Expand Up @@ -3268,6 +3272,9 @@ def plot(self,
lalpha=None,
bgalpha=1,
clabel=None,
cpad=.05,
caspect=30,
cshrink=.8,
surface_color=None,
submerged_color=None,
markersize=20,
Expand All @@ -3277,6 +3284,8 @@ def plot(self,
lscale=None,
fast=False,
hide_landmask=False,
xlocs = None,
ylocs = None,
**kwargs):
"""Basic built-in plotting function intended for developing/debugging.
Expand Down Expand Up @@ -3353,7 +3362,7 @@ def plot(self,
tlatmax)

fig, ax, crs, x, y, index_of_first, index_of_last = \
self.set_up_map(buffer=buffer, corners=corners, lscale=lscale, fast=fast, hide_landmask=hide_landmask, **kwargs)
self.set_up_map(buffer=buffer, corners=corners, lscale=lscale, fast=fast, hide_landmask=hide_landmask, xlocs = xlocs, ylocs = ylocs, **kwargs)

# x, y are longitude, latitude -> i.e. in a PlateCarree CRS
gcrs = ccrs.PlateCarree(globe=crs.globe)
Expand Down Expand Up @@ -3619,9 +3628,9 @@ def plot(self,
if mappable is not None and colorbar is True:
cb = fig.colorbar(mappable,
orientation='horizontal',
pad=.05,
aspect=30,
shrink=.8,
pad=cpad,
aspect=caspect,
shrink=cshrink,
drawedges=False)
# TODO: need better control of colorbar content
if clabel is not None:
Expand Down Expand Up @@ -3684,7 +3693,6 @@ def plot(self,
logger.warning(traceback.format_exc())

#plt.gca().tick_params(labelsize=14)

#fig.canvas.draw()
if filename is not None:
plt.savefig(filename)
Expand Down

0 comments on commit 38f858e

Please sign in to comment.