diff --git a/grizli/fitting.py b/grizli/fitting.py
index a681c98b..0394d640 100644
--- a/grizli/fitting.py
+++ b/grizli/fitting.py
@@ -2513,7 +2513,7 @@ def xfit_redshift(self, prior=None,
             
             
         """
-        from numpy import polyfit, polyval
+        from numpy.polynomial.polynomial import polyfit, polyval
         from scipy.stats import t as student_t
         from scipy.special import huber
         import peakutils
@@ -3480,7 +3480,7 @@ def objfun_scale(pscale, AxT, data, self, retval):
         spectra
         """
         import scipy.optimize
-        from numpy import polyval
+        from numpy.polynomial.polynomial import polyval
 
         scale = self.compute_scale_array(pscale, self.wavef[self.fit_mask])
         scale[-self.Nphot:] = 1.
diff --git a/grizli/multifit.py b/grizli/multifit.py
index 92828385..83f84706 100644
--- a/grizli/multifit.py
+++ b/grizli/multifit.py
@@ -2751,7 +2751,7 @@ def fit_redshift(self, prior=None, poly_order=1, fwhm=1200,
                      fsps_templates=False):
         """TBD
         """
-        from scipy import polyfit, polyval
+        from numpy.polynomial.polynomial import polyfit, polyval
 
         if zr is None:
             zr = [0.65, 1.6]
diff --git a/grizli/prep.py b/grizli/prep.py
index 168c92c4..91f1fb25 100644
--- a/grizli/prep.py
+++ b/grizli/prep.py
@@ -5144,7 +5144,7 @@ def tweak_align(direct_group={}, grism_group={}, max_dist=1., n_min=10, key=' ',
                                  'prep.tweak_align')
 
     from drizzlepac.astrodrizzle import AstroDrizzle
-    from scipy import polyfit, polyval
+    from numpy.polynomial.polynomial import polyfit, polyval
 
     if len(direct_group['files']) < 2:
         logstr = '# ! {0}: Only one direct image found, can\'t compute shifts'
diff --git a/grizli/stack.py b/grizli/stack.py
index fc4d79ca..7dd84eaf 100755
--- a/grizli/stack.py
+++ b/grizli/stack.py
@@ -561,7 +561,7 @@ def scale_AxT(p, Ax, spec_wave, Nphot, Next):
         """
         Scale spectrum templates by polynomial function
         """
-        from scipy import polyval
+        from numpy.polynomial.polynomial import polyval
 
         scale = np.ones(Ax.shape[1])
         scale[:-Nphot] = polyval(p[::-1]/10., (spec_wave-1.e4)/1000.)
@@ -578,7 +578,7 @@ def objective_scale(p, Ax, data, spec_wave, fit_mask, sivarf, Nphot, Next, retur
         spectra
         """
         import scipy.optimize
-        from scipy import polyval
+        from numpy.polynomial.polynomial import polyval
 
         scale = np.ones(Ax.shape[1])
         scale[:-Nphot] = polyval(p[::-1]/10., (spec_wave-1.e4)/1000.)