# API Reference

## models.plots

### Class

Open: load carbondrift simulations.

#### Class Methods

- __init__

                                  def __init__(self, filename):

                                      """Parameters
                                      -------------
                                      filename: str
                                          Name of .nc file to be imported."""

                                      self.data = Dataset(filename)
                                      self.steps_output = len(self.data["time"][:])

- get_property

                                  def get_property(self, property):
                                      return np.ma.transpose(self.data[property][:])

- get_time_array

                                  def get_time_array(self):
                                      time = self.data.variables["time"]
                                      units = time.units
                                      return num2date(time[:], units)

- get_properties_from_list

                                  def get_properties_from_list(self, properties):
                                      p = []
                                      for property in properties:
                                          if property in ["status", "origin_marker"]:
                                              prop = self.get_property(property)
                                              prop = np.ma.MaskedArray(prop.data, prop.mask, float)
                                              prop = np.ma.filled(prop, np.nan)
                                          else:
                                              prop = self.get_property(property)
                                              prop = np.ma.filled(prop, np.nan)
                                          p.append(prop)
                                      return p

Plot: plot carbondrift simulations.

#### Class Methods

- __init__

                                  def __init__(self, *files, cmap = None, color = None, linestyle = None,
                                      lons = None, lats = None, figsize = (20, 20), diffidx = None,
                                      fontsize = 17, title = None, depth = -200, add = False, suptitle = None,
                                      diff = True, absolute = False, fontweight = "normal", martincurve = None,
                                      outfile = None, shrink = 1, clip = None, locations = None, bins = None,
                                      prop1 = None, prop2 = None, colorbarlabel = None,
                                      xlabel = None, ylabel = None, xlim = None, ylim = None, linewidth = 2,
                                      legend = None, areagridpath = f"./supplementary_data/area_grid.npy",
                                      biomegridpath = f"./supplementary_data/biomegrid.npy", group = None,
                                      mhwintesitypath = None, dpi = 300):
                                      """
                                      Parameters
                                      -----------
                                      files: str
                                          Unlimited number of string simulation paths. At least one is required!
                                      cmap: str or None
                                          Matplolib cmap name (default None)
                                      color: list or None
                                          List of matplolib color names. E.g ["red","blue","orange"] (default None)
                                      linestyle: list or None
                                          List of linestyles. E.g. ["solid","dashed","dotted"] (default None)
                                      lons: numpy array or None
                                          Specific longitudes to use in plotting methods. If None the np.array(-180,180,1) is asummed.
                                      lons: numpy array or None
                                          Specific latitudes to use in plotting methods. If None the np.array(-90,90,1) is asummed.
                                      figsize: tuple
                                          Figure size (default(20, 20))
                                      diffidx: int or None
                                          Index of files which separates masses that are added with a plus sign and masses that are
                                          added with a minus sign. Only use if diff and add are both True. If None len(files) // 2 is asummed.
                                      fontsize: float
                                          Fontsize (default 17)
                                      title: str or None
                                          Figure title (default None)
                                      depth: float
                                          Depth in meters at which a calculation takes place. If positive a depth*(-1) preprocess is performed.
                                          If depth <= -5000, calculations at seafloor are asummed (default -200)
                                      add: bool
                                          Add file outputs (deafult False)
                                      suptitle: list or None
                                          List of suptitles (default None)
                                      diff: bool
                                          Take difference between two or more files if specified (default True)
                                      absolute: bool
                                          Show absolute value of mass/mass flux... (default False)
                                      fontweight: str
                                          Fontweight (default "normal")
                                      martincurve: float or None
                                          Martin curve coefficient. If specified the curve is plotted (default None)
                                      outfile: str or None
                                          Outfile path, if None plt.show() is called (default None)
                                      shrink: float
                                          Coorbar shrink (default 1)
                                      clip: list of length 2
                                          Len 2 list with min and max value for clipping (default None)
                                      locations: list of tuples or None
                                          Locations for dynamics analysis. Example [(longitude1, latitude1), (longitude2, latitude2),...]
                                          (default None)
                                      bins: int or None
                                          Number of bins for histogram (default None)
                                      prop1: str or None
                                          Property to plot on x axes. Must be the short name in the .nc file (default None)
                                      prop2: str or None
                                          Property to plot on y axes. Must be the short name in the .nc file (default None)
                                      xlabel: str or None
                                          xlabel (default None)
                                      ylabel: str or None
                                          ylabel (default None)
                                      xlim: list of length 2
                                          Len 2 list with min and max value for xlim (default None)
                                      ylim: list of length 2
                                          Len 2 list with min and max value for ylim (default None)
                                      linewidth: int
                                          linewidth (default 2)
                                      legend: list or None
                                          List of legend labels (default None)
                                      areagridpath: str
                                          Path to npy file containg area of each 1 by 1 grid cell (default f"./supplementary_data/area_grid.npy")
                                      biomegridpath: str
                                          Path to npy file containing biome int labels for each 1 by 1 grid cell (default f"./supplementary_data/biomegrid.npy")
                                      group: str or None
                                          Grouping of tracers. Choose between biome, lonmean, latmean or none (default None)
                                      mhwintensitypath: TODO EXPLAIN
                                      dpi: int
                                          Output image resolution (default 300)
                                      """

                                      logger.debug("Setting up figure.")
                                      plt.rcParams.update({'font.size': fontsize})
                                      fig, ax = plt.subplots(1, 1, figsize = figsize)
                                      self.fig = fig
                                      self.ax = ax

                                      if len(files) == 1 and (diff or add):
                                          logger.warning("Only one filepath given, changing diff and add to False.")
                                          diff = False
                                          add = False
                                      elif len(files) > 1 and not (diff or add):
                                          logger.warning("Multiple files are given but difference and add are set to False.")
                                      if diff and add and diffidx is None:
                                          logger.warning("Both diff and add are set to True, but diffidx is None. Setting diffidx to half the number of files.")
                                          diffidx = len(files) // 2

                                      logger.debug("Importing data.")

                                      self.objects =  []
                                      self.decay_type = []
                                      self.initial_velocities = []
                                      self.fragmentation = False
                                      for k in range(len(files)):
                                          obj = Open(files[k])
                                          setattr(self, "obj" + str(k + 1), obj)
                                          self.objects.append(getattr(self, "obj" + str(k + 1)))
                                          param_file = files[k].replace(".nc", ".json")
                                          if os.path.exists(param_file):
                                              with open(param_file, "r") as f:
                                                  param = json.load(f)
                                                  self.decay_type.append(param["decaytype"])
                                                  self.initial_velocities.append(param["initialvelocity"])
                                                  self.fragmentation = not param["fragmentation"]
                                              f.close()

                                      logger.debug("Adding attributes.")

                                      if depth > 0:
                                          depth*=-1
                                      if depth <= -5000:
                                          logger.debug("Setting depth to sea_floor.")
                                          depth = "sea_floor"

                                      self.depth = depth
                                      self.cmap = cmap
                                      self.bins = bins
                                      self.color = color
                                      if linestyle is None:
                                          self.linestyle = ["solid", "dashed", "dotted", "dashdot", (0, (3, 5, 1, 5, 1, 5))]
                                      else:
                                          self.linestyle = linestyle

                                      if lons is None:
                                          logger.warning("Lons is not given - asumming 1 degree resolution over all space.")
                                          lons = np.arange(-180, 180, 1)
                                      if lats is None:
                                          logger.warning("Lats is not given - asumming 1 degree resolution over all space.")
                                          lats = np.arange(-90, 90, 1)
                                      self.lons = lons
                                      self.lats = lats
                                      self.figsize = figsize
                                      self.shrink = shrink

                                      self.title = title
                                      self.suptitle = suptitle
                                      self.diff = diff
                                      self.add = add
                                      self.diffidx = diffidx
                                      self.abs = absolute
                                      self.outfile = outfile
                                      self.fontweight = fontweight

                                      logger.debug("Creating time array.")
                                      self.time = []
                                      for i in range(len(self.objects)):
                                          self.time.append(self.objects[i].get_time_array())

                                      logger.debug("Decrypting status numberings.")
                                      seafloor, stranded, decayed = [], [], []
                                      for i in range(len(self.objects)):
                                          stinfo = get_status_info(self.objects[i].data)
                                          seafloor.append(stinfo[0])
                                          stranded.append(stinfo[1])
                                          decayed.append(stinfo[2])
                                      # self.ice_idx = ice
                                      self.seafloor_idx = seafloor
                                      self.stranded_idx = stranded
                                      self.decayed_idx = decayed

                                      logger.debug("Setting up clipping.")
                                      if clip is None:
                                          self.clip = False
                                      else:
                                          self.clip = True
                                          self.Vmin, self.Vmax = clip

                                      self.areagridpath = areagridpath
                                      self.biomegridpath = biomegridpath
                                      self.mhwIpath = mhwintesitypath
                                      self.group = group
                                      self.loc = locations

                                      self.prop1 = prop1
                                      self.prop2 = prop2
                                      if martincurve is not None:
                                          self.plot_martin_curve = True
                                          self.martin_curve_coef = float(martincurve)
                                      else:self.plot_martin_curve = False
                                      self.cb_units = colorbarlabel
                                      self.xlabel = xlabel
                                      self.ylabel = ylabel
                                      self.xlim = xlim
                                      self.ylim = ylim

                                      self.lw = linewidth

                                      if legend is None:
                                          self.legend = False
                                      else:
                                          self.labels = legend
                                          self.legend = True

                                      self.dpi = dpi

- get_cmap

                                  @staticmethod
                                  def get_cmap(x, map):
                                      cmap = plt.get_cmap(map, len(x))
                                      c = np.arange(1., 1 + len(x))
                                      norm = mpl.colors.BoundaryNorm(np.arange(len(c)+1)+0.5,len(c))
                                      sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
                                      sm.set_array([])
                                      return cmap, sm, c

