diff --git a/doc/source/cookbook/halo_plotting.py b/doc/source/cookbook/halo_plotting.py
deleted file mode 100644
index d701fdd9b27..00000000000
--- a/doc/source/cookbook/halo_plotting.py
+++ /dev/null
@@ -1,12 +0,0 @@
-import yt
-
-# Load the dataset
-ds = yt.load("Enzo_64/RD0006/RedshiftOutput0006")
-
-# Load the halo list from a rockstar output for this dataset
-halos = yt.load("rockstar_halos/halos_0.0.bin")
-
-# Create a projection with the halos overplot on top
-p = yt.ProjectionPlot(ds, "x", "density")
-p.annotate_halos(halos)
-p.save()
diff --git a/doc/source/reference/api/api.rst b/doc/source/reference/api/api.rst
index e0f4a692ab8..393855875a4 100644
--- a/doc/source/reference/api/api.rst
+++ b/doc/source/reference/api/api.rst
@@ -683,7 +683,6 @@ See also :ref:`callbacks`.
    ~yt.visualization.plot_modifications.ContourCallback
    ~yt.visualization.plot_modifications.CuttingQuiverCallback
    ~yt.visualization.plot_modifications.GridBoundaryCallback
-   ~yt.visualization.plot_modifications.HaloCatalogCallback
    ~yt.visualization.plot_modifications.ImageLineCallback
    ~yt.visualization.plot_modifications.LinePlotCallback
    ~yt.visualization.plot_modifications.MagFieldCallback
diff --git a/doc/source/visualizing/callbacks.rst b/doc/source/visualizing/callbacks.rst
index ec5898cafbb..5ec74c46bc3 100644
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -390,60 +390,6 @@ Overplot Cell Edges
    slc.annotate_cell_edges()
    slc.save()
 
-.. _annotate-halos:
-
-Overplot Halo Annotations
-~~~~~~~~~~~~~~~~~~~~~~~~~
-
-.. function:: annotate_halos(self, halo_catalog, circle_args=None, \
-                             width=None, annotate_field=None, \
-                             radius_field='virial_radius', \
-                             center_field_prefix="particle_position", \
-                             text_args=None, factor=1.0)
-
-   (This is a proxy for
-   :class:`~yt.visualization.plot_modifications.HaloCatalogCallback`.)
-
-   Accepts a :class:`~yt_astro_analysis.halo_analysis.halo_catalog.HaloCatalog`
-   and plots a circle at the location of each halo with the radius of the
-   circle corresponding to the virial radius of the halo. Also accepts a
-   :ref:`loaded halo catalog dataset <halo-catalog-data>` or a data
-   container from a halo catalog dataset. If ``width`` is set
-   to None (default) all halos are plotted, otherwise it accepts a tuple in
-   the form (1.0, ‘Mpc’) to only display halos that fall within a slab with
-   width ``width`` centered on the center of the plot data.  The appearance of
-   the circles can be changed with the circle_kwargs dictionary, which is
-   supplied to the Matplotlib patch Circle.  One can label each of the halos
-   with the annotate_field, which accepts a field contained in the halo catalog
-   to add text to the plot near the halo (example: ``annotate_field=
-   'particle_mass'`` will write the halo mass next to each halo, whereas
-   ``'particle_identifier'`` shows the halo number). The size of the circles is
-   found from the field ``radius_field`` which is ``'virial_radius'`` by
-   default. If another radius has been found as part of your halo analysis
-   workflow, you can save that field and use it as the ``radius_field`` to
-   change the size of the halos. The position of each halo is determined using
-   ``center_field_prefix`` in the following way. If ``'particle_position'``
-   is the value of ``center_field_prefix`` as is the default, the x value of
-   the halo position is stored in the field ``'particle_position_x'``, y is
-   ``'particle_position_y'``, and z is ``'particle_position_z'``. If you have
-   stored another set of coordinates for each halo as part of your halo
-   analysis as fields such as ``'halo_position_x'``, you can use these fields
-   to determine halo position by passing ``'halo_position'`` to
-   ``center_field_prefix``. font_kwargs contains the arguments controlling the
-   text appearance of the annotated field. Factor is the number the virial
-   radius is multiplied by for plotting the circles. Ex: ``factor=2.0`` will
-   plot circles with twice the radius of each halo virial radius.
-
-.. python-script::
-
-   import yt
-
-   data_ds = yt.load("Enzo_64/RD0006/RedshiftOutput0006")
-   halos_ds = yt.load("rockstar_halos/halos_0.0.bin")
-
-   prj = yt.ProjectionPlot(data_ds, "z", ("gas", "density"))
-   prj.annotate_halos(halos_ds, annotate_field="particle_identifier")
-   prj.save()
 
 .. _annotate-image-line:
 
