diff --git a/examples/example_coastline_options.py b/examples/example_coastline_options.py index 487cbff8b..6708331e2 100755 --- a/examples/example_coastline_options.py +++ b/examples/example_coastline_options.py @@ -5,34 +5,89 @@ Example to illustrate stranding options using an artificial east-west oscillating current field -Knut-Frode Dagestad, Feb 2017 """ -from datetime import datetime +from datetime import datetime, timedelta from opendrift.readers import reader_oscillating from opendrift.models.oceandrift import OceanDrift +duration = timedelta(hours=36) +time_step = timedelta(hours=.5) +number = 500 -o = OceanDrift(loglevel=50) # Set loglevel to 0 for debug information - +# Oscillating east-west current to illustrate stranding options reader_osc = reader_oscillating.Reader('x_sea_water_velocity', amplitude=1, zero_time=datetime.utcnow()) -o.add_reader([reader_osc]) # Oscillating east-west current component -o.set_config('environment:fallback:y_sea_water_velocity', .2) # Adding northwards drift +#%% +# Coastline option "stranding" +# ============================ +# +# Note that particles are "jumping" a distance onto land, depending on calculation timestep +o = OceanDrift(loglevel=50) +o.add_reader([reader_osc]) +o.set_config('environment:fallback:y_sea_water_velocity', .2) +o.set_config('general:coastline_action', 'stranding') +o.seed_elements(lon=12.2, lat=67.7, radius=5000, number=number, time=reader_osc.zero_time) +o.run(duration=duration, time_step=time_step) +print(f'Calculation time: {o.timing["total time"]}') +o.animation() #%% -# Try different options: 'previous', 'stranding', 'none' -o.set_config('general:coastline_action', 'previous') +# .. image:: /gallery/animations/example_coastline_options_0.gif -o.seed_elements(lon=12.2, lat=67.7, radius=5000, number=25, time=reader_osc.zero_time) +#%% +# Coastline option "stranding" with higher precision +# ================================================== +# +# by setting config "general:coastline_approximation_precision" to desired accuracy in degrees. +# Note that with a (too) large compuation time step, particles may still "jump" over islands. +# An alternative to avoid this possibility is to use a smaller timestep for the simulation, though at a larger computational cost. +o = OceanDrift(loglevel=50) +o.add_reader([reader_osc]) +o.set_config('environment:fallback:y_sea_water_velocity', .2) +o.set_config('general:coastline_action', 'stranding') +o.set_config('general:coastline_approximation_precision', .001) # approx 100m +o.seed_elements(lon=12.2, lat=67.7, radius=5000, number=number, time=reader_osc.zero_time) +o.run(duration=duration, time_step=time_step) +print(f'Calculation time: {o.timing["total time"]}') +o.animation() -o.run(steps=36*4, time_step=900, time_step_output=1800) +#%% +# .. image:: /gallery/animations/example_coastline_options_1.gif + +#% +# Coastline option "previous" +# =========================== +# +# Particles hitting land are moved back to last position in water, +# and may move offshore if currents/winds/forcing changes direction. +# Here we see that several particles are jumping over small island. +# Reducing timestep to e.g. 15 minutes will reduce this problem. +o = OceanDrift(loglevel=50) +o.add_reader([reader_osc]) +o.set_config('environment:fallback:y_sea_water_velocity', .2) +o.set_config('general:coastline_action', 'previous') +o.seed_elements(lon=12.2, lat=67.7, radius=5000, number=number, time=reader_osc.zero_time) +o.run(duration=duration, time_step=time_step) +print(f'Calculation time: {o.timing["total time"]}') +o.animation() #%% -# Print and plot results -print(o) +# .. image:: /gallery/animations/example_coastline_options_2.gif + +#% +# Coastline option "none" +# ======================= +# +# No interaction with land, as used by e.g. WindBlow model (atmospheric drift) +o = OceanDrift(loglevel=50) +o.add_reader([reader_osc]) +o.set_config('environment:fallback:y_sea_water_velocity', .2) +o.set_config('general:coastline_action', 'none') +o.seed_elements(lon=12.2, lat=67.7, radius=5000, number=number, time=reader_osc.zero_time) +o.run(duration=duration, time_step=time_step) +print(f'Calculation time: {o.timing["total time"]}') o.animation() -o.plot() #%% -# .. image:: /gallery/animations/example_coastline_options_0.gif +# .. image:: /gallery/animations/example_coastline_options_3.gif diff --git a/examples/example_fjord.py b/examples/example_fjord.py index 188318019..09e60d248 100755 --- a/examples/example_fjord.py +++ b/examples/example_fjord.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ Fjord -================================== +===== """ from datetime import timedelta @@ -40,4 +40,3 @@ # .. image:: /gallery/animations/example_fjord_0.gif o.plot() - diff --git a/examples/example_leeway_backtrack.py b/examples/example_long_leeway_backtrack.py similarity index 84% rename from examples/example_leeway_backtrack.py rename to examples/example_long_leeway_backtrack.py index 12795ba87..fc31efa56 100755 --- a/examples/example_leeway_backtrack.py +++ b/examples/example_long_leeway_backtrack.py @@ -40,6 +40,8 @@ # Define domain of possible origin lons = np.arange(3.4, 5, .1/20) lats = np.arange(59.7, 60.8, .05/20) +#lons = lons[0::2] # using every second, due to memory limitation on CircleCI +#lats = lats[0::2] corners = [lons[0], lons[-1], lats[0], lats[-1]] lons, lats = np.meshgrid(lons, lats) @@ -53,7 +55,8 @@ density_backwards = density_backwards.where(density_backwards>0) density_backwards = density_backwards/density_backwards.sum()*100 vmax = density_backwards.max() -o.plot(background=density_backwards, clabel='Probability of origin [%]', text=text, corners=corners, fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax) +o.plot(background=density_backwards, clabel='Probability of origin [%]', text=text, corners=corners, + fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax) os.remove(outfile) #%% @@ -61,7 +64,7 @@ o = Leeway(loglevel=50) o.add_reader([reader_norkyst, reader_arome]) o.seed_elements(lon=lons, lat=lats, radius=0, - time=start_time, object_type=object_type) + time=start_time, object_type=object_type) o.run(duration=duration, time_step=900, time_step_output=3600, outfile=outfile) print(o) @@ -78,16 +81,20 @@ hit_start_lats = lat[hits, 0] o_hit = opendrift.open(outfile, elements=hits) -o.animation(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], fast=True, corners=corners, text=text) +o.animation(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], + fast=True, corners=corners, text=text) #%% # .. image:: /gallery/animations/example_leeway_backtrack_0.gif -o.plot(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], show_elements=False, fast=True, corners=corners, text=text) +o.plot(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], + show_elements=False, fast=True, corners=corners, text=text) #%% # Plot the initial density of elements that actually hit the target after 24 hours. To be compared with the density figure from backwards simulation (see top) of = opendrift.open_xarray(outfile, elements=hits) density_forwards = of.get_histogram(pixelsize_m=5000).isel(time=0).isel(origin_marker=0) density_forwards = density_forwards.where(density_forwards>0) -o_hit.plot(background=density_forwards/density_forwards.sum()*100, clabel='Probability of origin [%]', text=text, corners=corners, fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax) +ratio = density_forwards/density_forwards.sum()*100 +o_hit.plot(background=ratio, clabel='Probability of origin [%]', text=text, corners=corners, + fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax) diff --git a/examples/example_radionuclides.py b/examples/example_long_radionuclides.py similarity index 100% rename from examples/example_radionuclides.py rename to examples/example_long_radionuclides.py diff --git a/examples/example_openberg_new.py b/examples/example_openberg.py similarity index 93% rename from examples/example_openberg_new.py rename to examples/example_openberg.py index ed7882bb8..1a90010f6 100755 --- a/examples/example_openberg_new.py +++ b/examples/example_openberg.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ -Icebergs (openberg) -==================== +Icebergs (OpenBerg module) +========================== """ from opendrift.models.openberg import OpenBerg diff --git a/examples/example_openberg_det.py b/examples/example_openberg_det.py deleted file mode 100755 index 43c917a9e..000000000 --- a/examples/example_openberg_det.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python -""" -Ice berg (openberg) deterministic -================================== -""" - -from opendrift.readers import reader_netCDF_CF_generic -from opendrift.models.openberg_old import OpenBergOld - -#%% -# Initialize model -steps = 60 # This is the number of forecast steps - -o = OpenBergOld() # Basic drift model suitable for icebergs - -#%% -# Preparing Readers -reader_current = reader_netCDF_CF_generic.Reader(o.test_data_folder() + - '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc') - -reader_wind = reader_netCDF_CF_generic.Reader(o.test_data_folder() + - '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc') - -o.add_reader([reader_current, reader_wind]) - -#%% -# Seeding elements -# -# Icebergs are moved with the ocean current as per Barker et al (2004), -# in addition to a fraction of the wind speed (wind_drift_factor). -# This factor depends on the properties of the elements. -# Default empirical values are: -# - Wind drift fraction: 0.018 (1.8 %) (Garret 1985) -# - Iceberg size: Keel dept = 60m -# Waterline length = 90.5m - - -o.seed_elements(3.3, 60., radius=3000, - time=reader_current.start_time, - water_line_length=100, - keel_depth=90, number=500) - - -print('Starting free run .../n') - -print('Start time: ' + str(o.start_time)) - -#%% -# Running model - -o.run(time_step=3600, steps=steps) - -#%% -# Print and plot results -o.plot(fast=True) -o.animation(fast=True) - - -#%% -# .. image:: /gallery/animations/example_openberg_det_0.gif - -print('############## Latitudes:', o.history['lat']) -print('############## Longitudes:',o.history['lon']) - diff --git a/examples/example_openberg_stat.py b/examples/example_openberg_stat.py deleted file mode 100755 index 95ee9b2cf..000000000 --- a/examples/example_openberg_stat.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python -""" -OpenBergOld - statistical mode -============================== -""" - -from datetime import datetime, timedelta -from opendrift.readers import reader_netCDF_CF_generic -from opendrift.readers import reader_current_from_track -from opendrift.models.openberg_old import OpenBergOld - -#%% -# Create observation of iceberg track -obslon = [3.1, 3.123456] -obslat = [61.1, 61.132198] -obstime = [datetime(2015, 11, 16, 0), datetime(2015, 11, 16, 6)] - -#%% -# Initialize model -steps = 60 # This is the number of forecast steps -o = OpenBergOld(loglevel=30) # Basic drift model suitable for icebergs - -#%% -# Preparing Readers -reader_wind = reader_netCDF_CF_generic.Reader(o.test_data_folder() + - '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc',name='WIND') - -reader_current = reader_current_from_track.Reader( - obslon, obslat, obstime, wind_east=0, wind_north=0, - windreader=reader_wind, wind_factor=0.018) - -o.add_reader([reader_current, reader_wind]) - -#%% -# Seeding elements -# -# Icebergs are moved with the ocean current as per Barker et al (2004), -# in addition to a fraction of the wind speed (wind_drift_factor). -# This factor depends on the properties of the elements. -# Default empirical values are: -# - Wind drift fraction: 0.018 (1.8 %) (Garret 1985) -# - Iceberg size: Keel dept = 60m -# Waterline length = 90.5m -# NB! Iceberg size is irrelevant for current_reader with 1D z-profile - -o.seed_elements(3.3, 61.3, radius=3000, number=500, - time=reader_current.start_time) - -#%% -# Run model -print('Starting free run .../n') - -print('Start time: ' + str(o.start_time)) - -o.run(time_step=3600, steps=steps) - -#%% -# Print and plot results -o.plot(fast=True) -o.animation(fast=True) - -#%% -# .. image:: /gallery/animations/example_openberg_stat_0.gif - diff --git a/examples/example_wind_measurements.py b/examples/example_wind_measurements.py index e926a73be..12bba07b8 100644 --- a/examples/example_wind_measurements.py +++ b/examples/example_wind_measurements.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ Using wind measurements -=============== +======================= """ from datetime import datetime, timedelta, date diff --git a/history.rst b/history.rst index 99c4570df..d2b594d4e 100644 --- a/history.rst +++ b/history.rst @@ -1,8 +1,13 @@ History ======= +Next release +------------ +* New feature by TheSylex to move stranded particles closer to the actual coastline with a precision given by config setting `general:coastline_approximation_precision` in unit of degrees (1deg approx 111 km) +* Major updates to OpenBerg iceberg drift model (Achref Othmani, Lenny Hucher) + 2024-06-27 / Release v1.11.10 ----------------------------- +----------------------------- * Using now standard env variable names COPERNICUSMARINE_SERVICE_USERNAME and COPERNICUSMARINE_SERVICE_PASSWORD for reader_copernicus. Env variable COPERNICUSMARINE_CACHE_DIRECTORY can be set to empty string to disable caching. 2024-06-27 / Release v1.11.9 @@ -11,6 +16,7 @@ History * Hack in generic reader to make sure wind from ECMWF files is at 10m height * Raising now error if all elements are seeded on land + 2024-06-26 / Release v1.11.8 ---------------------------- * Raising now error if all elements are seeded on land diff --git a/opendrift/models/leeway.py b/opendrift/models/leeway.py index 665f9b132..b79e90bd5 100644 --- a/opendrift/models/leeway.py +++ b/opendrift/models/leeway.py @@ -234,11 +234,8 @@ def __init__(self, d=None, *args, **kwargs): 'processes:capsizing': { 'type': 'bool', 'default': False, - 'min': 0, - 'max': 1, 'description': 'If True, elements can be capsized when wind exceeds threshold given by config item capsize:wind_threshold', - 'units': 'fraction', 'level': CONFIG_LEVEL_BASIC }, 'capsizing:leeway_fraction': {