- shiftedColorMap

                                  @staticmethod
                                  def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shiftedcmap'):
                                      '''
                                      Function to offset the "center" of a colormap. Useful for
                                      data with a negative min and positive max and you want the
                                      middle of the colormap's dynamic range to be at zero.

                                      Input
                                      -----
                                      cmap : The matplotlib colormap to be altered
                                      start : Offset from lowest point in the colormap's range.
                                          Defaults to 0.0 (no lower offset). Should be between
                                          0.0 and `midpoint`.
                                      midpoint : The new center of the colormap. Defaults to
                                          0.5 (no shift). Should be between 0.0 and 1.0. In
                                          general, this should be  1 - vmax / (vmax + abs(vmin))
                                          For example if your data range from -15.0 to +5.0 and
                                          you want the center of the colormap at 0.0, `midpoint`
                                          should be set to  1 - 5/(5 + 15)) or 0.75
                                      stop : Offset from highest point in the colormap's range.
                                          Defaults to 1.0 (no upper offset). Should be between
                                          `midpoint` and 1.0.
                                      '''
                                      cdict = {
                                          'red': [],
                                          'green': [],
                                          'blue': [],
                                          'alpha': []
                                      }

                                      # regular index to compute the colors
                                      reg_index = np.linspace(start, stop, 257)

                                      # shifted index to match the data
                                      shift_index = np.hstack([
                                          np.linspace(0.0, midpoint, 128, endpoint=False),
                                          np.linspace(midpoint, 1.0, 129, endpoint=True)
                                      ])

                                      for ri, si in zip(reg_index, shift_index):
                                          r, g, b, a = cmap(ri)

                                          cdict['red'].append((si, r, r))
                                          cdict['green'].append((si, g, g))
                                          cdict['blue'].append((si, b, b))
                                          cdict['alpha'].append((si, a, a))

                                      newcmap = mpl.colors.LinearSegmentedColormap(name, cdict)
                                      # mpl.cm.register(cmap=newcmap)

                                      return newcmap

- get_colormap_midpoint

                                  @staticmethod
                                  def get_colormap_midpoint(arr):
                                      nans = np.invert(np.logical_or(np.isnan(arr), np.isinf(arr)))
                                      m = np.min(arr[nans])
                                      if m > 0:
                                          m = 0
                                      M = np.max(arr[nans])
                                      return m, 1 - M / (M +np.abs(m)), M

- reset_figure

                                  def reset_figure(self):
                                      fig, ax = plt.subplots(1, 1, figsize = self.figsize)
                                      self.fig = fig
                                      self.ax = ax

- clip_array

                                  def clip_array(self, array):
                                      m, M = np.min(array[np.invert(np.isnan(array))]), np.max(array[np.invert(np.isnan(array))])
                                      row, col = np.where(array < self.Vmin)
                                      array[row, col] = self.Vmin
                                      row, col = np.where(array > self.Vmax)
                                      array[row, col] = self.Vmax
                                      return array

- mass_map

                                  def mass_map(self):
                                      """Plot the mass at depth over a cartopy map."""

                                      plt.close("all")
                                      fig, ax = plt.subplots(1, 1, figsize=self.figsize, subplot_kw={'projection': ccrs.PlateCarree()})

                                      logger.debug("Initialize colormap.")
                                      if self.cmap is None:
                                          if self.diff:
                                              cmap = mpl.cm.RdBu_r
                                          else:
                                              cmap = mpl.cm.Reds
                                      else:
                                          try:
                                              cmap = getattr(mpl.cm, self.cmap)
                                          except AttributeError:
                                              logger.debug("Specified colormap doesn't exist, changing to Reds.")
                                              cmap = "Reds"

                                      h = self.depth

                                      logger.debug("Searching for bad trajectories.")
                                      bad = []
                                      for i in range(len(self.objects)):
                                          bad.append(self.clean_dataset(self.objects[i]))
                                      bad = np.ravel(bad)

                                      logger.debug("Start calculating mass at given depth.")
                                      mass1 = self.zone_crossing_event(self.obj1, self.lons, self.lats, h, bad,self.seafloor_idx[0])

                                      if self.clip and not self.diff and not self.add:
                                          mass1 = self.clip_array(np.copy(mass1))

                                      logger.debug("Finished calculating mass at given depth.")

                                      if self.diff:
                                          if not self.add:
                                              logger.debug("Start calculating mass at given depth for file2.")
                                              mass2 = self.zone_crossing_event(self.obj2, self.lons, self.lats, h, bad,self.seafloor_idx[1])
                                              logger.debug("Finished calculating mass at given depth for file2.")

                                              if self.abs:
                                                  mass = np.copy(mass1) - mass2
                                              else:
                                                  mass = (np.copy(mass1) - mass2) / np.copy(mass1)
                                              if self.clip:
                                                  mass = self.clip_array(np.copy(mass))
                                              m, mid, M = self.get_colormap_midpoint(mass)
                                          else:
                                              mass2 = np.zeros(mass1.shape)
                                              for i in range(1, self.diffidx):
                                                  mass1 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad,self.seafloor_idx[i])
                                              for i in range(self.diffidx, len(self.objects)):
                                                  mass2 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad,self.seafloor_idx[i])
                                              if self.abs:
                                                  mass = np.copy(mass1) - mass2
                                              else:
                                                  mass = (np.copy(mass1) - mass2) / np.copy(mass1)
                                              if self.clip:
                                                  mass = self.clip_array(np.copy(mass))
                                              m, mid, M = self.get_colormap_midpoint(mass)
                                      elif self.add:
                                          logger.debug("Start calculating mass at given depth for other files.")
                                          for i, obj in enumerate(self.objects[1:]):
                                              mass1 += self.zone_crossing_event(obj, self.lons, self.lats, h, bad,self.seafloor_idx[i])
                                          logger.debug("Finished calculating mass at given depth for other files.")
                                          mass = np.copy(mass1)
                                          if self.clip:
                                              mass = self.clip_array(np.copy(mass))
                                      if self.diff or self.add:
                                          NODATA = self.find_grid_cells_with_no_data(self.obj1)
                                          mass[np.isnan(NODATA)] = np.nan
                                      else:
                                          NODATA = self.find_grid_cells_with_no_data(self.obj1)
                                          mass1[np.isnan(NODATA)] = np.nan
                                      logger.debug("Start plotting")
                                      ax.coastlines(zorder = 3, resolution='10m')

                                      extend = 'both' if self.clip else None

                                      if not self.diff and not self.add:
                                          sm = ax.contourf(self.lons, self.lats, mass1.T, 20, transform=ccrs.PlateCarree(), cmap = cmap, zorder = 1,extend=extend, extendfrac='auto')
                                          cb = plt.colorbar(sm, ax = ax, orientation="horizontal", shrink = self.shrink)
                                      else:
                                          if self.diff:
                                              cmap = self.shiftedColorMap(cmap, midpoint = mid, name='shifted')
                                          sm = ax.contourf(self.lons, self.lats, mass.T, 20, transform=ccrs.PlateCarree(), cmap = cmap, zorder = 1,extend = extend, extendfrac='auto')
                                          cb = plt.colorbar(sm, ax = ax, orientation="horizontal", shrink = self.shrink)
                                      if self.cb_units is not None:
                                          cb.set_label(f"{self.cb_units}")
                                      self.axis_setup(ax)

                                      ax.add_feature(cartopy.feature.LAND, zorder=100, edgecolor='k', facecolor = "beige")
                                      gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                linewidth=2, color='k', alpha=0.8, linestyle='--')

                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- find_grid_cells_with_no_data

                                  def find_grid_cells_with_no_data(self, obj:Open):
                                      lon, lat = obj.get_properties_from_list(["lon", "lat"])
                                      lon, lat = lon[0, :], lat[0, :]
                                      X = np.zeros((len(self.lons), len(self.lats)))
                                      for i in range(len(self.lats)):
                                          for j in range(len(self.lons)):
                                              if not np.any((lon == self.lons[j]) & (lat == self.lats[i])):
                                                  X[j, i] = np.nan
                                      return X

- interpolate

                                  @staticmethod
                                  def interpolate(x, y, x0):
                                      if np.any(x == x0):
                                          idx2 = np.where(x == x0)[0][0]
                                          return y[idx2], idx2
                                      idx2 = np.where(x < x0)[0][0]
                                      y1, y2 = y[idx2 - 1], y[idx2]
                                      x1, x2 = x[idx2 - 1], x[idx2]
                                      k = (y2 - y1) / (x2 - x1)
                                      n = y2 - k * x2
                                      return k * x0 + n, idx2

- zone_crossing_event

                                  def zone_crossing_event(self, obj:Open, lons, lats, depth, bad, seafloor_idx):

                                      """Calculate particle properties when crossing a certain depth."""

                                      logger.debug("Reading required simulation properties.")
                                      z, mass, lon, lat, status = obj.get_properties_from_list([
                                                                      "z", "mass", "lon", "lat", "status"
                                                                      ])

                                      return self.mass_sum(lons, lats, depth, z, mass, lon, lat, status, seafloor_idx, bad)