diff --git a/yt/visualization/plot_modifications.py b/yt/visualization/plot_modifications.py
index 9fc3e140e78..6be50c71e39 100644
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -6,10 +6,8 @@
 import matplotlib
 import numpy as np
 
-from yt.data_objects.data_containers import YTDataContainer
 from yt.data_objects.level_sets.clump_handling import Clump
 from yt.data_objects.selection_objects.cut_region import YTCutRegion
-from yt.data_objects.static_output import Dataset
 from yt.frontends.ytdata.data_structures import YTClumpContainer
 from yt.funcs import is_sequence, mylog, validate_width_tuple
 from yt.geometry.geometry_handler import is_curvilinear
@@ -25,7 +23,6 @@
     pixelize_off_axis_cartesian,
 )
 from yt.utilities.math_utils import periodic_ray
-from yt.utilities.on_demand_imports import NotAModule
 from yt.visualization.image_writer import apply_colormap
 
 callback_registry = {}
@@ -1784,252 +1781,6 @@ def __call__(self, plot):
         super().__call__(plot)
 
 
-class HaloCatalogCallback(PlotCallback):
-    """
-    Plots circles at the locations of all the halos
-    in a halo catalog with radii corresponding to the
-    virial radius of each halo.
-
-    Note, this functionality requires the yt_astro_analysis
-    package. See https://yt-astro-analysis.readthedocs.io/
-    for more information.
-
-    Parameters
-    ----------
-    halo_catalog : Dataset, DataContainer,
-                   or ~yt.analysis_modules.halo_analysis.halo_catalog.HaloCatalog
-        The object containing halos to be overplotted. This can
-        be a HaloCatalog object, a loaded halo catalog dataset,
-        or a data container from a halo catalog dataset.
-    circle_args : list
-        Contains the arguments controlling the
-        appearance of the circles, supplied to the
-        Matplotlib patch Circle.
-    width : tuple
-        The width over which to select halos to plot,
-        useful when overplotting to a slice plot. Accepts
-        a tuple in the form (1.0, 'Mpc').
-    annotate_field : str
-        A field contained in the
-        halo catalog to add text to the plot near the halo.
-        Example: annotate_field = 'particle_mass' will
-        write the halo mass next to each halo.
-    radius_field : str
-        A field contained in the halo
-        catalog to set the radius of the circle which will
-        surround each halo. Default: 'virial_radius'.
-    center_field_prefix : str
-        Accepts a field prefix which will
-        be used to find the fields containing the coordinates
-        of the center of each halo. Ex: 'particle_position'
-        will result in the fields 'particle_position_x' for x
-        'particle_position_y' for y, and 'particle_position_z'
-        for z. Default: 'particle_position'.
-    text_args : dict
-        Contains the arguments controlling the text
-        appearance of the annotated field.
-    factor : float
-        A number the virial radius is multiplied by for
-        plotting the circles. Ex: factor = 2.0 will plot
-        circles with twice the radius of each halo virial radius.
-
-    Examples
-    --------
-
-    >>> import yt
-    >>> dds = yt.load("Enzo_64/DD0043/data0043")
-    >>> hds = yt.load("rockstar_halos/halos_0.0.bin")
-    >>> p = yt.ProjectionPlot(
-    ...     dds, "x", ("gas", "density"), weight_field=("gas", "density")
-    ... )
-    >>> p.annotate_halos(hds)
-    >>> p.save()
-
-    >>> # plot a subset of all halos
-    >>> import yt
-    >>> dds = yt.load("Enzo_64/DD0043/data0043")
-    >>> hds = yt.load("rockstar_halos/halos_0.0.bin")
-    >>> # make a region half the width of the box
-    >>> dregion = dds.box(
-    ...     dds.domain_center - 0.25 * dds.domain_width,
-    ...     dds.domain_center + 0.25 * dds.domain_width,
-    ... )
-    >>> hregion = hds.box(
-    ...     hds.domain_center - 0.25 * hds.domain_width,
-    ...     hds.domain_center + 0.25 * hds.domain_width,
-    ... )
-    >>> p = yt.ProjectionPlot(
-    ...     dds,
-    ...     "x",
-    ...     ("gas", "density"),
-    ...     weight_field=("gas", "density"),
-    ...     data_source=dregion,
-    ...     width=0.5,
-    ... )
-    >>> p.annotate_halos(hregion)
-    >>> p.save()
-
-    >>> # plot halos from a HaloCatalog
-    >>> import yt
-    >>> from yt.extensions.astro_analysis.halo_analysis.api import HaloCatalog
-    >>> dds = yt.load("Enzo_64/DD0043/data0043")
-    >>> hds = yt.load("rockstar_halos/halos_0.0.bin")
-    >>> hc = HaloCatalog(data_ds=dds, halos_ds=hds)
-    >>> p = yt.ProjectionPlot(
-    ...     dds, "x", ("gas", "density"), weight_field=("gas", "density")
-    ... )
-    >>> p.annotate_halos(hc)
-    >>> p.save()
-
-    """
-
-    _type_name = "halos"
-    region = None
-    _descriptor = None
-    _supported_geometries = ("cartesian", "spectral_cube")
-
-    def __init__(
-        self,
-        halo_catalog,
-        circle_args=None,
-        circle_kwargs=None,
-        width=None,
-        annotate_field=None,
-        radius_field="virial_radius",
-        center_field_prefix="particle_position",
-        text_args=None,
-        font_kwargs=None,
-        factor=1.0,
-    ):
-
-        try:
-            from yt_astro_analysis.halo_analysis.api import HaloCatalog
-        except ImportError:
-            HaloCatalog = NotAModule("yt_astro_analysis")
-
-        PlotCallback.__init__(self)
-        def_circle_args = {"edgecolor": "white", "facecolor": "None"}
-        def_text_args = {"color": "white"}
-
-        if isinstance(halo_catalog, YTDataContainer):
-            self.halo_data = halo_catalog
-        elif isinstance(halo_catalog, Dataset):
-            self.halo_data = halo_catalog.all_data()
-        elif isinstance(halo_catalog, HaloCatalog):
-            if halo_catalog.data_source.ds == halo_catalog.halos_ds:
-                self.halo_data = halo_catalog.data_source
-            else:
-                self.halo_data = halo_catalog.halos_ds.all_data()
-        else:
-            raise RuntimeError(
-                "halo_catalog argument must be a HaloCatalog object, "
-                + "a dataset, or a data container."
-            )
-
-        self.width = width
-        self.radius_field = radius_field
-        self.center_field_prefix = center_field_prefix
-        self.annotate_field = annotate_field
-        if circle_kwargs is not None:
-            circle_args = circle_kwargs
-            warnings.warn(
-                "The circle_kwargs keyword is deprecated.  Please "
-                "use the circle_args keyword instead."
-            )
-        if font_kwargs is not None:
-            text_args = font_kwargs
-            warnings.warn(
-                "The font_kwargs keyword is deprecated.  Please use "
-                "the text_args keyword instead."
-            )
-        if circle_args is None:
-            circle_args = def_circle_args
-        self.circle_args = circle_args
-        if text_args is None:
-            text_args = def_text_args
-        self.text_args = text_args
-        self.factor = factor
-
-    def __call__(self, plot):
-        from matplotlib.patches import Circle
-
-        data = plot.data
-        x0, x1, y0, y1 = self._physical_bounds(plot)
-        xx0, xx1, yy0, yy1 = self._plot_bounds(plot)
-
-        halo_data = self.halo_data
-        axis_names = plot.data.ds.coordinates.axis_name
-        xax = plot.data.ds.coordinates.x_axis[data.axis]
-        yax = plot.data.ds.coordinates.y_axis[data.axis]
-        field_x = f"{self.center_field_prefix}_{axis_names[xax]}"
-        field_y = f"{self.center_field_prefix}_{axis_names[yax]}"
-        field_z = f"{self.center_field_prefix}_{axis_names[data.axis]}"
-
-        # Set up scales for pixel size and original data
-        pixel_scale = self._pixel_scale(plot)[0]
-        units = plot.xlim[0].units
-
-        # Convert halo positions to code units of the plotted data
-        # and then to units of the plotted window
-        px = halo_data[("all", field_x)][:].in_units(units)
-        py = halo_data[("all", field_y)][:].in_units(units)
-
-        xplotcenter = (plot.xlim[0] + plot.xlim[1]) / 2
-        yplotcenter = (plot.ylim[0] + plot.ylim[1]) / 2
-
-        xdomaincenter = plot.ds.domain_center[xax]
-        ydomaincenter = plot.ds.domain_center[yax]
-
-        xoffset = xplotcenter - xdomaincenter
-        yoffset = yplotcenter - ydomaincenter
-
-        xdw = plot.ds.domain_width[xax].to(units)
-        ydw = plot.ds.domain_width[yax].to(units)
-
-        modpx = np.mod(px - xoffset, xdw) + xoffset
-        modpy = np.mod(py - yoffset, ydw) + yoffset
-
-        px[modpx != px] = modpx[modpx != px]
-        py[modpy != py] = modpy[modpy != py]
-
-        px, py = self._convert_to_plot(plot, [px, py])
-
-        # Convert halo radii to a radius in pixels
-        radius = halo_data[("all", self.radius_field)][:].in_units(units)
-        radius = np.array(radius * pixel_scale * self.factor)
-
-        if self.width:
-            pz = halo_data[("all", field_z)][:].in_units("code_length")
-            c = data.center[data.axis]
-
-            # I should catch an error here if width isn't in this form
-            # but I dont really want to reimplement get_sanitized_width...
-            width = data.ds.arr(self.width[0], self.width[1]).in_units("code_length")
-
-            indices = np.where((pz > c - 0.5 * width) & (pz < c + 0.5 * width))
-
-            px = px[indices]
-            py = py[indices]
-            radius = radius[indices]
-
-        for x, y, r in zip(px, py, radius):
-            plot._axes.add_artist(Circle(xy=(x, y), radius=r, **self.circle_args))
-
-        plot._axes.set_xlim(xx0, xx1)
-        plot._axes.set_ylim(yy0, yy1)
-
-        if self.annotate_field:
-            annotate_dat = halo_data[("all", self.annotate_field)]
-            texts = [f"{float(dat):g}" for dat in annotate_dat]
-            labels = []
-            for pos_x, pos_y, t in zip(px, py, texts):
-                labels.append(plot._axes.text(pos_x, pos_y, t, **self.text_args))
-
-            # Set the font properties of text from this callback to be
-            # consistent with other text labels in this figure
-            self._set_font_properties(plot, labels, **self.text_args)
-
-
 class ParticleCallback(PlotCallback):
     """
     Adds particle positions, based on a thick slab along *axis* with a