Update to LCS algorithms

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)

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,109 @@
 
 
 #%%
-# 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')
+
+<<<<<<< HEAD
+plt.figure()
+'''
+
+gcrs = ccrs.PlateCarree()
+ax.plot(str_lon.T, str_lat.T, '-r', ms=0.5, transform=gcrs)
+plt.show()
 
 o.plot(lcs=lcs, vmin=1e-7, vmax=1e-4, colorbar=True, show_elements=True)
 
+
+#o.animation(buffer=0, lcs=ftle, hide_landmask=True)

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
-

examples/example_double_gyre_LCS_snapshot.py

@@ -0,0 +1,75 @@
+#!/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=12)  # T
+time_step=timedelta(seconds=.5)
+time_step_output=timedelta(seconds=.5)
+delta=.01  # 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)
+
+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 = [.9]
+y = [.5]
+lon, lat = double_gyre.xy2lonlat(x, y)
+
+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)
+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')
+
+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])
+plt.title('Eigenvalues and Eigenvectors of Cauchy-Green Strain tensor')
+#o.animation(buffer=0, lcs=ftle, hide_landmask=True)
+
+#%%
+# .. image:: /gallery/animations/example_double_gyre_LCS_particles_0.gif