- mass_sum

                                  def mass_sum(self, lons, lats, depth, z, mass, lon, lat, status, sfidx, bad):
                                      '''
                                      Return mass at each grid cell at given depth.
                                      '''
                                      logger.debug("Start summing masses.")
                                      mass_at_depth = np.zeros((len(lons), len(lats)))
                                      for i in range(mass.shape[1]):

                                          #Remove bad particles.
                                          if i in bad:
                                              continue

                                          drifter_trajectory = z[:, i]
                                          drifter_mass = mass[:, i]
                                          m0 = drifter_mass[0]

                                          trajectory_nans = np.invert(np.isnan(drifter_trajectory))

                                          if type(depth) != str:
                                              if np.min(drifter_trajectory[trajectory_nans])>depth:
                                                  continue

                                              m, depth_idx = self.interpolate(drifter_trajectory[trajectory_nans], drifter_mass[trajectory_nans], depth)

                                          else:
                                              #Sea_floor.
                                              # bounce = np.where(drifter_trajectory[:-1] < np.roll(drifter_trajectory, -1)[:-1])[0]
                                              # if len(bounce) > 0:
                                              #     raise ValueError("Depth decreases.")

                                              s = status[:, i]
                                              sea_floor = np.where(s == sfidx)[0]
                                              if len(sea_floor) >0:
                                                  depth_idx = sea_floor[0]
                                              else:
                                                  depth_idx = len(drifter_trajectory)-1
                                              m = drifter_mass[depth_idx]

                                          if np.isnan(m):
                                              continue

                                          drifter_lon = lon[depth_idx, i]
                                          drifter_lat = lat[depth_idx, i]

                                          if np.isnan(drifter_lon) or np.isnan(drifter_lat):
                                              continue

                                          #Find closest position on grid.
                                          lon_grid = np.argmin(np.abs(lons - drifter_lon))
                                          lat_grid = np.argmin(np.abs(lats - drifter_lat))
                                          if not self.diff and not self.abs and not self.add:
                                              m/=m0
                                          mass_at_depth[lon_grid, lat_grid] = mass_at_depth[lon_grid, lat_grid] + m
                                      logger.debug("Finish summing masses.")
                                      return mass_at_depth

- create_timedelta_array

                                  def create_timedelta_array(self, n, objidx):
                                      dt = datetime.strptime(str(self.time[objidx][1]),'%Y-%m-%d %H:%M:%S') - datetime.strptime(str(self.time[objidx][0]),'%Y-%m-%d %H:%M:%S')
                                      dt = dt.total_seconds() / 3600
                                      return np.arange(0, n * dt, dt)

- clean_dataset

                                  def clean_dataset(self, obj:Open):
                                      """Get indicies of trajectories which are on land / have any other problems."""

                                      z, mass, T, lon, lat, status = obj.get_properties_from_list([
                                          "z", "mass", "sea_water_temperature", "lon", "lat", "status"
                                      ])

                                      bad_trajectories = []

                                      for i in range(mass.shape[1]):
                                          drifter_trajectory = z[:, i]
                                          drifter_mass = mass[:, i]

                                          if np.isnan(lat[0, i]):
                                              bad_trajectories.append(i)
                                              continue

                                          trajectory_nans = np.invert(np.isnan(drifter_trajectory))

                                          if (drifter_trajectory[trajectory_nans]>drifter_trajectory[0]).any() or np.any(drifter_mass[1:] >= drifter_mass[0]):
                                              bad_trajectories.append(i)
                                      return bad_trajectories

- current_strength

                                  def current_strength(self):
                                      """Plot horizontal distance traveled between intial position and position at some depth, for ecah grid cell."""
                                      plt.close()
                                      fig, ax = plt.subplots(1, 1, figsize=self.figsize, subplot_kw={'projection': ccrs.PlateCarree()})

                                      logger.debug("Initialize colormap.")
                                      if self.cmap is None:
                                          cmap = "Reds"
                                      else:
                                          try:
                                              cmap = getattr(mpl.cm, self.cmap)
                                          except AttributeError:
                                              logger.debug("Specified colormap doesn't exist, changing to Reds.")
                                              cmap = "Reds"

                                      logger.debug("Searching for bad trajectories.")
                                      bad = self.clean_dataset(self.obj1)

                                      logger.debug("Start calculating distance.")
                                      distance, in_cell, is_at_seafloor, is_stranded = self.calculate_horizontal_distance(self.obj1, self.lons, self.lats, bad)

                                      D = np.copy(distance)
                                      rows, cols = np.where(D == 0) #Land
                                      D[rows, cols] = np.nan
                                      # rows2, cols2 = np.where(in_cell == 0)
                                      # D[rows2, cols2] = np.nan
                                      logger.debug("Start plotting")
                                      ax.add_feature(cartopy.feature.LAND, zorder=2, edgecolor='k', facecolor = "beige")

                                      if self.clip:
                                          D = self.clip_array(np.copy(D))

                                      sm = ax.contourf(self.lons, self.lats, D.T, 20, transform=ccrs.PlateCarree(), cmap = cmap, zorder = 1, extend = "max")
                                      # ax.contour(self.lons, self.lats, in_cell.T, levels = 0, transform = ccrs.PlateCarree(), colors = "red")
                                      ax.contour(self.lons, self.lats, is_at_seafloor.T, levels = 1, transform = ccrs.PlateCarree(), colors = "red", label = "seafloor")
                                      # ax.contour(self.lons, self.lats, is_stranded.T, levels = 1, transform = ccrs.PlateCarree(), colors = "orange", label = "stranded")

                                      # incell_label = Line2D([0], [0], color='red', lw=2)
                                      seafloor_label = Line2D([0], [0], color='red', lw=2)
                                      # stranded_label = Line2D([0], [0], color='orange', lw=2)

                                      # ax.legend([incell_label, seafloor_label, stranded_label], ['moved out of cell', "seafloor", "stranded"], loc='upper right')
                                      ax.legend([seafloor_label], ["seafloor"], loc="upper right")

                                      cb = plt.colorbar(sm, ax = ax, orientation="horizontal", shrink = self.shrink)
                                      if self.cb_units is not None:
                                          cb.set_label(f"{self.cb_units}")
                                      self.axis_setup(ax)

                                      gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                linewidth=2, color='k', alpha=0.8, linestyle='--', zorder = 100)

                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- calculate_horizontal_distance

                                  def calculate_horizontal_distance(self, obj: Open, lons, lats, bad):
                                      '''
                                      Calculate horizontal distanace between initial position and position at given depth for each grid cell.

                                      Returns
                                      -------
                                      distance traveled
                                          ndarray of size (lon,lat) with distances
                                      is_in_cell
                                          boolean ndarray of size (lon, lat) of drifters which are still in the same cell as in the beegining
                                      is_at_seafloor
                                          boolean ndarray of size (lon, lat) of drifters which have reached the seafloor
                                      is_stranded
                                          boolean ndarray of size (lon, lat) of drifters which were stranded on the castline'''

                                      from haversine import haversine, Unit
                                      from tqdm import trange

                                      logger.debug("Reading required simulation properties.")

                                      lon, lat, z, status = obj.get_properties_from_list([
                                          "lon", "lat", "z", "status"
                                      ])

                                      distance_traveled = np.zeros((len(lons), len(lats)))
                                      is_in_cell = np.ones((len(lons), len(lats)))
                                      is_stranded = np.zeros((len(lons), len(lats)))
                                      is_at_seafloor = np.zeros((len(lons), len(lats)))
                                      for i in trange(lon.shape[1]):
                                          stranded = False
                                          reached_sea_floor = False
                                          #Remove bad particles.
                                          if i in bad:
                                              # progress.update(1)
                                              continue

                                          drifter_lon = lon[:, i]
                                          drifter_lat = lat[:, i]
                                          drifter_z = z[:, i]
                                          drifter_status = status[:, i]
                                          nans = np.invert(np.isnan(drifter_lat))
                                          if np.min(drifter_z[nans]) > self.depth:
                                              if len(np.where(drifter_status == self.stranded_idx[0])[0]) > 0:
                                                  idx = np.where(drifter_status == self.stranded_idx[0])[0][0]
                                                  stranded = True
                                              elif len(np.where(drifter_status == self.decayed_idx[0])[0]) > 0:
                                                  idx = np.where(drifter_status == self.decayed_idx[0])[0][0]
                                                  stranded = True
                                              else:
                                                  idx = np.where(drifter_status == self.seafloor_idx[0])[0]
                                                  if len(idx) == 0:
                                                      logger.warning("Drifter at initial position ({}, {}) has not been deactivated/passed the specified depth.".format(drifter_lat[0], drifter_lon[0]))
                                                      logger.warning("Taking the last position.")
                                                      idx = len(drifter_z[nans]) - 1
                                                  else:
                                                      idx = idx[0]
                                                      reached_sea_floor = True
                                              phi = drifter_lat[idx]
                                              lam = drifter_lon[idx]
                                          else:
                                              phi, idx = self.interpolate(drifter_z, drifter_lat, self.depth)
                                              lam, _ = self.interpolate(drifter_z, drifter_lon, self.depth)

                                          pos1 = (drifter_lat[0], drifter_lon[0])
                                          pos2 = (phi, lam)

                                          lon_grid = np.argmin(np.abs(lons - pos1[1]))
                                          lat_grid = np.argmin(np.abs(lats - pos1[0]))
                                          if np.abs(pos2[0] - pos1[0]) > (lats[1] - lats[0]) / 2 or np.abs(pos1[1] - pos2[1]) > (lons[1] - lons[0]) / 2:
                                              is_in_cell[lon_grid, lat_grid] = 0
                                              d = haversine(pos1, pos2, unit = Unit.KILOMETERS)
                                              distance_traveled[lon_grid, lat_grid] = d
                                          else:
                                              d = haversine(pos1, pos2, unit = Unit.KILOMETERS)
                                              distance_traveled[lon_grid, lat_grid] = d
                                          if stranded:
                                              is_stranded[lon_grid, lat_grid] = 1
                                          if reached_sea_floor:
                                              is_at_seafloor[lon_grid, lat_grid] = 1
                                      return distance_traveled, is_in_cell, is_at_seafloor, is_stranded

