Merge pull request #1628 from knutfrode/dev

[run-ex] Using TrajAn to calculate density in example_leeway, but plo…

Author: knutfrode

examples/example_leeway.py

@@ -7,6 +7,7 @@
 from datetime import timedelta
 import cmocean
 import xarray as xr
+import trajan as ta
 from opendrift import test_data_folder as tdf
 from opendrift.readers import reader_netCDF_CF_generic
 from opendrift.models.leeway import Leeway
@@ -36,7 +37,7 @@
 
 #%%
 # Running model
-o.run(duration=timedelta(hours=48), time_step=900, time_step_output=3600)
+ds = o.run(duration=timedelta(hours=48), time_step=1800, time_step_output=3600)
 
 #%%
 # Print and plot results
@@ -56,8 +57,9 @@
 
 #%%
 # Plot density of stranded elements
-d, dsub, dstr, lon, lat = o.get_density_array(pixelsize_m=3000)
-strand_density = xr.DataArray(dstr[-1,:,:], coords={'lon_bin': lon[0:-1], 'lat_bin': lat[0:-1]})
-o.plot(fast=True, background=strand_density.where(strand_density>0),
-       vmin=0, vmax=20, cmap=cmocean.cm.dense, clabel='Density of stranded elements',
+ds_stranded = ds.where(ds.status==1)
+grid = ds_stranded.traj.make_grid(dx=3000)
+ds_conc = ds_stranded.traj.concentration(grid).sum(dim='time')
+o.plot(fast=True, background=ds_conc.number.where(ds_conc.number>0),
+       cmap=cmocean.cm.thermal, clabel='Density of stranded elements',
        show_elements=False, linewidth=0)

opendrift/models/basemodel/__init__.py

@@ -3640,9 +3640,14 @@ def plot(self,
             else:
                 time = None
             if isinstance(background, xr.DataArray):
-                map_x = background.coords['lon_bin']
-                map_y = background.coords['lat_bin']
-                scalar = background
+                if 'lon' in background.coords:  # From TrajAn
+                    map_x = background.coords['lon']
+                    map_y = background.coords['lat']
+                    scalar = background.T
+                else:  # From get_density_array
+                    map_x = background.coords['lon_bin']
+                    map_y = background.coords['lat_bin']
+                    scalar = background
                 map_y, map_x = np.meshgrid(map_y, map_x)
             elif background == 'residence':
                 scalar, lon_res, lat_res = self.get_residence_time(