diff --git a/atlaselectrophysiology/compare_alignments.py b/atlaselectrophysiology/compare_alignments.py
index d09de30..e4e4298 100644
--- a/atlaselectrophysiology/compare_alignments.py
+++ b/atlaselectrophysiology/compare_alignments.py
@@ -4,7 +4,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 from pathlib import Path
 # Instantiate brain atlas and one
 brain_atlas = atlas.AllenAtlas(25)
diff --git a/atlaselectrophysiology/get_scale_factor.py b/atlaselectrophysiology/get_scale_factor.py
index a2ceb73..90fa288 100644
--- a/atlaselectrophysiology/get_scale_factor.py
+++ b/atlaselectrophysiology/get_scale_factor.py
@@ -12,7 +12,7 @@
 import pandas as pd
 import matplotlib.pyplot as plt
 import seaborn as sns
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 
 
 # Instantiate brain atlas and one
diff --git a/atlaselectrophysiology/load_data.py b/atlaselectrophysiology/load_data.py
index 05e7e5d..26c5a71 100644
--- a/atlaselectrophysiology/load_data.py
+++ b/atlaselectrophysiology/load_data.py
@@ -3,8 +3,10 @@
 from datetime import datetime
 import ibllib.pipes.histology as histology
 from neuropixel import trace_header
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 from ibllib.qc.alignment_qc import AlignmentQC
+from iblutil.numerical import ismember
+from iblutil.util import Bunch
 from one.api import ONE
 from one.remote import aws
 from pathlib import Path
@@ -259,16 +261,29 @@ def get_data(self):
         try:
             data['spikes'] = self.one.load_object(self.eid, 'spikes', collection=self.probe_collection,
                                                   attribute=['depths', 'amps', 'times', 'clusters'])
-            data['spikes']['exists'] = True
 
             data['clusters'] = self.one.load_object(self.eid, 'clusters', collection=self.probe_collection,
                                                     attribute=['metrics', 'peakToTrough', 'waveforms', 'channels'])
+
+            # Remove low firing rate clusters
+            min_firing_rate = 50. / 3600.
+            clu_idx = data['clusters'].metrics.firing_rate > min_firing_rate
+            data['clusters'] = Bunch({k: v[clu_idx] for k, v in data['clusters'].items()})
+            spike_idx, ib = ismember(data['spikes'].clusters, data['clusters'].metrics.index)
+            data['clusters'].metrics.reset_index(drop=True, inplace=True)
+            data['spikes'] = Bunch({k: v[spike_idx] for k, v in data['spikes'].items()})
+            data['spikes'].clusters = data['clusters'].metrics.index[ib].astype(np.int32)
+
+            data['spikes']['exists'] = True
             data['clusters']['exists'] = True
 
             data['channels'] = self.one.load_object(self.eid, 'channels', collection=self.probe_collection,
                                                     attribute=['rawInd', 'localCoordinates'])
             data['channels']['exists'] = True
 
+            # Set low firing rate clusters to bad
+
+
         except alf.exceptions.ALFObjectNotFound:
             logger.error(f'Could not load spike sorting for probe insertion {self.probe_id}, GUI'
                          f' will not work')
diff --git a/atlaselectrophysiology/load_data_local.py b/atlaselectrophysiology/load_data_local.py
index 78c5a6a..24e29a3 100644
--- a/atlaselectrophysiology/load_data_local.py
+++ b/atlaselectrophysiology/load_data_local.py
@@ -1,7 +1,7 @@
 import logging
 import numpy as np
 from datetime import datetime
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 from pathlib import Path
 import one.alf.io as alfio
 from one import alf
diff --git a/atlaselectrophysiology/plot_data.py b/atlaselectrophysiology/plot_data.py
index 97043dc..e3dd517 100644
--- a/atlaselectrophysiology/plot_data.py
+++ b/atlaselectrophysiology/plot_data.py
@@ -1,6 +1,6 @@
 from matplotlib import cm
 import numpy as np
-from brainbox.processing import bincount2D
+from iblutil.numerical import bincount2D
 from brainbox.io.spikeglx import Streamer
 from brainbox.population.decode import xcorr
 from brainbox.task import passive
@@ -596,7 +596,7 @@ def get_autocorr(self, clust_idx):
         autocorr = xcorr(self.data['spikes']['times'][idx], self.data['spikes']['clusters'][idx],
                          AUTOCORR_BIN_SIZE, AUTOCORR_WIN_SIZE)
 
-        return autocorr[0, 0, :], self.clust_id[clust_idx]
+        return autocorr[0, 0, :], self.data['clusters'].metrics.cluster_id[self.clust_id[clust_idx]]
 
     def get_template_wf(self, clust_idx):
         template_wf = (self.data['clusters']['waveforms'][self.clust_id[clust_idx], :, 0])