- drifter_locations

                                  def drifter_locations(self):
                                      '''Plot locations of specified drifters on a world map.'''
                                      from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
                                      import matplotlib.style as style

                                      style.use('seaborn-v0_8-bright')

                                      if self.loc is None:
                                          raise AttributeError("Locations are not specified.")

                                      plt.close()
                                      fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.PlateCarree()}, figsize = self.figsize)

                                      logger.debug("Searching for bad trajectories.")
                                      bad = self.clean_dataset(self.obj1)
                                      logger.debug("Reading required simulation properties.")

                                      lon, lat = self.obj1.get_properties_from_list(["lon", "lat"])

                                      k = 0
                                      markers = "s^vox"

                                      for i, j in self.loc:
                                          lon_idx = np.where(lon[0, :] == i)
                                          lat_idx = np.where(lat[0, :] == j)
                                          idx = np.intersect1d(lon_idx, lat_idx)
                                          if len(idx) == 0:
                                              logger.info("Could not find drifter with starting position ({}, {}).".format(i, j))
                                              continue
                                          else:
                                              idx = idx[0]

                                          if idx in bad:
                                              logger.info("({}, {}) is a bad trajectory and will not be included.".format(i, j))
                                              continue
                                          if self.color is not None:
                                              if len(self.color) >1:
                                                  if k >= len(self.color):
                                                      color = None
                                                  else:
                                                      color = self.color[k]
                                              else:
                                                  color = self.color[0]
                                          else:
                                              color = None

                                          ax.plot(lon[0, idx], lat[0, idx], markers[k%len(markers)], color = color, markersize = 10, label = self.labels[k] if self.legend else None, zorder = 100)
                                          k += 1

                                      ax.set_global()
                                      ax.add_feature(cartopy.feature.LAND, facecolor="beige",edgecolor='black', zorder = 1)
                                      ax.coastlines(zorder = 2)
                                      gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                linewidth=2, color='k', alpha=0.8, linestyle='--')
                                      self.axis_setup(ax)
                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- martin_curve

                                  def martin_curve(self, z, flux_100 = None):
                                      if flux_100 is None:
                                          return (np.abs(z) / 100) ** (-self.martin_curve_coef)
                                      else: return flux_100 * (np.abs(z) / 100) ** (-self.martin_curve_coef)

- drifter_properties

                                  def drifter_properties(self):
                                      '''Plot drifter prop1 vs prop2 at a specified location.'''
                                      import matplotlib.style as style

                                      style.use('tableau-colorblind10')

                                      plt.close()

                                      if self.loc is None:
                                          raise AttributeError("Locations are not specified.")

                                      self.fig, self.ax = plt.subplots(1, len(self.loc), figsize = self.figsize)
                                      try:
                                          self.ax[0]
                                      except TypeError:
                                          self.ax = [self.ax]
                                      logger.debug("Searching for bad trajectories.")
                                      bad = []
                                      for i in range(len(self.objects)):
                                          bad.append(self.clean_dataset(self.objects[i]))
                                      logger.debug("Reading required simulation properties.")

                                      if self.prop1 is None:
                                          logger.info("Property 1 not specified. Setting x to time.")
                                          self.prop1 = "time"
                                      if self.prop2 is None:
                                          logger.info("Property 2 not specified. Setting y to mass.")
                                          self.prop2 = "mass"

                                      lines = []

                                      m = 0
                                      for i, j in self.loc:
                                          for k, obj in enumerate(self.objects):
                                              lon, lat = obj.get_properties_from_list(["lon", "lat"])
                                              lon_idx = np.where(lon[0, :] == i)
                                              lat_idx = np.where(lat[0, :] == j)
                                              idx = np.intersect1d(lon_idx, lat_idx)

                                              if len(idx) == 0:
                                                  logger.info("Could not find drifter with starting position ({}, {}).".format(i, j))
                                                  continue
                                              else:
                                                  idx = idx[0]

                                              if idx in bad[k]:
                                                  logger.info("({}, {}) is a bad trajectory and will not be included.".format(i, j))
                                                  continue

                                              if self.prop1 == "total_horizontal_velocity":
                                                  u = obj.get_property("x_sea_water_velocity")
                                                  v = obj.get_property("y_sea_water_velocity")
                                                  x = np.sqrt(u ** 2 + v ** 2)
                                                  x = np.ma.filled(x, np.nan)
                                              elif self.prop1 == "tau":
                                                  T = obj.get_property("sea_water_temperature")
                                                  T = np.ma.filled(T, np.nan)
                                                  T[T< 0]=0
                                                  decay = self.decay_type[k]
                                                  if decay == "linear":
                                                      decayrate = 0.064 * T + 0.02
                                                  else:
                                                      decayrate = 0.140 * np.exp(0.145 * T)
                                                  x = 1 / decayrate
                                                  x *= 24
                                              elif self.prop1 != "time":
                                                  x = obj.get_property(self.prop1)
                                                  x = np.ma.filled(x, np.nan)
                                              else:
                                                  x = self.create_timedelta_array(lon.shape[0], k)

                                              if self.prop2 == "total_horizontal_velocity":
                                                  u = obj.get_property("x_sea_water_velocity")
                                                  v = obj.get_property("y_sea_water_velocity")
                                                  y = np.sqrt(u ** 2 + v ** 2)
                                                  Y = np.ma.filled(x, np.nan)
                                              elif self.prop2 == "tau":
                                                  T = obj.get_property("sea_water_temperature")
                                                  T = np.ma.filled(T, np.nan)
                                                  T[T< 0]=0
                                                  decay = self.decay_type[k]
                                                  if decay == "linear":
                                                      decayrate = 0.064 * T + 0.02
                                                  else:
                                                      decayrate = 0.140 * np.exp(0.145 * T)
                                                  y = 1 / decayrate
                                                  Y *= 24
                                              elif self.prop2 != "time":
                                                  y = obj.get_property(self.prop2)
                                                  y = np.ma.filled(y, np.nan)
                                              else:
                                                  y = self.create_timedelta_array(lon.shape[0], k)

                                              if self.prop1 == "time":
                                                  X = x
                                              else:
                                                  X = x[:, idx]
                                                  if not self.abs and self.prop1 == "mass":
                                                      X /= X[0]

                                              if self.prop2 == "time":
                                                  Y = y
                                              else:
                                                  Y = y[:, idx]
                                                  if not self.abs and self.prop2 == "mass":
                                                      Y /= Y[0]

                                              if self.color is not None:
                                                  try:
                                                      color = self.color[k % len(self.color)]
                                                  except IndexError:
                                                      logger.warning(f"Setting color to None for object number {k}", exc_info=IndexError)
                                                      color = None
                                              else:
                                                  color = None
                                              linestyle = self.linestyle[k % len(self.linestyle)]
                                              line, = self.ax[m].plot(X, Y, lw = self.lw, linestyle = linestyle, color = color)
                                              lines.append(line)

                                          if self.plot_martin_curve:
                                              if self.abs:
                                                  #ADD FLUXES
                                                  pass
                                              else:
                                                  line, = self.ax[m].plot(self.martin_curve(Y), Y, lw = self.lw, linestyle = "dotted")
                                              if m ==0 and self.legend:
                                                  lines.append(line)
                                                  self.labels.append("Martin Curve")

                                          self.axis_setup(self.ax[m], False, True)
                                          if self.suptitle is not None:
                                              self.ax[m].set_title(self.suptitle[m], fontweight = self.fontweight)

                                          # self.ax[m].tick_params(labelbottom=False,labeltop=True)
                                          self.ax[m].xaxis.set_ticks_position('top')
                                          self.ax[m].xaxis.set_label_position('top')
                                          m+=1
                                      if self.title is not None:
                                          self.fig.suptitle(self.title, fontweight = self.fontweight)
                                      if self.legend: self.ax[0].legend(lines, [f"{label}" for label in self.labels])
                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- get_mass_sum_at_depth

                                  def get_mass_sum_at_depth(self):

                                      '''Return sum of mass over all grid points at given depth.'''

                                      M = self.zone_crossing_event(self.obj1, self.lons, self.lats, self.depth, self.clean_dataset(self.obj1))
                                      if not self.abs:
                                          x = len(np.where(M > 0)[0])
                                      else:
                                          x = 1
                                      return np.sum(M[np.invert(np.isnan(M))]) / x

- get_biome_weighted_mass_at_depth

                                  def get_biome_weighted_mass_at_depth(self, obj = None):
                                      '''Return sum of mass over all grid points in the 4 biomes at given depth.'''
                                      logger.debug("Reading required simulation properties.")
                                      if obj is None:
                                          obj = self.obj1
                                      z, mass, lon, lat, biome, status = obj.get_properties_from_list([
                                          "z", "mass", "lon", "lat", "origin_marker", "status"
                                      ])

                                      logger.debug("Start summing masses.")

                                      depth = self.depth
                                      mass_by_biomes = np.zeros(4)
                                      for i in range(mass.shape[1]):

                                          drifter_trajectory = z[:, i]
                                          drifter_mass = mass[:, i]
                                          b = int(biome[0, i])

                                          trajectory_nans = np.invert(np.isnan(drifter_trajectory))

                                          if type(depth) != str:
                                              if np.min(drifter_trajectory[trajectory_nans])>depth:
                                                  continue

                                              m, depth_idx = self.interpolate(drifter_trajectory[trajectory_nans], drifter_mass[trajectory_nans], depth)

                                          else:

                                              s = status[:, i]
                                              sea_floor = np.where(s == self.seafloor_idx[0])[0]
                                              if len(sea_floor) >0:
                                                  depth_idx = sea_floor[0]
                                              else:
                                                  depth_idx = len(drifter_trajectory)-1
                                              m = drifter_mass[depth_idx]

                                          if np.isnan(m):
                                              continue

                                          mass_by_biomes[b] += m

                                      logger.debug("Finish summing masses.")
                                      return mass_by_biomes

