docs: def create_3D_grid(self, catalogue, completeness_table, t_f=1.0,
                         mag_inc=0.1):
          '''
          Counts the earthquakes observed in a three dimensional grid
  
  
          :param catalogue:
              Instance of the openquake.hmtk.seismicity.catalogue.Catalogue class
              catalogue.data dictionary containing the following -
              'year' - numpy.ndarray vector of years
              'longitude' - numpy.ndarray vector of longitudes
              'latitude' - numpy.ndarray vector of latitudes
              'depth' - numpy.ndarray vector of depths
  
          :param np.ndarray completeness_table:
              Completeness of the catalogue assuming evenly spaced magnitudes
              from most recent bin to oldest bin [year, magnitude]
  
          :param float t_f:
              Weichert adjustment factor
  
          :param float mag_inc:
              Increment of the completeness magnitude (rendered 0.1)
  
          :returns:
             Three-dimensional spatial grid of observed rates (or two dimensional
             if only one depth layer is considered)
  
          '''
          x_bins = np.arange(self.grid_limits['xmin'],
                             self.grid_limits['xmax'],
                             self.grid_limits['xspc'])
          if x_bins[-1] < self.grid_limits['xmax']:
              x_bins = np.hstack([x_bins, x_bins[-1] + self.grid_limits['xspc']])
  
          y_bins = np.arange(self.grid_limits['ymin'],
                             self.grid_limits['ymax'],
                             self.grid_limits['yspc'])
          if y_bins[-1] < self.grid_limits['ymax']:
              y_bins = np.hstack([y_bins, y_bins[-1] + self.grid_limits['yspc']])
  
          z_bins = np.arange(self.grid_limits['zmin'],
                             self.grid_limits['zmax'] + self.grid_limits['zspc'],
                             self.grid_limits['zspc'])
  
          if z_bins[-1] < self.grid_limits['zmax']:
              z_bins = np.hstack([z_bins, z_bins[-1] + self.grid_limits['zspc']])
  
          # Define centre points of grid cells
          gridx, gridy = np.meshgrid((x_bins[1:] + x_bins[:-1]) / 2.,
                                     (y_bins[1:] + y_bins[:-1]) / 2.)
  
          n_x, n_y = np.shape(gridx)
          gridx = np.reshape(gridx, [n_x * n_y, 1])
          gridy = np.reshape(np.flipud(gridy), [n_x * n_y, 1])
  
          # Only one depth range
          idx = np.logical_and(catalogue.data['depth'] >= z_bins[0],
                               catalogue.data['depth'] < z_bins[1])
          mid_depth = (z_bins[0] + z_bins[1]) / 2.
  
          data_grid = np.column_stack([
              gridx,
              gridy,
              mid_depth * np.ones(n_x * n_y, dtype=float),
              self.create_2D_grid_simple(catalogue.data['longitude'][idx],
                                         catalogue.data['latitude'][idx],
                                         catalogue.data['year'][idx],
                                         catalogue.data['magnitude'][idx],
                                         completeness_table,
                                         t_f,
                                         mag_inc)])
  
          if len(z_bins) < 3:
              # Only one depth range
              self.data = data_grid
              return
  
          # Multiple depth layers - append to grid
          for iloc in range(1, len(z_bins) - 1):
              idx = np.logical_and(catalogue.data['depth'] >= z_bins[iloc],
                                   catalogue.data['depth'] < z_bins[iloc + 1])
              mid_depth = (z_bins[iloc] + z_bins[iloc + 1]) / 2.
  
              temp_grid = np.column_stack([
                  gridx,
                  gridy,
                  mid_depth * np.ones(n_x * n_y, dtype=float),
                  self.create_2D_grid_simple(catalogue.data['longitude'][idx],
                                             catalogue.data['latitude'][idx],
                                             catalogue.data['year'][idx],
                                             catalogue.data['magnitude'][idx],
                                             completeness_table,
                                             t_f,
                                             mag_inc)])
  
              data_grid = np.vstack([data_grid, temp_grid])
          self.data = data_grid