diff --git a/examples/galaxies/plot_schechter.py b/examples/galaxies/plot_schechter.py
index 11b060b4..51501279 100644
--- a/examples/galaxies/plot_schechter.py
+++ b/examples/galaxies/plot_schechter.py
@@ -47,6 +47,7 @@
 from astropy.units import Quantity
 from matplotlib import pyplot as plt
 import numpy as np
+import scipy.integrate
 from skypy.galaxies import schechter_lf
 
 z_range = np.linspace(0.2, 1.0, 100)
@@ -81,7 +82,7 @@
     # SkyPy simulated galaxies
     z_mask = np.logical_and(redshift >= z_min, redshift < z_max)
     dV_dz = (cosmology.differential_comoving_volume(z) * sky_area).to_value('Mpc3')
-    dV = np.trapz(dV_dz, z)
+    dV = scipy.integrate.trapezoid(dV_dz, z)
     dM = (np.max(bins)-np.min(bins)) / (np.size(bins)-1)
     phi_skypy = np.histogram(magnitude[z_mask], bins=bins)[0] / dV / dM
 
diff --git a/skypy/galaxies/redshift.py b/skypy/galaxies/redshift.py
index f048a357..3173c63d 100644
--- a/skypy/galaxies/redshift.py
+++ b/skypy/galaxies/redshift.py
@@ -249,7 +249,7 @@ def redshifts_from_comoving_density(redshift, density, sky_area, cosmology, nois
     dN_dz *= density
 
     # integrate density to get expected number of galaxies
-    N = np.trapz(dN_dz, redshift)
+    N = scipy.integrate.trapezoid(dN_dz, redshift)
 
     # Poisson sample galaxy number if requested
     if noise:
diff --git a/skypy/galaxies/tests/test_redshift.py b/skypy/galaxies/tests/test_redshift.py
index d41e92ac..c974bab3 100644
--- a/skypy/galaxies/tests/test_redshift.py
+++ b/skypy/galaxies/tests/test_redshift.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+import scipy.integrate
 from scipy.stats import kstest
 
 
@@ -38,7 +39,7 @@ def test_schechter_lf_redshift():
     density *= (sky_area * cosmo.differential_comoving_volume(z)).to_value('Mpc3')
 
     # integrate total number
-    n_gal = np.trapz(density, z, axis=-1)
+    n_gal = scipy.integrate.trapezoid(density, z, axis=-1)
 
     # make sure noise-free sample has right size
     assert np.isclose(len(z_gal), n_gal, atol=1.0)
@@ -85,7 +86,7 @@ def test_schechter_smf_redshift():
     density *= (sky_area * cosmo.differential_comoving_volume(z)).to_value('Mpc3')
 
     # integrate total number
-    n_gal = np.trapz(density, z, axis=-1)
+    n_gal = scipy.integrate.trapezoid(density, z, axis=-1)
 
     # make sure noise-free sample has right size
     assert np.isclose(len(z_gal), n_gal, atol=1.0)