- get_mean_biome_mass_flux_at_depth

                                  def get_mean_biome_mass_flux_at_depth(self, obj = None):
                                      '''Return sum of mass flux (per area) over all grid points in the 4 biomes at a given depth.'''
                                      logger.debug("Reading required simulation properties.")
                                      if obj is None:
                                          obj = self.obj1
                                      z, mass, lon, lat, biome, status = obj.get_properties_from_list([
                                          "z", "mass", "lon", "lat", "origin_marker", "status"
                                      ])

                                      area = self.load_area()

                                      logger.debug("Start summing masses.")

                                      depth = self.depth
                                      mass_by_biomes = [[], [], [], []]
                                      for i in range(mass.shape[1]):

                                          drifter_trajectory = z[:, i]
                                          drifter_mass = mass[:, i]
                                          if drifter_mass[0] <=0:continue
                                          b = int(biome[0, i])

                                          trajectory_nans = np.invert(np.isnan(drifter_trajectory))

                                          if type(depth) != str:
                                              if np.min(drifter_trajectory[trajectory_nans])>depth:
                                                  continue

                                              m, depth_idx = self.interpolate(drifter_trajectory[trajectory_nans], drifter_mass[trajectory_nans], depth)

                                          else:
                                              #Sea_floor.
                                              # bounce = np.where(drifter_trajectory[:-1] < np.roll(drifter_trajectory, -1)[:-1])[0]
                                              # if len(bounce) > 0:
                                              #     raise ValueError("Depth decreases.")

                                              s = status[:, i]
                                              sea_floor = np.where(s == self.seafloor_idx[0])[0]
                                              if len(sea_floor) >0:
                                                  depth_idx = sea_floor[0]
                                              else:
                                                  depth_idx = len(drifter_trajectory)-1
                                              m = drifter_mass[depth_idx]

                                          if np.isnan(m):
                                              continue
                                          A = area[np.where(self.lats == lat[0, i])[0][0], np.where(self.lons == lon[0, i])[0][0]]
                                          if A == 0:
                                              continue
                                          mass_by_biomes[b].append(m / A)

                                      logger.debug("Finish summing masses.")
                                      return np.asarray([np.mean(i) for i in mass_by_biomes])

- load_area

                                  def load_area(self):
                                      return np.load(self.areagridpath)