diff --git a/atlasview/atlasview.py b/atlasview/atlasview.py
index 69b2675..07fee99 100644
--- a/atlasview/atlasview.py
+++ b/atlasview/atlasview.py
@@ -2,7 +2,7 @@
 TopView is the main Widget with the related ControllerTopView Class
 There are several SliceView windows (sagittal, coronal, possibly tilted etc...) that each have
 a SliceController object
-The underlying data model object is an ibllib.atlas.AllenAtlas object
+The underlying data model object is an iblatlas.atlas.AllenAtlas object
 
     TopView(QMainWindow)
     ControllerTopView(PgImageController)
@@ -20,7 +20,7 @@
 import pyqtgraph as pg
 import matplotlib
 
-from ibllib.atlas import AllenAtlas
+from iblatlas.atlas import AllenAtlas
 import qt
 
 
@@ -85,7 +85,7 @@ def add_image_layer(self, **kwargs):
         """
         :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None,
         'opacity': 1.0}
-        :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+        :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
         :return:
         """
         self.ctrl.fig_sagittal.add_image_layer(**kwargs)
@@ -160,7 +160,7 @@ def add_image_layer(self, **kwargs):
         """
         :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None,
         'opacity': 1.0}
-        :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+        :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
         :return:
         """
         il = ImageLayer(**kwargs)
@@ -354,7 +354,7 @@ class ImageLayer:
     Class for keeping track of image layers.
     :param image_item
     :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None, 'opacity': 1.0}
-    :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+    :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
     :param
     """
     image_item: pg.ImageItem = field(default_factory=pg.ImageItem)
diff --git a/ephysfeatures/features_across_region.py b/ephysfeatures/features_across_region.py
index 7111f84..929e9a0 100644
--- a/ephysfeatures/features_across_region.py
+++ b/ephysfeatures/features_across_region.py
@@ -7,7 +7,7 @@
 import pandas as pd
 import numpy as np
 
-from ibllib.atlas import AllenAtlas
+from iblatlas.atlas import AllenAtlas
 import atlaselectrophysiology.ColorBar as cb
 from ibllib.pipes.ephys_alignment import EphysAlignment
 from atlaselectrophysiology.AdaptedAxisItem import replace_axis
diff --git a/histology/atlas_mpl.py b/histology/atlas_mpl.py
index f51ec01..086f2bd 100644
--- a/histology/atlas_mpl.py
+++ b/histology/atlas_mpl.py
@@ -5,7 +5,7 @@
 
 from iblapps import qt
 from iblapps.qt_matplotlib import BaseMplCanvas
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 
 # Make sure that we are using QT5
 matplotlib.use('Qt5Agg')
diff --git a/launch_phy/phy_launcher.py b/launch_phy/phy_launcher.py
index b0d93f5..6ae4aa7 100644
--- a/launch_phy/phy_launcher.py
+++ b/launch_phy/phy_launcher.py
@@ -23,7 +23,7 @@ def launch_phy(probe_name=None, eid=None, pid=None, subj=None, date=None, sess_n
     # -------------------- #
 
     from one.api import ONE
-    from ibllib.atlas import AllenAtlas
+    from iblatlas.atlas import AllenAtlas
     from brainbox.io.one import SpikeSortingLoader
     from ibllib.io import spikeglx
     one = one or ONE(base_url='https://openalyx.internationalbrainlab.org')
diff --git a/needles2/probe_model.py b/needles2/probe_model.py
index 0552da2..feeec00 100644
--- a/needles2/probe_model.py
+++ b/needles2/probe_model.py
@@ -1,5 +1,6 @@
 import time
 import copy
+import re
 import numpy as np
 import pandas as pd
 from scipy.signal import fftconvolve
@@ -7,7 +8,8 @@
 from iblutil.numerical import ismember
 
 from ibllib.pipes import histology
-from ibllib.atlas import AllenAtlas, atlas
+from iblatlas.atlas import AllenAtlas
+from iblatlas import atlas
 from neuropixel import TIP_SIZE_UM, trace_header
 from ibllib.pipes.ephys_alignment import EphysAlignment
 from neurodsp.utils import fcn_cosine
@@ -317,7 +319,8 @@ def compute_coverage(self, trajs, dist_fcn=[50, 100], limit=True, coverage=None,
         """
 
         ba = self.ba
-        ACTIVE_LENGTH_UM = 3.84 * 1e3  # This is the length of the NP1 probe with electrodes
+        ACTIVE_LENGTH_UM_1shank = 3.84 * 1e3  # This is the length of the NP1 probe with electrodes
+        ACTIVE_LENGTH_UM_4shank = 705  # This is the length of the NP2 4 shank probe with electrodes
         MAX_DIST_UM = dist_fcn[1]  # max distance around the probe to be searched for
 
         # Covered_length_um are two values which indicate on the path from tip to entry of a given insertion
