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,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)

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

opendrift/models/basemodel.py

@@ -1753,11 +1753,11 @@ def seed_within_polygon(self, lons, lats, number=None, **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]
@@ -2043,7 +2043,7 @@ def set_fallback_values(self, refresh=False):
 
     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:
@@ -2294,14 +2294,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
@@ -4252,8 +4253,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),
@@ -4267,6 +4270,71 @@ 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):
+
+        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
+            self.reset()
+            self.seed_elements(lons.ravel(), lats.ravel(),
+                               time=t, z=z)
+            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))
+            lamba,xi =  cg_eigenvectors(x1-X, y1-Y, delta, T)
+            lcs['eigval'][i,:,:,:] = lamba
+            lcs['eigvec'][i,:,:,:,:] = xi
+
+        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
+
+
     def center_of_gravity(self, onlysurface=False):
         """
         calculate center of mass and variance of all elements