- mean_export_biome_flux

                                  def mean_export_biome_flux(self):
                                      '''
                                      Bar plot of mean fluxes of multiple groups (e.g. Cnidaria egestion, Ctenophora mortality...) for each biome.
                                      '''
                                      logger.debug("Setting up figure.")
                                      plt.close()
                                      fig = plt.figure(figsize = self.figsize)
                                      outer_grid = gridspec.GridSpec(1, 4, wspace=0.05)
                                      ax = []
                                      biome_titles = ["HCSS", "LC", "HCPS", "COAST"]
                                      for i in range(4):
                                          ax_outer = fig.add_subplot(outer_grid[i])
                                          ax_outer.axis('off')
                                          ax_outer.set_title(biome_titles[i], y = 1.05)
                                          innergrid = gridspec.GridSpecFromSubplotSpec(1, 2, subplot_spec=outer_grid[i], wspace=0, hspace=0)
                                          for j in range(2):
                                              if not (i == 0 and j == 0):
                                                  ax.append(plt.Subplot(fig, innergrid[j], sharey=ax[0]))
                                              else:
                                                  ax.append(plt.Subplot(fig, innergrid[j]))
                                      logger.debug("Figure set.")
                                      bottom = np.zeros((len(ax), len(self.objects) // 8+1))
                                      if self.color is None:
                                          colors = ["blue", "orange", "green"]
                                      else:
                                          colors = self.color
                                      if self.xlabel is not None:
                                          xticks = self.xlabel.split(",")
                                      else:
                                          xticks = np.arange(len(self.objects) // 8)
                                      j = 0
                                      for i, obj in enumerate(self.objects):
                                          flux = self.get_mean_biome_mass_flux_at_depth(obj)

                                          for k in range(4):
                                              y = np.zeros(len(xticks))
                                              y[j] = flux[k]
                                              ax[i%2 + 2*k].bar(xticks, y, bottom = bottom[i%2 + 2*k], color = colors[(i%6) // 2], alpha = 0.7)
                                              bottom[i%2 + 2*k]+=y
                                          if i % 6 == 5:
                                              j+=1
                                      if self.suptitle is None:
                                          suptitle = ["Carcasses", "Fecal Matter"]
                                      else:
                                          suptitle = self.suptitle
                                      for i in range(len(ax)):
                                          if i > 0:
                                              ax[i].yaxis.set_tick_params(labelleft=False)
                                          if i % 2 == 0:
                                              ax[i].set_title(suptitle[0], fontsize=15, fontweight=self.fontweight)
                                          else:
                                              ax[i].set_title(suptitle[1], fontsize=15, fontweight=self.fontweight)
                                          ax[i].tick_params(axis='x', labelrotation=90, labelsize = 15)

                                          fig.add_subplot(ax[i])

                                      if self.ylabel is not None:
                                          ax[0].set_ylabel(f"{self.ylabel}")

                                      if self.legend:
                                          lines = [mpl.patches.Patch(facecolor=colors[0]),
                                                   mpl.patches.Patch(facecolor=colors[1]),
                                                   mpl.patches.Patch(facecolor=colors[2])]
                                          fig.legend(lines, [f"{label}" for label in self.labels], loc = "center right")
                                      if self.xlim is not None: ax[0].set_xlim(*self.xlim)
                                      if self.ylim is not None: ax[0].set_ylim(*self.ylim)
                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- mass_flux_map

                                  def mass_flux_map(self):
                                      """Plot the mass flux reached at given depth on world map.
                                      """

                                      plt.close("all")
                                      fig, ax = plt.subplots(1, 1, figsize=self.figsize, subplot_kw={'projection': ccrs.PlateCarree()})

                                      logger.debug("Initialize colormap.")
                                      if self.cmap is None:
                                          if self.diff:
                                              cmap = mpl.cm.RdBu_r
                                          else:
                                              cmap = mpl.cm.Reds
                                      else:
                                          try:
                                              cmap = getattr(mpl.cm, self.cmap)
                                          except AttributeError:
                                              logger.debug("Specified colormap doesn't exist, changing to Reds.")
                                              cmap = "Reds"

                                      h = self.depth
                                      area = np.load(self.areagridpath).T
                                      i = np.where(np.arange(-180, 180, 1) == np.min(self.lons))[0][0]
                                      j = np.where(np.arange(-180, 180, 1) == np.max(self.lons))[0][0]
                                      k = np.where(np.arange(-90, 90, 1) == np.min(self.lats))[0][0]
                                      l = np.where(np.arange(-90, 90, 1) == np.max(self.lats))[0][0]
                                      area = area[i:j+1, k:l+1]
                                      logger.debug("Searching for bad trajectories.")
                                      bad = []
                                      for i in range(len(self.objects)):
                                          bad.append(self.clean_dataset(self.objects[i]))
                                          break

                                      logger.debug("Start calculating mass at given depth.")
                                      flux1 = self.zone_crossing_event(self.obj1, self.lons, self.lats, h, bad, self.seafloor_idx[0]) / area

                                      if self.clip and not self.diff and not self.add:
                                          flux1 = self.clip_array(np.copy(flux1))

                                      logger.debug("Finished calculating mass at given depth.")

                                      if self.diff:
                                          if not self.add:
                                              logger.debug("Start calculating mass at given depth for file2.")
                                              flux2 = self.zone_crossing_event(self.obj2, self.lons, self.lats, h, bad,self.seafloor_idx[1]) / area
                                              logger.debug("Finished calculating mass at given depth for file2.")

                                              if self.abs:
                                                  flux = np.copy(flux1) - flux2
                                              else:
                                                  flux= (np.copy(flux1) - flux2) / np.copy(flux1)
                                              if self.clip:
                                                  flux = self.clip_array(np.copy(flux))
                                              m, mid, M = self.get_colormap_midpoint(flux)
                                          else:
                                              flux2 = np.zeros(flux1.shape)
                                              for i in range(1, self.diffidx):
                                                  flux1 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad, self.seafloor_idx[i]) / area
                                              for i in range(self.diffidx, len(self.objects)):
                                                  flux2 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad, self.seafloor_idx[i]) / area
                                              if self.abs:
                                                  flux = np.copy(flux1) - flux2
                                              else:
                                                  flux= (np.copy(flux1) - flux2) / np.copy(flux1)
                                              if self.clip:
                                                  flux = self.clip_array(np.copy(flux))
                                              m, mid, M = self.get_colormap_midpoint(flux)
                                      elif self.add:
                                          logger.debug("Start calculating mass at given depth for other files.")
                                          for i, obj in enumerate(self.objects[1:]):
                                              flux1 += self.zone_crossing_event(obj, self.lons, self.lats, h, bad, self.seafloor_idx[i]) / area
                                          logger.debug("Finished calculating mass at given depth for other files.")
                                          flux = np.copy(flux1)
                                          if self.clip:
                                              flux = self.clip_array(np.copy(flux))
                                      if self.diff or self.add:
                                          NODATA = self.find_grid_cells_with_no_data(self.obj1)
                                          flux[np.isnan(NODATA)] = np.nan
                                      else:
                                          NODATA = self.find_grid_cells_with_no_data(self.obj1)
                                          flux1[np.isnan(NODATA)] = np.nan
                                      logger.debug("Start plotting")
                                      ax.coastlines(zorder = 3, resolution='10m')

                                      extend = 'both' if self.clip else None

                                      if not self.diff and not self.add:
                                          sm = ax.contourf(self.lons, self.lats, flux1.T, 20, transform=ccrs.PlateCarree(), cmap = cmap, zorder = 1,extend=extend, extendfrac='auto')
                                          cb = plt.colorbar(sm, ax = ax, orientation="horizontal", shrink = self.shrink)
                                      else:
                                          if self.diff:
                                              cmap = self.shiftedColorMap(cmap, midpoint = mid, name='shifted')
                                          sm = ax.contourf(self.lons, self.lats, flux.T, 20, transform=ccrs.PlateCarree(), cmap = cmap, zorder = 1,extend = extend, extendfrac='auto')
                                          cb = plt.colorbar(sm, ax = ax, orientation="horizontal", shrink = self.shrink)
                                      if self.cb_units is not None:
                                          cb.set_label(f"{self.cb_units}")
                                      self.axis_setup(ax)

                                      ax.add_feature(cartopy.feature.LAND, zorder=2, edgecolor='k', facecolor = "beige")
                                      gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                linewidth=2, color='k', alpha=0.8, linestyle='--', zorder = 100)

                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- animate_3D

                                  def animate_3D(self):
                                      """A 3d animation of the simulation."""

                                      plt.close()

                                      fig = plt.figure(figsize=self.figsize)
                                      ax = fig.add_subplot(projection='3d')

                                      logger.debug("Reading required simulation properties.")
                                      lon, lat, z, mass = self.obj1.get_properties_from_list([
                                          "lon", "lat", "z", "mass"
                                      ])

                                      logger.debug("Initialize colormap.")
                                      if self.cmap is None:
                                          logger.warning("Colormap hasn't been provided, using jet.")
                                          cmap = mpl.cm.jet
                                      else:
                                          try:
                                              cmap = getattr(mpl.cm, self.cmap)
                                          except AttributeError:
                                              logger.debug("Specified colormap doesn't exist, changing to jet.")
                                              cmap = "jet"

                                      if self.fragmentation:
                                          k = 100
                                      else:
                                          k = 1 / mass[0, :]
                                      dt = (self.time[0][1]-self.time[0][0]).total_seconds()
                                      idx = np.where(np.invert(np.logical_or(np.isnan(mass[1, :]), np.isnan(z[1, :]))))[0][0]
                                      try:
                                          w0 = abs(self.initial_velocities[0]) * 24 * 3600
                                      except Exception:
                                          w0 = abs((z[1, idx] - z[0, idx]) * (mass[0, idx]/mass[1, idx]) ** (1/6) / dt)*24*3600
                                      sc = ax.scatter(lon[0, :], lat[0, :], z[0, :], s=mass[0, :] * k, c = w0 * np.ones(len(mass[0, :])), cmap = cmap, vmin=0, vmax = w0)

                                      logger.debug("Setting up axes properties.")
                                      ax.yaxis.labelpad=30
                                      ax.xaxis.labelpad=30
                                      ax.zaxis.labelpad=30
                                      ax.set_xlabel("lon")
                                      ax.set_ylabel("lat")
                                      ax.set_zlabel("z")
                                      xlim = [np.min(lon[np.invert(np.isnan(lon))]), np.max(lon[np.invert(np.isnan(lon))])]
                                      ylim = [np.min(lat[np.invert(np.isnan(lat))]), np.max(lat[np.invert(np.isnan(lat))])]
                                      zlim = [np.min(z[np.invert(np.isnan(z))]), np.max(z[np.invert(np.isnan(z))])]
                                      ax.set_xlim(*xlim)
                                      ax.set_ylim(*ylim)
                                      ax.set_zlim(*zlim)
                                      ax.set_title(self.time[0][0])
                                      ax.view_init(elev=10, azim=45, roll=0)

                                      cb = plt.colorbar(sc, shrink = self.shrink)
                                      cb.set_label(r"$|w(t)|\,\mathrm{[m\,day^{-1}]}$")

                                      def update(i):
                                          ax.cla()
                                          ax.set_xlim(*xlim)
                                          ax.set_ylim(*ylim)
                                          ax.set_zlim(*zlim)
                                          ax.yaxis.labelpad=30
                                          ax.xaxis.labelpad=30
                                          ax.zaxis.labelpad=30
                                          ax.set_xlabel("lon")
                                          ax.set_ylabel("lat")
                                          ax.set_zlabel("z")
                                          idx = np.invert(np.isnan(lon[i, :]))
                                          w = abs((z[i - 1, idx] - z[i, idx]) / dt)*24*3600
                                          if self.fragmentation:
                                              q = k
                                          else:
                                              q = k[idx]
                                          sc = ax.scatter(lon[i, idx], lat[i, idx], z[i, idx], s=mass[i, idx] * q, c = w, cmap = cmap, vmin=0, vmax = w0)
                                          cb.update_normal(sc)
                                          ax.set_title(self.time[0][i])
                                          ax.view_init(elev=10, azim=45, roll=0)

                                      frames = len(self.time[0])

                                      logger.debug("Initializing animation.")
                                      anim=animation.FuncAnimation(fig, update, blit=False, frames = frames, interval=100)

                                      if self.outfile is not None:
                                          logger.debug("Adding PillowWriter to animation.")
                                          writer = animation.PillowWriter(fps=10)
                                          logger.debug("Saving animation.")
                                          anim.save(self.outfile, writer = writer)
                                      else:
                                          plt.show()

- mean_lon_mass_flux

                                  def mean_lon_mass_flux(self):
                                      """Plot latitude vs longitude averaged mass fluxes"""
                                      area = np.load(self.areagridpath).T
                                      logger.debug("Searching for bad trajectories.")
                                      bad = self.clean_dataset(self.objects[0])
                                      bad = np.ravel(bad)

                                      logger.debug("Start calculating mass at given depth.")
                                      if self.legend is not None:
                                          model_num = len(self.labels)
                                      else:
                                          model_num = 1
                                      sim_num = len(self.objects) // model_num
                                      for mnum in range(model_num):
                                          flux1 = self.zone_crossing_event(self.objects[mnum*sim_num], self.lons, self.lats, self.depth, bad,self.seafloor_idx[mnum*sim_num]) / area
                                          if self.add:
                                              logger.debug("Start calculating mass at given depth for other files.")
                                              for i, obj in enumerate(self.objects[mnum*sim_num + 1:mnum*sim_num+sim_num]):
                                                  flux1 += self.zone_crossing_event(obj, self.lons, self.lats, self.depth, bad, self.seafloor_idx[mnum*sim_num + 1 + i]) / area
                                              logger.debug("Finished calculating mass at given depth for other files.")
                                          flux = np.copy(flux1)
                                          flux[np.isnan(flux)] = 0
                                          flux_mask = flux == 0
                                          flux = np.ma.array(data = flux, mask = flux_mask)
                                          mean_lon_flux = np.ma.mean(flux, axis = 0)
                                          # print(mean_lon_flux.shape)
                                          if self.color is not None:
                                              try:
                                                  color = self.color[mnum]
                                              except IndexError:
                                                  color = "black"
                                                  logger.warning(f"Setting color to None.", exc_info=IndexError)
                                          else:
                                              color = "black"
                                          if self.linestyle is not None:
                                              try:
                                                  ls = self.linestyle[mnum]
                                              except IndexError:
                                                  ls = "solid"
                                                  logger.warning(f"Setting linestyle to solid.", exc_info=IndexError)
                                          else:
                                              ls = "solid"
                                          if self.legend is not None:
                                              label = self.labels[mnum]
                                          else:
                                              label = None
                                          plt.plot(self.lats, mean_lon_flux, color = color, linewidth = self.lw, label = label, ls = ls)
                                      if self.xlabel is not None: plt.xlabel(f"{self.xlabel}")
                                      if self.ylabel is not None: plt.ylabel(f"{self.ylabel}")
                                      if self.xlim is not None:
                                          plt.xlim(*self.xlim)
                                      else:
                                          plt.xlim((-90, 90))
                                          plt.xticks([-90, -45, 0, 45, 90], [r"$-90^\circ$", r"$-45^\circ$",r"$0^\circ$", r"$45^\circ$", r"$90^\circ$"])
                                      if self.ylim is not None: plt.ylim(*self.ylim)
                                      if self.legend is not None: plt.legend()
                                      plt.grid()
                                      if self.title is not None:
                                          plt.title(self.title)
                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- mass_flux_distribution

                                  def mass_flux_distribution(self):
                                      """Mass flux distributions at some depth for different grouping of tarcers."""
                                      area = np.load(self.areagridpath).T
                                      logger.debug("Searching for bad trajectories.")
                                      bad = []
                                      for i in range(len(self.objects)):
                                          bad.append(self.clean_dataset(self.objects[i]))
                                      bad = np.ravel(bad)

                                      logger.debug("Start calculating mass at given depth.")
                                      flux1 = self.zone_crossing_event(self.obj1, self.lons, self.lats, self.depth, bad,self.seafloor_idx[0]) / area
                                      if self.add:
                                          logger.debug("Start calculating mass at given depth for other files.")
                                          for i, obj in enumerate(self.objects[1:]):
                                              flux1 += self.zone_crossing_event(obj, self.lons, self.lats, self.depth, bad,self.seafloor_idx[i]) / area
                                          logger.debug("Finished calculating mass at given depth for other files.")
                                      flux = np.copy(flux1)
                                      flux[np.isnan(flux)] = 0
                                      flux_mask = flux == 0
                                      if self.group is None or self.group == "none":
                                          if self.color is not None:
                                              color = self.color[0]
                                          else:
                                              color = None
                                          self.ax.hist(flux[np.invert(flux_mask)], bins=self.bins, color = color, density = True)
                                      elif self.group == "biome":
                                          biome_names = ["HCSS", "LC", "HCPS", "COAST"]
                                          biomes = np.load(self.biomegridpath).T
                                          for i in range(4):
                                              if self.color is not None:
                                                  color = self.color[i % len(self.color)]
                                              else:
                                                  color = None
                                              mask = np.invert(flux_mask) & (biomes == i)
                                              self.ax.hist(flux[mask], bins = self.bins, color = color, alpha = 0.8, label = biome_names[i], density=True)
                                      elif self.group in ["lonmean", "latmean"]:
                                          if self.color is not None:
                                              color = self.color[0]
                                          else:
                                              color = None
                                          flux = np.ma.array(data = flux, mask = flux_mask)
                                          self.ax.hist(np.ma.mean(flux, axis = 0 if self.group == "lonmean" else 1), bins = self.bins, color = color, density=True)
                                      self.ax.axvline(np.ma.mean(flux[mask]), color = "black", lw = self.lw, linestyle = "dashed", label = f"povprečje; {round(np.ma.mean(flux[mask]), 2)}")
                                      self.ax.legend()
                                      self.axis_setup(self.ax)

                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- mhw_correlations

                                      def mhw_correlations(self):
                                      logger.debug("Searching for bad trajectories.")
                                      bad = []
                                      for i in range(len(self.objects)):
                                          bad.append(self.clean_dataset(self.objects[i]))
                                      logger.debug("Reading required simulation properties.")
                                      h = self.depth
                                      area = np.load(self.areagridpath).T
                                      I = np.load(self.mhwIpath).T
                                      flux1 = np.zeros((360, 180))
                                      flux2 = np.zeros((360, 180))
                                      for i in range(self.diffidx):
                                          flux1 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad,self.seafloor_idx[i]) / area
                                      for i in range(self.diffidx, len(self.objects)):
                                          flux2 += self.zone_crossing_event(self.objects[i], self.lons, self.lats, h, bad,self.seafloor_idx[i]) / area
                                      flux = (np.copy(flux1) - flux2) / np.copy(flux1)
                                      self.ax.scatter(I.ravel(), flux.ravel(), color = "k")
                                      mask = np.invert((flux1 == 0) | (flux2==0) | np.isnan(flux1) | np.isnan(flux2) | np.isnan(I))
                                      if self.color is not None:
                                          color = self.color[0]
                                      else:
                                          color = "k"
                                      self.ax.scatter(I.ravel()[mask.ravel()], flux.ravel()[mask.ravel()], color = color)
                                      cov = np.cov(I.ravel()[mask.ravel()], flux.ravel()[mask.ravel()])
                                      pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
                                      ell_radius_x = np.sqrt(1 + pearson)
                                      ell_radius_y = np.sqrt(1 - pearson)
                                      ellipse = Ellipse((0, 0), width=ell_radius_x * 2, height=ell_radius_y * 2, edgecolor = "red", facecolor = 'none', linewidth = 3)
                                      scale_x = np.sqrt(cov[0, 0]) * 5
                                      mean_x = np.mean(I.ravel()[mask.ravel()])

                                      # calculating the standard deviation of y ...
                                      scale_y = np.sqrt(cov[1, 1]) * 5
                                      mean_y = np.mean(flux.ravel()[mask.ravel()])

                                      transf = transforms.Affine2D()
                                      transf.rotate_deg(45)
                                      transf.scale(scale_x, scale_y)
                                      transf.translate(mean_x, mean_y)

                                      ellipse.set_transform(transf + self.ax.transData)
                                      self.ax.add_patch(ellipse)
                                      self.axis_setup(self.ax, None, True)
                                      plt.tight_layout()

                                      if self.outfile is not None:
                                          logger.debug("Saving output file.")
                                          plt.savefig(self.outfile, dpi = self.dpi, bbox_inches = "tight")
                                      else:
                                          plt.show()

- get_mass_over_mhw_area

                                  def get_mass_over_mhw_area(self, obj = None):
                                      '''Return sum of mass over all grid points in the 4 biomes at given depth.'''
                                      logger.debug("Reading required simulation properties.")
                                      if obj is None:
                                          obj = self.obj1
                                      loc = self.loc
                                      if loc is None:
                                          raise AttributeError("Missing mhw Area. Locations must be added.")
                                      lonmin = min([i[0] for i in loc])
                                      lonmax = max([i[0] for i in loc])
                                      latmin = min([i[1] for i in loc])
                                      latmax = max([i[1] for i in loc])
                                      z, mass, lon, lat, biome, status = obj.get_properties_from_list([
                                          "z", "mass", "lon", "lat", "origin_marker", "status"
                                      ])

                                      logger.debug("Start summing masses.")

                                      depth = self.depth
                                      mass_by_biomes = 0
                                      for i in range(mass.shape[1]):

                                          drifter_trajectory = z[:, i]
                                          drifter_mass = mass[:, i]
                                          if lon[0, i]lonmax or lat[0, i]latmax:
                                              continue

                                          trajectory_nans = np.invert(np.isnan(drifter_trajectory))

                                          if type(depth) != str:
                                              if np.min(drifter_trajectory[trajectory_nans])>depth:
                                                  continue

                                              m, depth_idx = self.interpolate(drifter_trajectory[trajectory_nans], drifter_mass[trajectory_nans], depth)

                                          else:
                                              s = status[:, i]
                                              sea_floor = np.where(s == self.seafloor_idx)[0]
                                              if len(sea_floor) >0:
                                                  depth_idx = sea_floor[0]
                                              else:
                                                  depth_idx = len(drifter_trajectory)-1
                                              m = drifter_mass[depth_idx]

                                          if np.isnan(m):
                                              continue

                                          mass_by_biomes += m

                                      logger.debug("Finish summing masses.")
                                      return mass_by_biomes * 10 ** (-15)

- animate_current_3D

                                  def animate_current_3D(self):
                                      """A 3d and 2d current animation of the simulation.
                                      Sizes of particles are relative to their initial mass."""

                                      plt.close()

                                      fig = plt.figure(figsize=self.figsize)
                                      ax = fig.add_subplot(121, projection='3d')
                                      ax2 = fig.add_subplot(122, projection = ccrs.PlateCarree())

                                      logger.debug("Reading required simulation properties.")
                                      z, mass, lon, lat, status = self.obj1.get_properties_from_list([
                                          "z", "mass", "lon", "lat", "status"
                                      ])

                                      seafloor_tracker = {"lon":[], "lat":[], "z":[], "m":[]}
                                      u = self.obj1.get_property("x_sea_water_velocity")
                                      v = self.obj1.get_property("y_sea_water_velocity")
                                      U = np.sqrt(u ** 2 + v ** 2)
                                      U = np.ma.filled(U, np.nan)

                                      logger.debug("Initialize colormap.")
                                      if self.cmap is None:
                                          logger.warning("Colormap hasn't been provided, using jet.")
                                          cmap = mpl.cm.jet
                                      else:
                                          try:
                                              cmap = getattr(mpl.cm, self.cmap)
                                          except AttributeError:
                                              logger.debug("Specified colormap doesn't exist, changing to jet.")
                                              cmap = "jet"

                                      dt = (self.time[0][1]-self.time[0][0]).total_seconds()
                                      idx = np.where(np.invert(np.logical_or(np.isnan(mass[1, :]), np.isnan(z[1, :]))))[0][0]
                                      try:
                                          w0 = abs(self.initial_velocities[0]) * 24 * 3600
                                      except Exception:
                                          w0 = abs((z[1, idx] - z[0, idx]) * (mass[0, idx]/mass[1, idx]) ** (1/6) / dt)*24*3600

                                      lonlim = np.where((lon[0,:]>=self.xlim[0]) & (lon[0,:]<=self.xlim[1]))[0]
                                      latlim = np.where((lat[0,:]>=self.ylim[0]) & (lat[0,:]<=self.ylim[1]))[0]
                                      intersectlim = np.intersect1d(lonlim, latlim)
                                      xlim = [np.min(lon[:, intersectlim][np.invert(np.isnan(lon[:, intersectlim]))]),
                                              np.max(lon[:, intersectlim][np.invert(np.isnan(lon[:, intersectlim]))])]
                                      ylim = [np.min(lat[:, intersectlim][np.invert(np.isnan(lat[:, intersectlim]))]),
                                              np.max(lat[:, intersectlim][np.invert(np.isnan(lat[:, intersectlim]))])]
                                      zlim = [np.min(z[:, intersectlim][np.invert(np.isnan(z[:, intersectlim]))]),
                                              np.max(z[:, intersectlim][np.invert(np.isnan(z[:, intersectlim]))])]
                                      U0 = np.max(U[:, intersectlim][np.invert(np.isnan(U[:, intersectlim]))])
                                      zmin = zlim[0]
                                      k = 1 / mass[0, intersectlim]*100
                                      sc = ax.scatter(lon[0, intersectlim], lat[0, intersectlim], z[0, intersectlim],
                                                      s=mass[0, intersectlim] * k, c = U[0, intersectlim], cmap = cmap, vmin=0, vmax = U0)
                                      sc2 = ax2.scatter(lon[0, intersectlim], lat[0, intersectlim], s=mass[0, intersectlim] * k,
                                                         c=z[0, intersectlim], cmap = cmap, vmin=zmin, vmax = 0)
                                      ax2.quiver(lon[0, :], lat[0, :], u[0, :], v[0, :], color='b', units='xy', scale=1/5)
                                      # ax.scatter(lon[0, intersectlim], lat[0, intersectlim], z[0, intersectlim], s=mass[0, intersectlim] * k, c = U[0, intersectlim], cmap = cmap, vmin=0, vmax = U0)
                                      logger.debug("Setting up axes properties.")
                                      ax.yaxis.labelpad=30
                                      ax.xaxis.labelpad=30
                                      ax.zaxis.labelpad=30
                                      ax.set_xlabel("lon")
                                      ax.set_ylabel("lat")
                                      ax.set_zlabel("z")
                                      ax.set_xlim(*xlim)
                                      ax.set_ylim(*ylim)
                                      ax.set_zlim(*zlim)
                                      fig.suptitle(self.time[0][0])
                                      ax.view_init(elev = 8, azim=45, roll=0)

                                      cb = plt.colorbar(sc, shrink = self.shrink, orientation = "horizontal")
                                      cb.set_label(r"$|\vec{u}_H(t)|\,\mathrm{[m\,s^{-1}]}$")

                                      ax2.add_feature(cartopy.feature.LAND, zorder=2, edgecolor='k', facecolor = "beige")
                                      gl = ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                        linewidth=2, color='k', alpha=0.8, linestyle='--', zorder = 100)
                                      ax2.set_xlabel("lon")
                                      ax2.set_ylabel("lat")
                                      ax2.set_xlim(*xlim)
                                      ax2.set_ylim(*ylim)
                                      ax2.grid()

                                      cb2 = plt.colorbar(sc2, shrink = self.shrink, orientation = "horizontal")
                                      cb2.set_label(r"$z\,\mathrm{[m]}$")

                                      def update(i):
                                          ax.cla()
                                          ax2.cla()
                                          ax.set_xlim(*xlim)
                                          ax.set_ylim(*ylim)
                                          ax.set_zlim(*zlim)
                                          ax2.set_xlim(*xlim)
                                          ax2.set_ylim(*ylim)
                                          ax.yaxis.labelpad=30
                                          ax.xaxis.labelpad=30
                                          ax.zaxis.labelpad=30
                                          ax.set_xlabel("lon")
                                          ax.set_ylabel("lat")
                                          ax.set_zlabel("z")
                                          ax2.set_xlabel("lon")
                                          ax2.set_ylabel("lat")
                                          idx = np.where(np.invert(np.isnan(lon[i, :])))[0]
                                          idx = np.intersect1d(idx, intersectlim)
                                          sf = np.where(status[i, :] == self.seafloor_idx)[0]
                                          sf = np.intersect1d(sf, intersectlim)
                                          for da_tracer in sf:
                                              seafloor_tracker["lon"].append(lon[i, da_tracer])
                                              seafloor_tracker["lat"].append(lat[i, da_tracer])
                                              seafloor_tracker["z"].append(z[i, da_tracer])
                                              seafloor_tracker["m"].append(mass[i, da_tracer]/mass[0, da_tracer]*100)
                                          k = 1 / mass[0, idx]*100
                                          sc = ax.scatter(lon[i, idx], lat[i, idx], z[i, idx], s=mass[i, idx] * k, c = U[i, idx], cmap = cmap, vmin=0, vmax = U0)
                                          sc2 = ax2.scatter(lon[i, idx], lat[i, idx], s=mass[i, idx] * k, c=z[i, idx], cmap = cmap, vmin=zmin, vmax = 0)
                                          ax2.quiver(lon[i, idx], lat[i, idx], u[i, idx], v[i, idx], color='b', units='xy', scale=1/5)
                                          ax.scatter(seafloor_tracker["lon"], seafloor_tracker["lat"], seafloor_tracker["z"], color = "black", s=seafloor_tracker["m"])
                                          ax2.scatter(seafloor_tracker["lon"], seafloor_tracker["lat"], color = "black", s=seafloor_tracker["m"])
                                          ax2.grid()
                                          ax2.add_feature(cartopy.feature.LAND, zorder=2, edgecolor='k', facecolor = "beige")
                                          gl = ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                                                        linewidth=2, color='k', alpha=0.8, linestyle='--', zorder = 100)
                                          cb.update_normal(sc)
                                          cb2.update_normal(sc2)
                                          fig.suptitle(self.time[0][i])
                                          ax.view_init(elev=8, azim=45, roll=0)

                                      frames = len(self.time[0])

                                      logger.debug("Initializing animation.")
                                      anim=animation.FuncAnimation(fig, update, blit=False, frames = frames, interval=100)

                                      if self.outfile is not None:
                                          logger.debug("Adding PillowWriter to animation.")
                                          writer = animation.PillowWriter(fps=10)
                                          logger.debug("Saving animation.")
                                          anim.save(self.outfile, writer = writer)
                                      else:
                                          plt.show()

