From 2f69512b648a0d156b699b6b43ef6318f7056c45 Mon Sep 17 00:00:00 2001 From: johannesro Date: Thu, 8 Oct 2020 08:18:40 +0000 Subject: [PATCH 01/11] update to FVCOM reader --- opendrift/readers/reader_FVCOM.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opendrift/readers/reader_FVCOM.py b/opendrift/readers/reader_FVCOM.py index f67a0da71..642afe323 100644 --- a/opendrift/readers/reader_FVCOM.py +++ b/opendrift/readers/reader_FVCOM.py @@ -24,9 +24,9 @@ import numpy as np from netCDF4 import Dataset, MFDataset, num2date from scipy.interpolate import LinearNDInterpolator +import pyproj -from basereader import BaseReader, pyproj - +from opendrift.readers.basereader import BaseReader class Reader(BaseReader): From 913f8562454630bb242d3a8a58ddacab67a901a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Fri, 23 Oct 2020 15:02:26 +0200 Subject: [PATCH 02/11] fix bug in lcs computation (timing of forward flowmap) and function for returning CC-eigenvector --- examples/example_double_gyre_LCS_particles.py | 5 +- examples/example_double_gyre_LCS_snapshot.py | 65 +++++++++++++++++++ opendrift/models/basemodel.py | 4 +- opendrift/models/physics_methods.py | 37 ++++++++++- 4 files changed, 105 insertions(+), 6 deletions(-) create mode 100644 examples/example_double_gyre_LCS_snapshot.py diff --git a/examples/example_double_gyre_LCS_particles.py b/examples/example_double_gyre_LCS_particles.py index 99ac6dc41..3108003ca 100755 --- a/examples/example_double_gyre_LCS_particles.py +++ b/examples/example_double_gyre_LCS_particles.py @@ -3,7 +3,7 @@ Double gyre - LCS with particles ============================================ -Drift of particles in an idealised (analytical) eddy current field, +Drift of particles in an idealised (analytical) eddy current field, plotted on top of the LCS. This takes some minutes to calculate. """ @@ -53,8 +53,7 @@ o.disable_vertical_motion() o.run(duration=duration, time_step=time_step, time_step_output=time_step_output) -o.animation(buffer=0, lcs=lcs, hide_landmask=True) +o.animation(buffer=0, lcs=lcs, cmap='cividis', hide_landmask=True) #%% # .. image:: /gallery/animations/example_double_gyre_LCS_particles_0.gif - diff --git a/examples/example_double_gyre_LCS_snapshot.py b/examples/example_double_gyre_LCS_snapshot.py new file mode 100644 index 000000000..6801a0cd1 --- /dev/null +++ b/examples/example_double_gyre_LCS_snapshot.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python +""" +Double gyre - LCS with particles +============================================ + +Drift of particles in an idealised (analytical) eddy current field, +plotted on top of the LCS. This takes some minutes to calculate. +""" + +from datetime import datetime, timedelta +import matplotlib.pyplot as plt +import numpy as np + +from opendrift.readers import reader_double_gyre +from opendrift.models.oceandrift import OceanDrift + +#%% +# Setting some parameters +plot_time = timedelta(seconds=12) +duration = timedelta(seconds=6) # T +time_step=timedelta(seconds=.5) +time_step_output=timedelta(seconds=.5) +delta=.02 # spatial resolution +#steps = int(duration.total_seconds()/ time_step_output.total_seconds() + 1) + +o = OceanDrift(loglevel=20) + +#%% +# Note that Runge-Kutta here makes a difference to Euler scheme +o.set_config('drift:scheme', 'runge-kutta4') +o.disable_vertical_motion() +o.fallback_values['land_binary_mask'] = 0 + +double_gyre = reader_double_gyre.Reader(epsilon=.25, omega=0.628, A=0.1) +print(double_gyre) +o.add_reader(double_gyre) + + +#%% +# Calculate Lyapunov exponents +#times = [double_gyre.initial_time + n*time_step_output for n in range(steps)] +ftle = o.calculate_ftle(time=double_gyre.initial_time + plot_time, time_step=time_step, + duration=duration, delta=delta, RLCS=False) + +#%% +# Make run with particles for the same period +o.reset() +x = np.linspace(0.1,1.9,1000) +y = np.linspace(0.1,0.9,1000) +#X,Y = np.meshgrid(x,y) +lon, lat = double_gyre.xy2lonlat(x,y) + +o.seed_elements(lon, lat, + time=double_gyre.initial_time) +o.disable_vertical_motion() +o.run(duration=plot_time, time_step=time_step, + time_step_output=time_step_output) +lonmin, latmin = double_gyre.xy2lonlat(0.,0.) +lonmax, latmax = double_gyre.xy2lonlat(2.,1.) +o.plot(lcs=ftle, show_initial=False, linewidth = 0, corners =[lonmin, lonmax, latmin, latmax], cmap='cividis') + +#o.animation(buffer=0, lcs=ftle, hide_landmask=True) + +#%% +# .. image:: /gallery/animations/example_double_gyre_LCS_particles_0.gif diff --git a/opendrift/models/basemodel.py b/opendrift/models/basemodel.py index 69fcf2945..36755778b 100644 --- a/opendrift/models/basemodel.py +++ b/opendrift/models/basemodel.py @@ -4002,8 +4002,10 @@ def calculate_ftle(self, reader=None, delta=None, domain=None, # Backwards if ALCS is True: self.reset() + #self.seed_elements(lons.ravel(), lats.ravel(), + # time=t+duration, z=z) self.seed_elements(lons.ravel(), lats.ravel(), - time=t+duration, z=z) + time=t, z=z) self.run(duration=duration, time_step=-time_step) b_x1, b_y1 = proj( self.history['lon'].T[-1].reshape(X.shape), diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py index 85d073476..b96621210 100644 --- a/opendrift/models/physics_methods.py +++ b/opendrift/models/physics_methods.py @@ -25,7 +25,7 @@ def oil_wave_entrainment_rate_li2017(dynamic_viscosity, oil_density, interfacial # Z. Li, M.L. Spaulding, D. French McCay, J. Mar. Pollut. Bull. (2016): # An algorithm for modeling entrainment and naturally and chemically dispersed # oil droplet size distribution under surface breaking wave conditions - + if wave_breaking_fraction is None: if wind_speed is None: raise ValueError('wave_breaking_fraction or wind_speed must be provided') @@ -216,6 +216,39 @@ def ftle(X, Y, delta, duration): return FTLE +def cauchygreen_eigvals(X, Y, delta, duration): + """Calculate eigenvector and eigenvalues of cauchy-green strain tensor + returns eigenvalues, eigenvectors of the flow map X,Y + """ + + # From Johannes Rohrs + nx = X.shape[0] + ny = X.shape[1] + J = np.empty([nx,ny,2,2],np.float) + lamba = np.empty([nx,ny,2],np.float) + xi = np.empty([nx,ny,2,2],np.float) + + # gradient + dx = np.gradient(X) + dy = np.gradient(Y) + + # Jacobian + J[:,:,0,0] = dx[0] / (2*delta) + J[:,:,1,0] = dy[0] / (2*delta) + J[:,:,0,1] = dx[1] / (2*delta) + J[:,:,1,1] = dy[1] / (2*delta) + + for i in range(0,nx): + for j in range(0,ny): + # Green-Cauchy tensor + D = np.dot(np.transpose(J[i,j]), J[i,j]) + # its largest eigenvalue + lamda[i,j], xi[i,j] = np.linalg.eigh(D) + + return lamba, xi + + + class PhysicsMethods(object): """Physics methods to be inherited by OpenDriftSimulation class""" @@ -341,7 +374,7 @@ def advect_with_sea_ice(self, factor=1): ice_velocity_y = self.environment.y_sea_water_velocity + \ 0.015*self.environment.y_wind self.update_positions( - factor*ice_velocity_x, + factor*ice_velocity_x, factor*ice_velocity_y) def advect_wind(self, factor=1): From b91d298e4593414d410867733a931193583bb8d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Fri, 23 Oct 2020 19:45:23 +0200 Subject: [PATCH 03/11] developing more lcs tools --- examples/example_double_gyre_LCS_snapshot.py | 22 +++++-- opendrift/models/basemodel.py | 69 ++++++++++++++++++++ opendrift/models/physics_methods.py | 8 +-- 3 files changed, 89 insertions(+), 10 deletions(-) diff --git a/examples/example_double_gyre_LCS_snapshot.py b/examples/example_double_gyre_LCS_snapshot.py index 6801a0cd1..8a8a3240a 100644 --- a/examples/example_double_gyre_LCS_snapshot.py +++ b/examples/example_double_gyre_LCS_snapshot.py @@ -17,10 +17,10 @@ #%% # Setting some parameters plot_time = timedelta(seconds=12) -duration = timedelta(seconds=6) # T +duration = timedelta(seconds=12) # T time_step=timedelta(seconds=.5) time_step_output=timedelta(seconds=.5) -delta=.02 # spatial resolution +delta=.01 # spatial resolution #steps = int(duration.total_seconds()/ time_step_output.total_seconds() + 1) o = OceanDrift(loglevel=20) @@ -42,13 +42,15 @@ ftle = o.calculate_ftle(time=double_gyre.initial_time + plot_time, time_step=time_step, duration=duration, delta=delta, RLCS=False) +lcs = o.calculate_lcs(time=double_gyre.initial_time + plot_time, time_step=-time_step, + duration=duration, delta=delta) #%% # Make run with particles for the same period o.reset() -x = np.linspace(0.1,1.9,1000) -y = np.linspace(0.1,0.9,1000) -#X,Y = np.meshgrid(x,y) -lon, lat = double_gyre.xy2lonlat(x,y) +x = np.linspace(0.1,1.9,100) +y = np.linspace(0.1,0.9,100) +X,Y = np.meshgrid(x.ravel(),y.ravel()) +lon, lat = double_gyre.xy2lonlat(X,Y) o.seed_elements(lon, lat, time=double_gyre.initial_time) @@ -59,6 +61,14 @@ lonmax, latmax = double_gyre.xy2lonlat(2.,1.) o.plot(lcs=ftle, show_initial=False, linewidth = 0, corners =[lonmin, lonmax, latmin, latmax], cmap='cividis') +fig = plt.figure() +fig.add_subplot(211) +plt.pcolormesh(np.log(np.sqrt(lcs['eigval'][0,:,:,0])), cmap='cividis'),plt.colorbar() +plt.quiver(lcs['eigvec'][0,:,:,0,0], lcs['eigvec'][0,:,:,0,1]) +fig.add_subplot(212) +plt.pcolormesh(np.log(np.sqrt(lcs['eigval'][0,:,:,1])), cmap='cividis'),plt.colorbar() +plt.quiver(lcs['eigvec'][0,:,:,0,0], lcs['eigvec'][0,:,:,0,1]) + #o.animation(buffer=0, lcs=ftle, hide_landmask=True) #%% diff --git a/opendrift/models/basemodel.py b/opendrift/models/basemodel.py index 36755778b..307a075da 100644 --- a/opendrift/models/basemodel.py +++ b/opendrift/models/basemodel.py @@ -4019,6 +4019,75 @@ def calculate_ftle(self, reader=None, delta=None, domain=None, return lcs + def calculate_lcs(self, reader=None, delta=None, domain=None, + time=None, time_step=None, duration=None, z=0, + forward=False): + + if reader is None: + self.logger.info('No reader provided, using first available:') + reader = list(self.readers.items())[0][1] + self.logger.info(reader.name) + if isinstance(reader, pyproj.Proj): + proj = reader + elif isinstance(reader,str): + proj = pyproj.Proj(reader) + else: + proj = reader.proj + + import scipy.ndimage as ndimage + from opendrift.models.physics_methods import cg_eigenvectors + + if not isinstance(duration, timedelta): + duration = timedelta(seconds=duration) + + if domain==None: + xs = np.arange(reader.xmin, reader.xmax, delta) + ys = np.arange(reader.ymin, reader.ymax, delta) + else: + xmin, xmax, ymin, ymax = domain + xs = np.arange(xmin, xmax, delta) + ys = np.arange(ymin, ymax, delta) + + X, Y = np.meshgrid(xs, ys) + lons, lats = proj(X, Y, inverse=True) + + if time is None: + time = reader.start_time + if not isinstance(time, list): + time = [time] + # dictionary to hold LCS calculation + lcs = {'time': time, 'lon': lons, 'lat':lats} + lcs['LCS'] = np.zeros((len(time), len(ys), len(xs))) + lcs['eigvec'] = np.zeros((len(time), len(ys), len(xs), 2, 2 )) + lcs['eigval'] = np.zeros((len(time), len(ys), len(xs), 2)) + + T = np.abs(duration.total_seconds()) + for i, t in enumerate(time): + self.logger.info('Calculating LCS for ' + str(t)) + # Forward or bachward flow map + dir = -1 + if forward==True: + dir = 1. + self.reset() + self.seed_elements(lons.ravel(), lats.ravel(), + time=t, z=z) + self.run(duration=duration, time_step=time_step) + x1, y1 = proj( + self.history['lon'].T[-1].reshape(X.shape), + self.history['lat'].T[-1].reshape(X.shape)) + lamba,xi = cg_eigenvectors(x1-X, y1-Y, delta, T) + lcs['eigval'][i,:,:,:] = lamba + lcs['eigvec'][i,:,:,:,:] = xi + + #lcs['RLCS'] = np.ma.masked_invalid(lcs['RLCS']) + #lcs['ALCS'] = np.ma.masked_invalid(lcs['ALCS']) + # Flipping ALCS left-right. Not sure why this is needed + #lcs['ALCS'] = lcs['ALCS'][:,::-1,::-1] + lcs['eigval'] = lcs['eigval'][:,::-1,::-1] + lcs['eigvec'] = lcs['eigvec'][:,::-1,::-1] + return lcs + + def center_of_gravity(self, onlysurface=False): """ calculate center of mass and variance of all elements diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py index b96621210..7a30182a9 100644 --- a/opendrift/models/physics_methods.py +++ b/opendrift/models/physics_methods.py @@ -189,7 +189,7 @@ def stokes_drift_profile_breivik(stokes_u_surface, stokes_v_surface, def ftle(X, Y, delta, duration): - """Calculate Finite Time Lyapunov Exponents""" + """Calculate Finite Time Lyapunov Exponents from flow map""" # From Johannes Rohrs nx = X.shape[0] ny = X.shape[1] @@ -216,9 +216,9 @@ def ftle(X, Y, delta, duration): return FTLE -def cauchygreen_eigvals(X, Y, delta, duration): +def cg_eigenvectors(X, Y, delta, duration): """Calculate eigenvector and eigenvalues of cauchy-green strain tensor - returns eigenvalues, eigenvectors of the flow map X,Y + returns eigenvalues, eigenvectors of the flow map X,Y """ # From Johannes Rohrs @@ -243,7 +243,7 @@ def cauchygreen_eigvals(X, Y, delta, duration): # Green-Cauchy tensor D = np.dot(np.transpose(J[i,j]), J[i,j]) # its largest eigenvalue - lamda[i,j], xi[i,j] = np.linalg.eigh(D) + lamba[i,j], xi[i,j] = np.linalg.eigh(D) return lamba, xi From 8bac376c36ee465dd3f0ed74e75a8611913be913 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Tue, 22 Dec 2020 13:33:35 +0100 Subject: [PATCH 04/11] FTLE example for norkyst working as expected --- examples/example_FTLE_norkyst.py | 43 ++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 examples/example_FTLE_norkyst.py diff --git a/examples/example_FTLE_norkyst.py b/examples/example_FTLE_norkyst.py new file mode 100644 index 000000000..b58f66941 --- /dev/null +++ b/examples/example_FTLE_norkyst.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +""" +LCS Norkyst +================================== +""" + +from datetime import datetime, timedelta +import numpy as np +from opendrift.readers import reader_netCDF_CF_generic +from opendrift.models.oceandrift import OceanDrift + +o = OceanDrift(loglevel=20) # Set loglevel to 0 for debug information + +# Norkyst ocean model +reader_norkyst = reader_netCDF_CF_generic.Reader(o.test_data_folder() + '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc') + +#o.add_reader([reader_norkyst, reader_arome]) +o.add_reader([reader_norkyst]) + + +#%% +# Calculating attracting/backwards FTLE at 20 hours +lcs = o.calculate_ftle( + time=reader_norkyst.start_time + timedelta(hours=24), + time_step=timedelta(minutes=15), + duration=timedelta(hours=3), delta=800, + RLCS=False) + +#%% +# Simulation from beginning and up to 30 hours (time of LCS) +o.reset() +lons = np.linspace(3.2, 5.0, 100) +lats = np.linspace(59.8, 61, 100) +lons, lats = np.meshgrid(lons, lats) +lons = lons.ravel() +lats = lats.ravel() +o.seed_elements(lons, lats, radius=0, number=10000, + time=reader_norkyst.start_time) + +o.run(end_time=reader_norkyst.start_time+timedelta(hours=24), + time_step=timedelta(minutes=30)) + +o.plot(lcs=lcs, vmin=1e-7, vmax=1e-4, cmap='Reds', markersize=1, colorbar=True, show_particles=True, show_initial=False, linewidth=0) From 1dbb439f288bd755e244ffc4fc676d1d038faa1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Tue, 22 Dec 2020 13:35:51 +0100 Subject: [PATCH 05/11] example for LCS in norkyst now based on strainlines, not FTLE --- examples/example_LCS_norkyst.py | 106 +++++++++++++++++++++++++++++--- 1 file changed, 96 insertions(+), 10 deletions(-) diff --git a/examples/example_LCS_norkyst.py b/examples/example_LCS_norkyst.py index 06a1daad2..41e9ce2f1 100755 --- a/examples/example_LCS_norkyst.py +++ b/examples/example_LCS_norkyst.py @@ -5,9 +5,11 @@ """ from datetime import datetime, timedelta - from opendrift.readers import reader_netCDF_CF_generic from opendrift.models.oceandrift import OceanDrift +import matplotlib.pyplot as plt +import numpy as np +import cartopy.crs as ccrs o = OceanDrift(loglevel=20) # Set loglevel to 0 for debug information @@ -19,20 +21,104 @@ #%% -# Calculating attracting/backwards FTLE/LCS at 20 hours -lcs = o.calculate_ftle( - time=reader_norkyst.start_time + timedelta(hours=20), - time_step=timedelta(minutes=30), - duration=timedelta(hours=5), delta=1000, - RLCS=False) + +x0,y0 = reader_norkyst.lonlat2xy(4.5,60.5) +d = 30000 +lcsproj = reader_norkyst.proj +lcs = o.calculate_lcs( + time=reader_norkyst.start_time + timedelta(hours=24), + time_step=timedelta(minutes=15), reader=lcsproj, + duration=timedelta(hours=3), delta=400, domain=[x0-d,x0+d,y0-d,y0+d]) + +l1 = lcs['eigval'][0,:,:,0] +l2 = lcs['eigval'][0,:,:,1] +mask = l2>l1 + +lmax = l1.copy() +lmax[mask] = l2[mask] + +lmin = l2.copy() +lmin[mask] = l1[mask] + +xi1 = lcs['eigvec'][0,:,:,0] +xi2 = lcs['eigvec'][0,:,:,1] + +stretch = xi1.copy() +stretch[mask] = xi2[mask] + +shrink = xi2.copy() +shrink[mask] = xi1[mask] + + +# find local maxima of largest eigenvalue +from skimage.feature import peak_local_max +#peaks = peak_local_max(lcs['eigval'][0,:,:,1],5) +peaks = peak_local_max(-lmin,5) +plon = [lcs['lon'][peaks[i,0],peaks[i,1]] for i in range(peaks.shape[0])] +plat = [lcs['lat'][peaks[i,0],peaks[i,1]] for i in range(peaks.shape[0])] + + +from opendrift.readers.reader_constant_2d import Reader +x,y = reader_norkyst.lonlat2xy(lcs['lon'],lcs['lat']) +r = Reader(x=x, y=y, proj4=reader_norkyst.proj4, array_dict = {'x_sea_water_velocity': stretch[:,:,0], 'y_sea_water_velocity': stretch[:,:,1]}) + + +xi = OceanDrift(loglevel=20) +xi.add_reader(r) +xi.seed_elements(plon, plat, time=datetime.now()) +xi.run(duration=timedelta(hours=3), time_step=300) +str_lon = xi.history['lon'] +str_lat = xi.history['lat'] + +#xi.plot(linewidth=2, linecolor='r',show_particles=False) + +xi.reset() +xi.seed_elements(plon, plat, time=datetime.now()) +xi.run(duration=timedelta(hours=3), time_step=-300) +str_lon = np.concatenate((str_lon, xi.history['lon'][:,::-1]), axis=0) +str_lat = np.concatenate((str_lat, xi.history['lat'][:,::-1]), axis=0) + #%% # Simulation from beginning and up to 30 hours (time of LCS) + o.reset() -o.seed_elements(lon=4.4, lat=60.2, number=1000, radius=1000, +lons = np.linspace(4.0, 5.0, 100) +lats = np.linspace(60., 61., 100) +lons, lats = np.meshgrid(lons, lats) +lons = lons.ravel() +lats = lats.ravel() +o.seed_elements(lons, lats, radius=0, number=10000, time=reader_norkyst.start_time) -o.run(end_time=reader_norkyst.start_time+timedelta(hours=20), + +o.run(end_time=reader_norkyst.start_time+timedelta(hours=24), time_step=timedelta(minutes=30)) +ax, plt = o.plot(cmap='Reds', vmax=2, markersize=1, colorbar=True, show_particles=True, show_initial=False, linewidth=0, show=False) + +''' +fig = plt.figure() +#proj=reader_norkyst.proj + +lon, lat = lcs['lon'], lcs['lat'] +x,y = reader_norkyst.lonlat2xy(lon, lat) +px,py = reader_norkyst.lonlat2xy(plon, plat) +#strx, stry = reader_norkyst.lonlat2xy(str_lon, str_lat) +fig.add_subplot(121) +plt.pcolormesh(x, y, np.log(np.sqrt(lmin)),vmin=-1,vmax=3, cmap='cividis'),plt.colorbar() +plt.quiver(x[::5,::5], y[::5,::5], shrink[::5,::5,0], shrink[::5,::5,1]) +plt.plot(px, py, '.r') +plt.title('shrink lines') +fig.add_subplot(122) + +plt.pcolormesh(x,y , np.log(np.sqrt(lmax)), cmap='cividis'),plt.colorbar() +plt.quiver(x[::5,::5], y[::5,::5], stretch[::5,::5,0], stretch[::5,::5,1]) + +plt.title('stretch lines') -o.plot(lcs=lcs, vmin=1e-7, vmax=1e-4, colorbar=True, show_particles=True) +plt.figure() +''' +gcrs = ccrs.PlateCarree() +ax.plot(str_lon.T, str_lat.T, '-r', ms=0.5, transform=gcrs) +plt.show() +#o.animation(buffer=0, lcs=ftle, hide_landmask=True) From eadc7e9ef44ee3199312dcc2f118426186f2c021 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Tue, 22 Dec 2020 13:38:54 +0100 Subject: [PATCH 06/11] display one frame of double gyre for testing algorithm --- examples/example_double_gyre_LCS_snapshot.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/example_double_gyre_LCS_snapshot.py b/examples/example_double_gyre_LCS_snapshot.py index 8a8a3240a..59039611a 100644 --- a/examples/example_double_gyre_LCS_snapshot.py +++ b/examples/example_double_gyre_LCS_snapshot.py @@ -47,13 +47,13 @@ #%% # Make run with particles for the same period o.reset() -x = np.linspace(0.1,1.9,100) -y = np.linspace(0.1,0.9,100) -X,Y = np.meshgrid(x.ravel(),y.ravel()) -lon, lat = double_gyre.xy2lonlat(X,Y) +x = [.9] +y = [.5] +lon, lat = double_gyre.xy2lonlat(x, y) -o.seed_elements(lon, lat, +o.seed_elements(lon, lat, radius=.15, number=2000, time=double_gyre.initial_time) + o.disable_vertical_motion() o.run(duration=plot_time, time_step=time_step, time_step_output=time_step_output) @@ -68,7 +68,7 @@ fig.add_subplot(212) plt.pcolormesh(np.log(np.sqrt(lcs['eigval'][0,:,:,1])), cmap='cividis'),plt.colorbar() plt.quiver(lcs['eigvec'][0,:,:,0,0], lcs['eigvec'][0,:,:,0,1]) - +plt.title('Eigenvalues and Eigenvectors of Cauchy-Green Strain tensor') #o.animation(buffer=0, lcs=ftle, hide_landmask=True) #%% From 9363cdd6903050481784239980f117aae73baac4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B8hrs?= Date: Tue, 22 Dec 2020 13:39:42 +0100 Subject: [PATCH 07/11] return eigenvalue and eigenvector in LCS method --- opendrift/models/basemodel.py | 41 +++++++++++++---------------- opendrift/models/physics_methods.py | 4 +-- 2 files changed, 21 insertions(+), 24 deletions(-) diff --git a/opendrift/models/basemodel.py b/opendrift/models/basemodel.py index 3c6e1a4c6..e1d1b7cae 100644 --- a/opendrift/models/basemodel.py +++ b/opendrift/models/basemodel.py @@ -1638,7 +1638,7 @@ def seed_repeated_segment(self, lons, lats, 2) total_number, in which case the number of elements per segment is: total_number / len(times). Any extra elements are duplicated along at the first segment. - + """ numtimes = int((end_time-start_time).total_seconds()/ @@ -1711,11 +1711,11 @@ def seed_within_polygon(self, lons, lats, number, **kwargs): area = 0.0 for i in range(-1, len(x)-1): area += x[i] * (y[i+1] - y[i-1]) - area = abs(area) / 2 # Make points, evenly distributed deltax = np.sqrt(area/number) lonpoints = np.array([]) + area = abs(area) / 2 latpoints = np.array([]) lonlat = poly.get_xy() lon = lonlat[:, 0] @@ -1990,7 +1990,7 @@ def remove_deactivated_elements(self): def run(self, time_step=None, steps=None, time_step_output=None, duration=None, end_time=None, outfile=None, export_variables=None, - export_buffer_length=100, stop_on_error=False): + export_buffer_length=100, stop_on_error=False, move_from_land=True): """Start a trajectory simulation, after initial configuration. Performs the main loop: @@ -2233,14 +2233,15 @@ def run(self, time_step=None, steps=None, time_step_output=None, self.timer_end('preparing main loop:making dynamical landmask') # Move point seed on land to ocean - if self.get_config('seed:ocean_only') is True and \ - ('land_binary_mask' not in self.fallback_values) and \ - ('land_binary_mask' in self.required_variables): - self.timer_start('preparing main loop:moving elements to ocean') - self.elements_scheduled.lon, self.elements_scheduled.lat = \ + if move_from_land==True: + if self.get_config('seed:ocean_only') is True and \ + ('land_binary_mask' not in self.fallback_values) and \ + ('land_binary_mask' in self.required_variables): + self.timer_start('preparing main loop:moving elements to ocean') + self.elements_scheduled.lon, self.elements_scheduled.lat = \ self.closest_ocean_points(self.elements_scheduled.lon, self.elements_scheduled.lat) - self.timer_end('preparing main loop:moving elements to ocean') + self.timer_end('preparing main loop:moving elements to ocean') #################################################################### # Preparing history array for storage in memory and eventually file @@ -2840,7 +2841,7 @@ def animation(self, buffer=.2, corners=None, filename=None, compare=None, self.get_density_array(pixelsize_m=density_pixelsize_m, weight=density_weight) H = H + H_submerged + H_stranded - + # x, y are longitude, latitude -> i.e. in a PlateCarree CRS gcrs = ccrs.PlateCarree() @@ -4211,8 +4212,7 @@ def calculate_ftle(self, reader=None, delta=None, domain=None, return lcs def calculate_lcs(self, reader=None, delta=None, domain=None, - time=None, time_step=None, duration=None, z=0, - forward=False): + time=None, time_step=None, duration=None, z=0): if reader is None: self.logger.info('No reader provided, using first available:') @@ -4256,13 +4256,10 @@ def calculate_lcs(self, reader=None, delta=None, domain=None, for i, t in enumerate(time): self.logger.info('Calculating LCS for ' + str(t)) # Forward or bachward flow map - dir = -1 - if forward==True: - dir = 1. self.reset() self.seed_elements(lons.ravel(), lats.ravel(), time=t, z=z) - self.run(duration=duration, time_step=time_step) + self.run(duration=duration, time_step=time_step, move_from_land=False) x1, y1 = proj( self.history['lon'].T[-1].reshape(X.shape), self.history['lat'].T[-1].reshape(X.shape)) @@ -4270,12 +4267,12 @@ def calculate_lcs(self, reader=None, delta=None, domain=None, lcs['eigval'][i,:,:,:] = lamba lcs['eigvec'][i,:,:,:,:] = xi - #lcs['RLCS'] = np.ma.masked_invalid(lcs['RLCS']) - #lcs['ALCS'] = np.ma.masked_invalid(lcs['ALCS']) - # Flipping ALCS left-right. Not sure why this is needed - #lcs['ALCS'] = lcs['ALCS'][:,::-1,::-1] - lcs['eigval'] = lcs['eigval'][:,::-1,::-1] - lcs['eigvec'] = lcs['eigvec'][:,::-1,::-1] + lcs['eigvec'] = np.ma.masked_invalid(lcs['eigvec']) + lcs['eigval'] = np.ma.masked_invalid(lcs['eigval']) + # Flipping ALCS left-right if using backward time flow map; Not sure why this is needed + if time_step.total_seconds() < 0: + lcs['eigval'] = lcs['eigval'][:,::-1,::-1] + lcs['eigvec'] = lcs['eigvec'][:,::-1,::-1] return lcs diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py index 7a30182a9..5f256bf39 100644 --- a/opendrift/models/physics_methods.py +++ b/opendrift/models/physics_methods.py @@ -243,8 +243,8 @@ def cg_eigenvectors(X, Y, delta, duration): # Green-Cauchy tensor D = np.dot(np.transpose(J[i,j]), J[i,j]) # its largest eigenvalue - lamba[i,j], xi[i,j] = np.linalg.eigh(D) - + lamba[i,j], xi[i,j] = np.linalg.eig(D) + lamba[lamba>1000]=np.nan # mask out land artifacts return lamba, xi From e43234d27f45145c0db6a120f82768081f52a45b Mon Sep 17 00:00:00 2001 From: vilhelm_fagerstrom Date: Thu, 18 Nov 2021 17:56:21 +0100 Subject: [PATCH 08/11] new module for L. Pertusa --- opendrift/models/lophelia.py | 215 +++++++++++++++++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 opendrift/models/lophelia.py diff --git a/opendrift/models/lophelia.py b/opendrift/models/lophelia.py new file mode 100644 index 000000000..28e13d0bb --- /dev/null +++ b/opendrift/models/lophelia.py @@ -0,0 +1,215 @@ +# This file is part of OpenDrift. +# +# OpenDrift is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, version 2 +# +# OpenDrift is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with OpenDrift. If not, see . +# +# Copyright 2015, Knut-Frode Dagestad, MET Norway + +import numpy as np +import logging; logger = logging.getLogger(__name__) + +from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray +#from opendrift.elements import LagrangianArray + + +# Defining element properties +class LopheliaLarvae(Lagrangian3DArray): + """Extending Lagrangian3DArray with specific properties for Lophelia larvae + """ + + variables = Lagrangian3DArray.add_variables([ + ('diameter', {'dtype': np.float32, + 'units': 'm', + 'default': 0.0002}), # For Lophelia + ('density', {'dtype': np.float32, + 'units': 'kg/m^3', + 'default': 1028.}), + ('age_seconds', {'dtype': np.float32, + 'units': 's', + 'default': 0.}), + ('competence', {'dtype': np.float32, + 'units': 's', + 'default': 0.}), + ('hatched', {'dtype': np.float32, + 'units': '', + 'default': 0.})]) + +class LopheliaLarvaeDrift(OceanDrift): + """Buoyant particle trajectory model based on the OpenDrift framework. + + Developed at MET Norway + + Generic module for particles that are subject to vertical turbulent + mixing with the possibility for positive or negative buoyancy + + Lophelia larvae assumed to ber neutrally buoyant throughout + + Particles could be e.g. oil droplets, plankton, or sediments + + Under construction. + """ + + + + ElementType = LopheliaLarvae + + required_variables = { + 'x_sea_water_velocity': {'fallback': 0}, + 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0}, + 'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0}, + 'sea_surface_wave_significant_height': {'fallback': 0}, + 'sea_ice_area_fraction': {'fallback': 0}, + 'x_wind': {'fallback': 0}, + 'y_wind': {'fallback': 0}, + 'land_binary_mask': {'fallback': None}, + 'sea_floor_depth_below_sea_level': {'fallback': 100}, + 'ocean_vertical_diffusivity': {'fallback': 0.02, 'profiles': True}, + 'ocean_mixed_layer_thickness': {'fallback': 50}, + 'sea_water_temperature': {'fallback': 10, 'profiles': True}, + 'sea_water_salinity': {'fallback': 34, 'profiles': True}, + 'surface_downward_x_stress': {'fallback': 0}, + 'surface_downward_y_stress': {'fallback': 0}, + 'turbulent_kinetic_energy': {'fallback': 0}, + 'turbulent_generic_length_scale': {'fallback': 0}, + 'upward_sea_water_velocity': {'fallback': 0}, + } + + # Vertical profiles of the following parameters will be available in + # dictionary self.environment.vertical_profiles + # E.g. self.environment_profiles['x_sea_water_velocity'] + # will be an array of size [vertical_levels, num_elements] + # The vertical levels are available as + # self.environment_profiles['z'] or + # self.environment_profiles['sigma'] (not yet implemented) + required_profiles = ['sea_water_temperature', + 'sea_water_salinity', + 'ocean_vertical_diffusivity'] + # The depth range (in m) which profiles shall cover + required_profiles_z_range = [-3000, 0] + + + # Default colors for plotting + status_colors = {'initial': 'green', 'active': 'blue', + 'hatched': 'red', 'eaten': 'yellow', 'died': 'magenta'} + + + def __init__(self, *args, **kwargs): + + # Calling general constructor of parent class + super(LopheliaLarvaeDrift, self).__init__(*args, **kwargs) + + # By default, elements do not strand towards coastline + self.set_config('general:coastline_action', 'previous') + + # Vertical mixing is enabled by default + self.set_config('drift:vertical_mixing', True) + + def update_terminal_velocity(self, Tprofiles=None, + Sprofiles=None, z_index=None): + """Calculate terminal velocity for larvae + + according to + S. Sundby (1983): A one-dimensional model for the vertical + distribution of pelagic fish eggs in the mixed layer + Deep Sea Research (30) pp. 645-661 + + Method copied from ibm.f90 module of LADIM: + Vikebo, F., S. Sundby, B. Aadlandsvik and O. Otteraa (2007), + Fish. Oceanogr. (16) pp. 216-228 + """ + g = 9.81 # ms-2 + + # Properties that determine element buoyancy + larvaesize = self.elements.diameter # 0.0002 for Lophelia + #eggsalinity = self.elements.neutral_buoyancy_salinity + + # Prepare interpolation of temperature and salinity + # if not (Tprofiles is None and Sprofiles is None): + # if z_index is None: + # z_i = range(Tprofiles.shape[0]) # Evtl. move out of loop + # z_index = interp1d(-self.environment_profiles['z'], + # z_i, bounds_error=False) + # zi = z_index(-self.elements.z) + # upper = np.maximum(np.floor(zi).astype(np.int), 0) + # lower = np.minimum(upper+1, Tprofiles.shape[0]-1) + # weight_upper = 1 - (zi - upper) + # + # # Do interpolation of temp, salt if profiles were passed into + # # this function, if not, use reader by calling self.environment + # if Tprofiles is None: + # T0 = self.environment.sea_water_temperature + # else: + # T0 = Tprofiles[upper, range(Tprofiles.shape[1])] * \ + # weight_upper + \ + # Tprofiles[lower, range(Tprofiles.shape[1])] * \ + # (1-weight_upper) + # if Sprofiles is None: + # S0 = self.environment.sea_water_salinity + # else: + # S0 = Sprofiles[upper, range(Sprofiles.shape[1])] * \ + # weight_upper + \ + # Sprofiles[lower, range(Sprofiles.shape[1])] * \ + # (1-weight_upper) + # + # # The larvae have the same temperature and salinity as the ambient water (can be modified). + # DENSw = self.sea_water_density(T=T0, S=S0) + # DENSegg = self.sea_water_density(T=T0, S=S0) + # dr = DENSw-DENSegg # density difference + # + # W = np.zeros(len(self.elements.ID)) + + # Set vertical velocities depending on stage (age) in larval phase: + + scenario = 'short' # Set scenario, 'short' or 'long' depending on length of pre-competency period + sday = 24*60*60 # Seconds in one day + reftemp = 8 # degree C + + if scenario == 'short': + dup2dwn = 21 + elif scenario == 'long': + dup2dwn = 42 + + idd = np.arange(15, dtype=float) + idd = np.insert(idd, [7, 8], [6.001, 7.001]) + idd = np.append(idd, [dup2dwn, dup2dwn+1]) + idd = idd*sday*reftemp + + w_idd = np.array([0, 0, 0.00005, 0.00005, 0.00005, 0.0002, 0.00025, 0, 0, 0, 0.0003, 0.00035, 0.0004, 0.00045, 0.0005, 0.00055, 0.0006, 0.0006, -0.0006]) # Vertical velocities at specified stages + + W = np.interp(self.elements.competence, idd, w_idd) + + + self.elements.terminal_velocity = W + + def update(self): + """Update positions and properties of particles.""" + + # Stokes drift + self.stokes_drift() + + # Update element age + self.elements.age_seconds += self.time_step.total_seconds() + + # Update competence + self.elements.competence += self.time_step.total_seconds() * self.environment.sea_water_temperature + + # Turbulent Mixing + self.update_terminal_velocity() + self.vertical_mixing() + + # Horizontal advection + self.advect_ocean_current() + + # Vertical advection + if self.get_config('drift:vertical_advection') is True: + self.vertical_advection() From caebf1253693ba0da91ea8607306868456a1fd6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B6hrs?= Date: Fri, 19 Nov 2021 11:48:31 +0100 Subject: [PATCH 09/11] Update lophelia.py Include deativation of larvae due to spawning, settling and death --- opendrift/models/lophelia.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/opendrift/models/lophelia.py b/opendrift/models/lophelia.py index 28e13d0bb..f8083fa6e 100644 --- a/opendrift/models/lophelia.py +++ b/opendrift/models/lophelia.py @@ -100,7 +100,7 @@ class LopheliaLarvaeDrift(OceanDrift): # Default colors for plotting status_colors = {'initial': 'green', 'active': 'blue', - 'hatched': 'red', 'eaten': 'yellow', 'died': 'magenta'} + 'spawned': 'red', 'settled': 'yellow', 'died': 'magenta'} def __init__(self, *args, **kwargs): @@ -172,7 +172,7 @@ def update_terminal_velocity(self, Tprofiles=None, scenario = 'short' # Set scenario, 'short' or 'long' depending on length of pre-competency period sday = 24*60*60 # Seconds in one day - reftemp = 8 # degree C + self.reftemp = 8 # degree C if scenario == 'short': dup2dwn = 21 @@ -182,7 +182,7 @@ def update_terminal_velocity(self, Tprofiles=None, idd = np.arange(15, dtype=float) idd = np.insert(idd, [7, 8], [6.001, 7.001]) idd = np.append(idd, [dup2dwn, dup2dwn+1]) - idd = idd*sday*reftemp + idd = idd*sday*self.reftemp w_idd = np.array([0, 0, 0.00005, 0.00005, 0.00005, 0.0002, 0.00025, 0, 0, 0, 0.0003, 0.00035, 0.0004, 0.00045, 0.0005, 0.00055, 0.0006, 0.0006, -0.0006]) # Vertical velocities at specified stages @@ -213,3 +213,11 @@ def update(self): # Vertical advection if self.get_config('drift:vertical_advection') is True: self.vertical_advection() + + if seft.time_step.total_seconds() < 0: + self.deactivate_elements(self.elements.competence <= 0., reason='spawned') + + if seft.time_step.total_seconds() > 0: + self.deactivate_elements(self.elements.competence > 21.*24*3600*ref_temp & self.elements.z < -sea_floor_depth, reason='settled') + self.deactivate_elements(self.elements.competence > 60.*24*3600* ref_temp , reason='died') + From e0bf2bbbb3f37c6a63e6003b4b0869adac61dce7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20R=C3=B6hrs?= Date: Fri, 19 Nov 2021 11:52:45 +0100 Subject: [PATCH 10/11] Update lophelia.py bugfix --- opendrift/models/lophelia.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opendrift/models/lophelia.py b/opendrift/models/lophelia.py index f8083fa6e..6535f0bd2 100644 --- a/opendrift/models/lophelia.py +++ b/opendrift/models/lophelia.py @@ -218,6 +218,6 @@ def update(self): self.deactivate_elements(self.elements.competence <= 0., reason='spawned') if seft.time_step.total_seconds() > 0: - self.deactivate_elements(self.elements.competence > 21.*24*3600*ref_temp & self.elements.z < -sea_floor_depth, reason='settled') - self.deactivate_elements(self.elements.competence > 60.*24*3600* ref_temp , reason='died') + self.deactivate_elements(self.elements.competence > 21.*24*3600* self.reftemp & self.elements.z < -sea_floor_depth, reason='settled') + self.deactivate_elements(self.elements.competence > 60.*24*3600* self.reftemp , reason='died') From 26dfad303e4ffcc772dc48613eb5f9337113f9de Mon Sep 17 00:00:00 2001 From: vilhelm_fagerstrom Date: Fri, 19 Nov 2021 14:58:43 +0100 Subject: [PATCH 11/11] update --- opendrift/models/lophelia.py | 63 ++++++------------------------------ 1 file changed, 10 insertions(+), 53 deletions(-) diff --git a/opendrift/models/lophelia.py b/opendrift/models/lophelia.py index 6535f0bd2..c5d7a5c8a 100644 --- a/opendrift/models/lophelia.py +++ b/opendrift/models/lophelia.py @@ -114,18 +114,9 @@ def __init__(self, *args, **kwargs): # Vertical mixing is enabled by default self.set_config('drift:vertical_mixing', True) - def update_terminal_velocity(self, Tprofiles=None, - Sprofiles=None, z_index=None): + def update_terminal_velocity(self, z_index=None): """Calculate terminal velocity for larvae - according to - S. Sundby (1983): A one-dimensional model for the vertical - distribution of pelagic fish eggs in the mixed layer - Deep Sea Research (30) pp. 645-661 - - Method copied from ibm.f90 module of LADIM: - Vikebo, F., S. Sundby, B. Aadlandsvik and O. Otteraa (2007), - Fish. Oceanogr. (16) pp. 216-228 """ g = 9.81 # ms-2 @@ -133,40 +124,6 @@ def update_terminal_velocity(self, Tprofiles=None, larvaesize = self.elements.diameter # 0.0002 for Lophelia #eggsalinity = self.elements.neutral_buoyancy_salinity - # Prepare interpolation of temperature and salinity - # if not (Tprofiles is None and Sprofiles is None): - # if z_index is None: - # z_i = range(Tprofiles.shape[0]) # Evtl. move out of loop - # z_index = interp1d(-self.environment_profiles['z'], - # z_i, bounds_error=False) - # zi = z_index(-self.elements.z) - # upper = np.maximum(np.floor(zi).astype(np.int), 0) - # lower = np.minimum(upper+1, Tprofiles.shape[0]-1) - # weight_upper = 1 - (zi - upper) - # - # # Do interpolation of temp, salt if profiles were passed into - # # this function, if not, use reader by calling self.environment - # if Tprofiles is None: - # T0 = self.environment.sea_water_temperature - # else: - # T0 = Tprofiles[upper, range(Tprofiles.shape[1])] * \ - # weight_upper + \ - # Tprofiles[lower, range(Tprofiles.shape[1])] * \ - # (1-weight_upper) - # if Sprofiles is None: - # S0 = self.environment.sea_water_salinity - # else: - # S0 = Sprofiles[upper, range(Sprofiles.shape[1])] * \ - # weight_upper + \ - # Sprofiles[lower, range(Sprofiles.shape[1])] * \ - # (1-weight_upper) - # - # # The larvae have the same temperature and salinity as the ambient water (can be modified). - # DENSw = self.sea_water_density(T=T0, S=S0) - # DENSegg = self.sea_water_density(T=T0, S=S0) - # dr = DENSw-DENSegg # density difference - # - # W = np.zeros(len(self.elements.ID)) # Set vertical velocities depending on stage (age) in larval phase: @@ -188,7 +145,7 @@ def update_terminal_velocity(self, Tprofiles=None, W = np.interp(self.elements.competence, idd, w_idd) - + print(w) self.elements.terminal_velocity = W def update(self): @@ -206,6 +163,7 @@ def update(self): # Turbulent Mixing self.update_terminal_velocity() self.vertical_mixing() + #self.vertical_buoyancy() # Horizontal advection self.advect_ocean_current() @@ -213,11 +171,10 @@ def update(self): # Vertical advection if self.get_config('drift:vertical_advection') is True: self.vertical_advection() - - if seft.time_step.total_seconds() < 0: - self.deactivate_elements(self.elements.competence <= 0., reason='spawned') - - if seft.time_step.total_seconds() > 0: - self.deactivate_elements(self.elements.competence > 21.*24*3600* self.reftemp & self.elements.z < -sea_floor_depth, reason='settled') - self.deactivate_elements(self.elements.competence > 60.*24*3600* self.reftemp , reason='died') - + + if self.time_step.total_seconds() < 0: + self.deactivate_elements(self.elements.competence <= 0., reason='spawned') + + if self.time_step.total_seconds() > 0: + self.deactivate_elements(self.elements.competence > 21.*24*3600* self.reftemp & self.elements.z < -sea_floor_depth, reason='settled') + self.deactivate_elements(self.elements.competence > 60.*24*3600* self.reftemp , reason='died')