diff --git a/docs/conf.py b/docs/conf.py
index 55f64c9..2eee0cf 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -374,7 +374,6 @@
     'contextlib.closing',
     'h5py',
     'hashlib',
-    'healpy',
     'numpy',
     'PIL',
     'PIL.Image',
diff --git a/docs/installation.rst b/docs/installation.rst
index 8cd240a..439b89b 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -90,8 +90,8 @@ satisfied:
 * :code:`numpy`
 * :code:`scipy`
 * :code:`astropy`
+* :code:`astropy_healpix`
 * :code:`h5py`
-* :code:`healpy`
 * :code:`requests`
 * :code:`six`
 * :code:`progressbar2`
diff --git a/dustmaps/bayestar.py b/dustmaps/bayestar.py
index aabcdcf..1634ded 100644
--- a/dustmaps/bayestar.py
+++ b/dustmaps/bayestar.py
@@ -30,7 +30,7 @@
 import astropy.coordinates as coordinates
 import astropy.units as units
 import h5py
-import healpy as hp
+from astropy_healpix.healpy import ang2pix
 
 from .std_paths import *
 from .map_base import DustMap, WebDustMap, ensure_flat_galactic
@@ -62,14 +62,14 @@ def lb2pix(nside, l, b, nest=True):
         if (b < -90.) or (b > 90.):
             return -1
 
-        pix_idx = hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)
+        pix_idx = ang2pix(nside, theta, phi, nest=nest)
 
         return pix_idx
 
     idx = (b >= -90.) & (b <= 90.)
 
     pix_idx = np.empty(l.shape, dtype='i8')
-    pix_idx[idx] = hp.pixelfunc.ang2pix(nside, theta[idx], phi[idx], nest=nest)
+    pix_idx[idx] = ang2pix(nside, theta[idx], phi[idx], nest=nest)
     pix_idx[~idx] = -1
 
     return pix_idx
diff --git a/dustmaps/csfd.py b/dustmaps/csfd.py
index 682a742..90e6915 100644
--- a/dustmaps/csfd.py
+++ b/dustmaps/csfd.py
@@ -24,7 +24,6 @@
 
 import os
 import numpy as np
-import healpy as hp
 import astropy.io.fits as fits
 import astropy.units as units
 
diff --git a/dustmaps/edenhofer2023.py b/dustmaps/edenhofer2023.py
index 6066704..5933d7e 100644
--- a/dustmaps/edenhofer2023.py
+++ b/dustmaps/edenhofer2023.py
@@ -169,7 +169,7 @@ def _get_sphere(filepath):
 def _interp_hpxr2lbd(data, radii, nside, nest, lon, lat, dist):
     """Interpolate a 3D map of HEALPix times radii to arbitrary longitude,
     latitude, distance positions."""
-    from healpy.pixelfunc import get_interp_weights
+    from astropy_healpix.healpy import get_interp_weights
 
     assert lon.shape == lat.shape == dist.shape
     final_shape = data.shape[:-2] + lon.shape
diff --git a/dustmaps/gaia_tge.py b/dustmaps/gaia_tge.py
index 65c5abd..5a0805f 100644
--- a/dustmaps/gaia_tge.py
+++ b/dustmaps/gaia_tge.py
@@ -24,7 +24,6 @@
 
 import os
 import numpy as np
-import healpy as hp
 from astropy.table import Table
 import astropy.units as units
 
diff --git a/dustmaps/healpix_map.py b/dustmaps/healpix_map.py
index 6b176b8..3b4f8fe 100644
--- a/dustmaps/healpix_map.py
+++ b/dustmaps/healpix_map.py
@@ -24,7 +24,7 @@
 import six
 
 import numpy as np
-import healpy as hp
+from astropy_healpix.healpy import npix2nside
 import astropy.io.fits as fits
 
 from .map_base import DustMap, coord2healpix
@@ -46,7 +46,7 @@ def __init__(self, pix_val, nest, coord_frame, flags=None):
             coord_frame (str): The coordinate system that the HEALPix map is in.
                 Should be one of the frames supported by `astropy.coordinates`.
         """
-        self._nside = hp.pixelfunc.npix2nside(len(pix_val))
+        self._nside = npix2nside(len(pix_val))
         self._pix_val = pix_val
         self._nest = nest
         self._frame = coord_frame
diff --git a/dustmaps/lenz2017.py b/dustmaps/lenz2017.py
index a8e68d3..8322f64 100644
--- a/dustmaps/lenz2017.py
+++ b/dustmaps/lenz2017.py
@@ -26,7 +26,6 @@
 
 import os
 import numpy as np
-import healpy as hp
 import astropy.io.fits as fits
 import astropy.units as units
 
diff --git a/dustmaps/map_base.py b/dustmaps/map_base.py
index b08fe28..de39014 100755
--- a/dustmaps/map_base.py
+++ b/dustmaps/map_base.py
@@ -23,7 +23,7 @@
 from __future__ import print_function, division
 
 import numpy as np
-import healpy as hp
+from astropy_healpix.healpy import ang2pix, vec2pix
 import astropy.coordinates as coordinates
 import astropy.units as units
 
@@ -65,17 +65,17 @@ def coord2healpix(coords, frame, nside, nest=True):
     if hasattr(c, 'ra'):
         phi = c.ra.rad
         theta = 0.5*np.pi - c.dec.rad
-        return hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)
+        return ang2pix(nside, theta, phi, nest=nest)
     elif hasattr(c, 'l'):
         phi = c.l.rad
         theta = 0.5*np.pi - c.b.rad
-        return hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)
+        return ang2pix(nside, theta, phi, nest=nest)
     elif hasattr(c, 'x'):
         x,y,z = [v.to('kpc').value for v in (c.x,c.y,c.z)]
-        return hp.pixelfunc.vec2pix(nside, x, y, z, nest=nest)
+        return vec2pix(nside, x, y, z, nest=nest)
     elif hasattr(c, 'w'):
         x,y,z = [v.to('kpc').value for v in (c.u,c.v,c.w)]
-        return hp.pixelfunc.vec2pix(nside, x, y, z, nest=nest)
+        return vec2pix(nside, x, y, z, nest=nest)
     else:
         raise dustexceptions.CoordFrameError(
             'No method to transform from coordinate frame "{}" to HEALPix.'.format(
diff --git a/dustmaps/planck.py b/dustmaps/planck.py
index 8faf439..544ace2 100644
--- a/dustmaps/planck.py
+++ b/dustmaps/planck.py
@@ -24,7 +24,6 @@
 
 import os
 import numpy as np
-import healpy as hp
 import astropy.io.fits as fits
 import astropy.units as units
 
diff --git a/setup.py b/setup.py
index 535b411..3561917 100644
--- a/setup.py
+++ b/setup.py
@@ -183,7 +183,7 @@ def readme():
         'scipy',
         'astropy',
         'h5py',
-        'healpy',
+        'astropy_healpix',
         'requests',
         'progressbar2',
         'six'