- vertical_particle_distribution

                                      def vertical_particle_distribution(self):
                                      """Create an interactive slider vertical distribution of particle count or save it as an animation."""

                                      plt.close()
                                      fig = plt.figure(figsize = self.figsize)

                                      if self.outfile is not None:
                                          logger.warning(
                                              "Outfile is given, the plot will be saved as an animation. If you want slider control, pelase remove the outfile argument."
                                              )
                                          mainplot = fig.add_subplot(111)
                                      else:
                                          logger.warning(
                                              "Outfile is not given. If you want to save the plot as an animation, please add the outfile argument."
                                              )
                                          mainplot = fig.add_axes([.15, .3, .8, .5])
                                          sliderax = fig.add_axes([.15, .08, .75, .05])

                                      maxdepth = np.min(self.obj1.get_property("depth")) if self.depth is None else self.depth

                                      z_dist = np.linspace(maxdepth, 0, self.bins)

                                      logger.debug("Retrieving time information.")
                                      maxtime = 0
                                      time_steps = []
                                      for obj in self.objects:
                                          t = obj.data["time"].shape[0]
                                          maxtime = max(t, maxtime)
                                          time_step = datetime.strptime(obj.data.time_step_output, "%H:%M:%S")
                                          time_steps.append(timedelta(hours = time_step.hour, minutes=time_step.minute, seconds=time_step.second).total_seconds() / 3600)
                                      time_steps = np.asarray(np.asarray(time_steps) // np.min(time_steps), dtype=np.int64)

                                      if self.outfile is None:
                                          logger.debug("Setting up slider.")
                                          tslider = Slider(sliderax, 'time', 0, maxtime-1,
                                                          valinit = 0, valfmt='%0.0f', color = "black",
                                                          valstep = 1)

                                      logger.debug("Precalculating distribution.")
                                      particle_count = np.zeros((maxtime, len(self.objects), self.bins))
                                      sea_floor_particles = np.zeros((maxtime, self.bins))

                                      for i, obj in enumerate(self.objects):
                                          z, status = obj.get_properties_from_list([
                                              "z", "status"
                                          ])
                                          time_step = time_steps[i]
                                          for t in range(obj.data["time"].shape[0]):
                                              active = np.invert(np.isnan(z[t, :]))
                                              hist = np.histogram(z[t, active], bins = self.bins,
                                                                                  range=[maxdepth, 0])
                                              particle_count[t*time_step:(t+1)*time_step, i, :] = hist[0]
                                              sea_floor_particles = self.update_sea_floor_particle_count(sea_floor_particles, self.bins,z, active,
                                                                                                         status, hist[1], t, self.seafloor_idx[i], time_step)

                                      if not self.abs:
                                          totnum = np.sum(particle_count[0, :, :])
                                          particle_count /= totnum
                                          sea_floor_particles /= totnum

                                      if self.outfile is not None:
                                          mainplot.barh(z_dist, sea_floor_particles[0, :], color = "black", height = maxdepth/self.bins, label = "seafloor")
                                          for i in range(particle_count.shape[1]):
                                              color = self.color[i] if self.color is not None else None
                                              label = self.labels[i] if self.legend else None
                                              mainplot.barh(z_dist, particle_count[0, i, :],
                                                           left = sea_floor_particles[0, :] + np.sum(particle_count[0, :i, :], axis=0),
                                                           height = maxdepth/self.bins, color = color, label = label)
                                          self.axis_setup(mainplot, "lower right", True)

                                      def update(val):
                                          tindex = int(tslider.val) if self.outfile is None else val
                                          mainplot.cla()
                                          mainplot.barh(z_dist, sea_floor_particles[tindex, :], color = "black", height = maxdepth/self.bins, label = "seafloor")
                                          for i in range(particle_count.shape[1]):
                                              color = self.color[i] if self.color is not None else None
                                              label = self.labels[i] if self.legend else None
                                              mainplot.barh(z_dist, particle_count[tindex, i, :],
                                                           left = sea_floor_particles[tindex, :] + np.sum(particle_count[tindex, :i, :], axis=0),
                                                           height = maxdepth/self.bins, color = color, label = label)
                                          if self.outfile is None: fig.canvas.draw_idle()
                                          self.axis_setup(mainplot, "lower right", True)

                                      if self.outfile is not None:
                                          logger.debug("Initializing animation.")
                                          anim=animation.FuncAnimation(fig, update, blit=False, frames = maxtime, interval=100)
                                          logger.debug("Adding PillowWriter to animation.")
                                          writer = animation.PillowWriter(fps=20)
                                          logger.debug("Saving animation.")
                                          anim.save(self.outfile, writer = writer)
                                      else:
                                          update(0)
                                          tslider.on_changed(update)
                                          plt.show()

- update_sea_floor_particle_count

                                  @staticmethod
                                  @jit(nopython=True)
                                  def update_sea_floor_particle_count(sea_floor_particles, bins, z, active, status, dist, t, sf, time_step):
                                      for j in range(bins):
                                          if j == 0:
                                              zj = z[t, active]<= dist[j]
                                          else:
                                              zj = np.logical_and(z[t, active]<= dist[j], z[t, active]> dist[j-1])
                                          sea_floor_particles[t*time_step:, j] += len(np.where(status[t, active][zj] == sf)[0])
                                      return sea_floor_particles

- axis_setup

                                  def axis_setup(self, ax, legendloc = None, grid=False):
                                      if legendloc!=False:
                                          if self.legend:ax.legend(loc = legendloc)
                                      if self.xlim is not None:ax.set_xlim(*self.xlim)
                                      if self.ylim is not None: ax.set_ylim(*self.ylim)
                                      if self.xlabel is not None: ax.set_xlabel(f"{self.xlabel}")
                                      if self.ylabel is not None: ax.set_ylabel(f"{self.ylabel}")
                                      if self.title is not None: ax.set_title(self.title, fontweight = self.fontweight)
                                      if grid:ax.grid()