@@ -325,8 +328,11 @@ def compute_coverage(self, trajs, dist_fcn=[50, 100], limit=True, coverage=None,
         # Note that the second value is negative, because the covered regions extends beyond the tip of the probe
         # We multiply by sqrt(2) to translate the radius given by dist_fcn[1] into the side length of a square
         # that is contained in the circle with radius dist_fcn[1]
-        covered_length_um = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM + MAX_DIST_UM * np.sqrt(2),
-                                                    -MAX_DIST_UM * np.sqrt(2)])
+        covered_length_um_1shank = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM_1shank + MAX_DIST_UM * np.sqrt(2),
+                                                           -MAX_DIST_UM * np.sqrt(2)])
+        covered_length_um_4shank = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM_4shank + MAX_DIST_UM * np.sqrt(2),
+                                                           -MAX_DIST_UM * np.sqrt(2)])
+
 
         # Horizontal slice of voxels to be considered around each trajectory is only dependent on MAX_DIST_UM
         # and the voxel resolution, so can be defined here. We translate max dist in voxels and add 1 for safety
@@ -361,6 +367,19 @@ def crawl_up_from_tip(ins, covered_length):
             if len(trajs) > 20 and self.verbose is True:
                 if p % 20 == 0:
                     print(p / len(trajs))
+
+            # Here we find out if this is trajectory is from a 1shank or 4shank probe
+            # In an ideal world we would read in the metadata and find this out, but that would require
+            # downloading this dataset for all trajectories. Instead we go based on naming convention. We know
+            # probes with the label probe00a, probe00b etc. are 4 shank probes that have been split.
+            pname = traj['probe_name']
+            if len(re.findall("[a-d]", pname[-1])) == 1:
+                covered_length_um = covered_length_um_4shank
+                ACTIVE_LENGTH_UM = ACTIVE_LENGTH_UM_4shank
+            else:
+                covered_length_um = covered_length_um_1shank
+                ACTIVE_LENGTH_UM = ACTIVE_LENGTH_UM_1shank
+
             # Get one trajectory from the list and create an insertion in the brain atlas
             # x and y coordinates of entry are translated to the atlas voxel space
             # z is locked to surface of the brain at these x,y coordinates (disregarding actual z value of trajectory)
@@ -513,7 +532,8 @@ def compute_coverage_dict(self, trajs, dist_fcn=[50, 100], limit=True, coverage=
         """
 
         ba = self.ba
-        ACTIVE_LENGTH_UM = 3.84 * 1e3  # This is the length of the NP1 probe with electrodes
+        ACTIVE_LENGTH_UM_1shank = 3.84 * 1e3  # This is the length of the NP1 probe with electrodes
+        ACTIVE_LENGTH_UM_4shank = 705  # This is the length of the NP2 4 shank probe with electrodes
         MAX_DIST_UM = dist_fcn[1]  # max distance around the probe to be searched for
 
         # Covered_length_um are two values which indicate on the path from tip to entry of a given insertion
@@ -521,8 +541,10 @@ def compute_coverage_dict(self, trajs, dist_fcn=[50, 100], limit=True, coverage=
         # Note that the second value is negative, because the covered regions extends beyond the tip of the probe
         # We multiply by sqrt(2) to translate the radius given by dist_fcn[1] into the side length of a square
         # that is contained in the circle with radius dist_fcn[1]
-        covered_length_um = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM + MAX_DIST_UM * np.sqrt(2),
-                                                    -MAX_DIST_UM * np.sqrt(2)])
+        covered_length_um_1shank = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM_1shank + MAX_DIST_UM * np.sqrt(2),
+                                                           -MAX_DIST_UM * np.sqrt(2)])
+        covered_length_um_4shank = TIP_SIZE_UM + np.array([ACTIVE_LENGTH_UM_4shank + MAX_DIST_UM * np.sqrt(2),
+                                                           -MAX_DIST_UM * np.sqrt(2)])
 
         # Horizontal slice of voxels to be considered around each trajectory is only dependent on MAX_DIST_UM
         # and the voxel resolution, so can be defined here. We translate max dist in voxels and add 1 for safety
@@ -542,10 +564,25 @@ def crawl_up_from_tip(ins, covered_length):
             if len(trajs) > 20 and self.verbose is True:
                 if p % 20 == 0:
                     print(p / len(trajs))
-            # Get one trajectory from the list and create an insertion in the brain atlas
+
+            # Get one trajectory from the list and
+            traj = trajs[p]
+
+            # Here we find out if this is trajectory is from a 1shank or 4shank probe
+            # In an ideal world we would read in the metadata and find this out, but that would require
+            # downloading this dataset for all trajectories. Instead we go based on naming convention. We know
+            # probes with the label probe00a, probe00b etc. are 4 shank probes that have been split.
+            pname = traj['probe_name']
+            if len(re.findall("[a-d]", pname[-1])) == 1:
+                covered_length_um = covered_length_um_4shank
+                ACTIVE_LENGTH_UM = ACTIVE_LENGTH_UM_4shank
+            else:
+                covered_length_um = covered_length_um_1shank
+                ACTIVE_LENGTH_UM = ACTIVE_LENGTH_UM_1shank
+
+            # Create an insertion in the brain atlas
             # x and y coordinates of entry are translated to the atlas voxel space
             # z is locked to surface of the brain at these x,y coordinates (disregarding actual z value of trajectory)
-            traj = trajs[p]
             ins = atlas.Insertion.from_dict(traj, brain_atlas=ba)
             # Don't use probes that have same entry and tip, something is wrong
             set_nan = False
diff --git a/needles2/run_needles2.py b/needles2/run_needles2.py
index f7ba1d5..2f4337f 100644
--- a/needles2/run_needles2.py
+++ b/needles2/run_needles2.py
@@ -9,7 +9,7 @@
 import pyqtgraph as pg
 import matplotlib
 
-from ibllib.atlas import AllenAtlas, Insertion
+from iblatlas.atlas import AllenAtlas, Insertion
 
 import qt
 
@@ -1219,7 +1219,7 @@ def add_image_layer(self, idx=None, **kwargs):
         :param name: name of the image item to keep track of layers
         :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None,
         'opacity': 1.0}
-        :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+        :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
         :return:
         """
         il = ImageLayer(**kwargs)
@@ -1342,7 +1342,7 @@ def add_scatter(self, idx=None, **kwargs):
         :param name: name of the image item to keep track of layers
         :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None,
         'opacity': 1.0}
