diff --git a/sedpy/observate.py b/sedpy/observate.py
index 98b36fe..893e1aa 100644
--- a/sedpy/observate.py
+++ b/sedpy/observate.py
@@ -7,6 +7,11 @@
 """
 
 import numpy as np
+try:
+    from numpy import trapezoid
+except(ImportError):
+    from numpy import trapz as trapezoid
+
 import os
 try:
     from pkg_resources import resource_filename, resource_listdir
@@ -238,13 +243,13 @@ def get_properties(self):
         of many of these quantities.
         """
         # Calculate some useful integrals
-        i0 = np.trapz(self.transmission * np.log(self.wavelength),
+        i0 = trapezoid(self.transmission * np.log(self.wavelength),
                       np.log(self.wavelength))
-        i1 = np.trapz(self.transmission,
+        i1 = trapezoid(self.transmission,
                       np.log(self.wavelength))
-        i2 = np.trapz(self.transmission * self.wavelength,
+        i2 = trapezoid(self.transmission * self.wavelength,
                       self.wavelength)
-        i3 = np.trapz(self.transmission,
+        i3 = trapezoid(self.transmission,
                       self.wavelength)
 
         tmax = self.transmission.max()
@@ -262,14 +267,14 @@ def get_properties(self):
         rind = np.argmin(self.wavelength[sel])
         self.red_edge = self.wavelength[sel][rind]
 
-        i4 = np.trapz(self.transmission *
+        i4 = trapezoid(self.transmission *
                       (np.log(self.wavelength / self.wave_effective))**2.0,
                       np.log(self.wavelength))
         self.gauss_width = (i4 / i1)**(0.5)
         self.effective_width = (2.0 * np.sqrt(2. * np.log(2.)) *
                                 self.gauss_width *
                                 self.wave_effective)
-        # self.norm  = np.trapz(transmission,wavelength)
+        # self.norm  = trapezoid(transmission,wavelength)
 
         # Get zero points and AB to Vega conversion
         self.ab_zero_counts = self.obj_counts(self.wavelength,
@@ -340,7 +345,7 @@ def obj_counts_hires(self, sourcewave, sourceflux, sourceflux_unc=0):
             # assert ~edge, "Source spectrum does not span filter."
             ind = slice(max(positive.min() - 1, 0),
                         min(positive.max() + 2, len(sourcewave)))
-            counts = np.trapz(sourcewave[ind] * newtrans[ind] *
+            counts = trapezoid(sourcewave[ind] * newtrans[ind] *
                               sourceflux[..., ind],
                               sourcewave[ind], axis=-1)
             return np.squeeze(counts)
@@ -376,7 +381,7 @@ def obj_counts_lores(self, sourcewave, sourceflux, sourceflux_unc=0):
 
         # Integrate lambda*f_lambda*R
         if True in (newflux > 0.):
-            counts = np.trapz(self.wavelength * self.transmission * newflux,
+            counts = trapezoid(self.wavelength * self.transmission * newflux,
                               self.wavelength)
             return np.squeeze(counts)
         else:
@@ -564,7 +569,7 @@ def _build_super_trans(self):
         """
         # TODO: investigate using arrays that cover different wavelength
         # intervals (but same number of elements) for each filter. Then when
-        # integrating (with dot or vectorized trapz) the source spectrum must be
+        # integrating (with dot or vectorized trapezoid) the source spectrum must be
         # offset for each filter. This would remove the vast majority of zeros
 
         ab_counts = np.array([f.ab_zero_counts for f in self.filters])
@@ -664,7 +669,7 @@ def get_sed_maggies_trapz(self, sourceflux, sourcewave=None):
             flux = sourceflux
         y = self.trans_dw * flux[self.finds]
         x = self.lam[self.finds]
-        maggies = np.trapz(y, x, axis=-1)
+        maggies = trapezoid(y, x, axis=-1)
         return maggies
 
     def get_sed_maggies_numba(self, sourceflux, sourcewave=None):
@@ -854,7 +859,7 @@ def Lbol(wave, spec, wave_min=90, wave_max=1e6):
        of length (...nsource)
     """
     inds = np.where(np.logical_and(wave < wave_max, wave >= wave_min))
-    return np.trapz(spec[..., inds[0]], wave[inds])
+    return trapezoid(spec[..., inds[0]], wave[inds])
 
 
 def air2vac(air):
diff --git a/sedpy/smoothing.py b/sedpy/smoothing.py
index 4a54d85..62c4fd8 100644
--- a/sedpy/smoothing.py
+++ b/sedpy/smoothing.py
@@ -6,6 +6,11 @@
 #  * more testing of the smooth_lsf* methods.
 
 import numpy as np
+try:
+    from numpy import trapezoid
+except(ImportError):
+    from numpy import trapz as trapezoid
+
 from numpy.fft import fft, ifft, fftfreq, rfftfreq
 
 __all__ = ["smoothspec", "smooth_wave", "smooth_vel", "smooth_lsf",
@@ -229,7 +234,7 @@ def smooth_vel(wave, spec, outwave, sigma, nsigma=10, inres=0, **extras):
         else:
             _spec = spec
         f = np.exp(-0.5 * x**2)
-        flux[i] = np.trapz(f * _spec, x) / np.trapz(f, x)
+        flux[i] = trapezoid(f * _spec, x) / trapezoid(f, x)
     return flux
 
 
@@ -352,7 +357,7 @@ def smooth_wave(wave, spec, outwave, sigma, nsigma=10, inres=0, in_vel=False,
         else:
             _spec = spec
         f = np.exp(-0.5 * x**2)
-        flux[i] = np.trapz(f * _spec, x) / np.trapz(f, x)
+        flux[i] = trapezoid(f * _spec, x) / trapezoid(f, x)
     return flux
 
 
@@ -463,8 +468,8 @@ def smooth_lsf(wave, spec, outwave, sigma=None, lsf=None, return_kernel=False,
     # should this be axis=0 or axis=1?
     kernel = kernel / kernel.sum(axis=1)[:, None]
     newspec = np.dot(kernel, spec)
-    # kernel /= np.trapz(kernel, wave, axis=1)[:, None]
-    # newspec = np.trapz(kernel * spec[None, :], wave, axis=1)
+    # kernel /= trapezoid(kernel, wave, axis=1)[:, None]
+    # newspec = trapezoid(kernel * spec[None, :], wave, axis=1)
     if return_kernel:
         return newspec, kernel
     return newspec