docs: def plot_magnitude_time_density(
          catalogue, mag_int, time_int, completeness=None,
          normalisation=False, logscale=True, bootstrap=None, xlim=[], ylim=[],
          filename=None, figure_size=(8, 6), filetype='png', dpi=300, ax=None):
      """
      Creates a plot of magnitude-time density
  
      :param catalogue:
          Earthquake catalogue as instance of :class:
          openquake.hmtk.seismicity.catalogue.Catalogue
      :param float mag_int:
          Width of the histogram for the magnitude bins
      :param float time_int:
          Width of the histogram for the time bin (in decimal years)
      :param bool normalisation:
          Normalise the histogram to give output as PMF (True) or count (False)
      :param int bootstrap:
          To sample magnitude and depth uncertainties choose number of samples
      """
      if ax is None:
          fig, ax = plt.subplots(figsize=figure_size)
      else:
          fig = ax.get_figure()
  
      # Create the magnitude bins
      if isinstance(mag_int, (np.ndarray, list)):
          mag_bins = mag_int
      else:
          mag_bins = np.arange(
              np.min(catalogue.data['magnitude']),
              np.max(catalogue.data['magnitude']) + mag_int / 2.,
              mag_int)
      # Creates the time bins
      if isinstance(time_int, (np.ndarray, list)):
          time_bins = time_int
      else:
          time_bins = np.arange(
              float(np.min(catalogue.data['year'])),
              float(np.max(catalogue.data['year'])) + 1.,
              float(time_int))
      # Get magnitude-time distribution
      mag_time_dist = catalogue.get_magnitude_time_distribution(
          mag_bins,
          time_bins,
          normalisation,
          bootstrap)
      # Get smallest non-zero value
      vmin_val = np.min(mag_time_dist[mag_time_dist > 0.])
      # Create plot
      if logscale:
          norm_data = LogNorm(vmin=vmin_val, vmax=np.max(mag_time_dist))
      else:
          if normalisation:
              norm_data = Normalize(vmin=vmin_val, vmax=np.max(mag_time_dist))
          else:
              norm_data = Normalize(vmin=1.0, vmax=np.max(mag_time_dist))
  
      im = ax.pcolor(time_bins[:-1],
                     mag_bins[:-1],
                     mag_time_dist.T,
                     norm=norm_data)
      ax.set_xlabel('Time (year)')
      ax.set_ylabel('Magnitude')
      if len(xlim) == 2:
          ax.set_xlim(xlim[0], xlim[1])
      else:
          ax.set_xlim(time_bins[0], time_bins[-1])
      if len(ylim) == 2:
          ax.set_ylim(ylim[0], ylim[1])
      else:
          ax.set_ylim(mag_bins[0], mag_bins[-1] + (mag_bins[-1] - mag_bins[-2]))
      # Fix the title
      if normalisation:
          fig.colorbar(im, label='Event Density', shrink=0.9, ax=ax)
      else:
          fig.colorbar(im, label='Event Count', shrink=0.9, ax=ax)
      ax.grid(True)
      # Plot completeness
      if completeness is not None:
          _plot_completeness(ax, completeness, time_bins[0], time_bins[-1])
  
      _save_image(fig, filename, filetype, dpi)