-        :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+        :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
         :return:
         """
         sc = ScatterLayer(**kwargs)
@@ -1429,7 +1429,7 @@ class ImageLayer:
     Class for keeping track of image layers.
     :param image_item
     :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None, 'opacity': 1.0}
-    :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+    :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
     :param
     """
     name: str = field(default='base')
diff --git a/needles2/spike_features.py b/needles2/spike_features.py
index 6c40826..a769168 100644
--- a/needles2/spike_features.py
+++ b/needles2/spike_features.py
@@ -458,7 +458,7 @@ class SpikeSortFigures:
     Class for keeping track of image layers.
     :param image_item
     :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None, 'opacity': 1.0}
-    :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+    :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
     :param
     """
     name: str = field(default='base')
@@ -473,7 +473,7 @@ class SpikeSortRawData:
     Class for keeping track of image layers.
     :param image_item
     :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None, 'opacity': 1.0}
-    :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+    :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
     :param
     """
     name: str = field(default='base')
@@ -488,7 +488,7 @@ class SpikeSortPlotData:
     Class for keeping track of image layers.
     :param image_item
     :param pg_kwargs: pyqtgraph setImage arguments: {'levels': None, 'lut': None, 'opacity': 1.0}
-    :param slice_kwargs: ibllib.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
+    :param slice_kwargs: iblatlas.atlas.slice arguments: {'volume': 'image', 'mode': 'clip'}
     :param
     """
     name: str = field(default='base')
diff --git a/tests/test_alignment_qc_gui.py b/tests/test_alignment_qc_gui.py
index 7977bd5..e61e1bb 100644
--- a/tests/test_alignment_qc_gui.py
+++ b/tests/test_alignment_qc_gui.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 from one.api import ONE
-from ibllib.atlas import AllenAtlas
+from iblatlas.atlas import AllenAtlas
 from atlaselectrophysiology.load_data import LoadData
 from ibllib.pipes.ephys_alignment import EphysAlignment
 from ibllib.pipes.misc import create_alyx_probe_insertions
diff --git a/viewspikes/examples_local.py b/viewspikes/examples_local.py
index dfe3815..bd2ab40 100644
--- a/viewspikes/examples_local.py
+++ b/viewspikes/examples_local.py
@@ -89,7 +89,7 @@
  'clusters.channels',
  'clusters.mlapdv']
 
-from ibllib.atlas import atlas
+from iblatlas import atlas
 from ibllib.pipes import histology
 from ibllib.ephys import neuropixel
 
diff --git a/viewspikes/plots.py b/viewspikes/plots.py
index 611e92d..2601d34 100644
--- a/viewspikes/plots.py
+++ b/viewspikes/plots.py
@@ -4,7 +4,7 @@
 import scipy.signal
 import pyqtgraph as pg
 
-import ibllib.atlas as atlas
+import iblatlas.atlas as atlas
 from neuropixel import SITES_COORDINATES
 from ibllib.pipes.ephys_alignment import EphysAlignment
 from ibllib.plots import wiggle, color_cycle