diff --git a/CHANGES.rst b/CHANGES.rst
index 7500ad652f..b178c76ab9 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -43,7 +43,7 @@ associations
 - Excluded nearby background candidates from NIRSpec fixed slit associations
   for S1600A1 with 5 point dithers, to reduce overlap between background nods
   and science exposure. [#8744]
-  
+
 - Added association rule for level 3 image mosaic candidates. [#8798]
 
 badpix_selfcal
@@ -274,7 +274,7 @@ pipeline
   optional `on_disk` parameter to govern whether models in the library should be stored
   in memory or on disk. [#8683]
 
-- Updated ``calwebb_spec2`` to run ``nsclean`` on NIRSpec imprint and background 
+- Updated ``calwebb_spec2`` to run ``nsclean`` on NIRSpec imprint and background
   association members. [#8786, #8809]
 
 - Updated `calwebb_spec3` to not save the `pixel_replacement` output by default.[#8765]
diff --git a/changes/8866.resample.rst b/changes/8866.resample.rst
new file mode 100644
index 0000000000..f31024d0ee
--- /dev/null
+++ b/changes/8866.resample.rst
@@ -0,0 +1 @@
+Updated resample code to support the new ``drizzle`` API, see https://github.com/spacetelescope/drizzle/pull/134 for more details.
diff --git a/jwst/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py
index 63f0a84ad7..22450fc87e 100644
--- a/jwst/resample/gwcs_drizzle.py
+++ b/jwst/resample/gwcs_drizzle.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from drizzle import cdrizzle, utils
+from drizzle import cdrizzle
 from . import resample_utils
 
 import logging
@@ -380,13 +380,8 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon,
 
     # Compute the mapping between the input and output pixel coordinates
     # for use in drizzle.cdrizzle.tdriz
-
-    # pixmap = utils.calc_pixmap(input_wcs, output_wcs)
     pixmap = resample_utils.calc_gwcs_pixmap(input_wcs, output_wcs, insci.shape)
 
-    # inwht[np.isnan(pixmap[:,:,0])] = 0.
-
-    print(f"Pixmap shape: {pixmap.shape}")
     log.debug(f"Pixmap shape: {pixmap[:,:,0].shape}")
     log.debug(f"Input Sci shape: {insci.shape}")
     log.debug(f"Output Sci shape: {outsci.shape}")
diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py
index a0d1b53d1f..c7501cc8ba 100644
--- a/jwst/resample/resample.py
+++ b/jwst/resample/resample.py
@@ -192,7 +192,6 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True,
             output_pix_area * np.rad2deg(3600)**2
         )
 
-
     def do_drizzle(self, input_models):
         """Pick the correct drizzling mode based on self.single
         """
@@ -357,7 +356,7 @@ def resample_many_to_many(self, input_models):
             # build ModelLibrary as an association from the output files
             # this saves memory if there are multiple groups
             asn = asn_from_list(output_models, product_name='outlier_i2d')
-            asn_dict = json.loads(asn.dump()[1]) # serializes the asn and converts to dict
+            asn_dict = json.loads(asn.dump()[1])  # serializes the asn and converts to dict
             return ModelLibrary(asn_dict, on_disk=True)
         # otherwise just build it as a list of in-memory models
         return ModelLibrary(output_models, on_disk=False)
@@ -397,9 +396,11 @@ def resample_many_to_one(self, input_models):
                 log.debug(f'Using intensity scale iscale={iscale}')
                 img.meta.iscale = iscale
 
-                inwht = resample_utils.build_driz_weight(img,
-                                                        weight_type=self.weight_type,
-                                                        good_bits=self.good_bits)
+                inwht = resample_utils.build_driz_weight(
+                    img,
+                    weight_type=self.weight_type,
+                    good_bits=self.good_bits,
+                )
                 # apply sky subtraction
                 blevel = img.meta.background.level
                 if not img.meta.background.subtracted and blevel is not None:
@@ -446,14 +447,13 @@ def resample_many_to_one(self, input_models):
 
         return ModelLibrary([output_model,], on_disk=False)
 
-
     def resample_variance_arrays(self, output_model, input_models):
         """Resample variance arrays from input_models to the output_model.
 
         Variance images from each input model are resampled individually and
-        added to a weighted sum. If weight_type is 'ivm', the inverse of the
+        added to a weighted sum. If ``weight_type`` is 'ivm', the inverse of the
         resampled read noise variance is used as the weight for all the variance
-        components. If weight_type is 'exptime', the exposure time is used.
+        components. If ``weight_type`` is 'exptime', the exposure time is used.
 
         The output_model is modified in place.
         """
@@ -493,8 +493,10 @@ def resample_variance_arrays(self, output_model, input_models):
                 if rn_var is not None:
                     mask = (rn_var >= 0) & np.isfinite(rn_var) & (weight > 0)
                     weighted_rn_var[mask] = np.nansum(
-                        [weighted_rn_var[mask],
-                        rn_var[mask] * weight[mask] * weight[mask]],
+                        [
+                            weighted_rn_var[mask],
+                            rn_var[mask] * weight[mask] * weight[mask]
+                        ],
                         axis=0
                     )
                     total_weight_rn_var[mask] += weight[mask]
@@ -506,8 +508,10 @@ def resample_variance_arrays(self, output_model, input_models):
                 if pn_var is not None:
                     mask = (pn_var >= 0) & np.isfinite(pn_var) & (weight > 0)
                     weighted_pn_var[mask] = np.nansum(
-                        [weighted_pn_var[mask],
-                        pn_var[mask] * weight[mask] * weight[mask]],
+                        [
+                            weighted_pn_var[mask],
+                            pn_var[mask] * weight[mask] * weight[mask]
+                        ],
                         axis=0
                     )
                     total_weight_pn_var[mask] += weight[mask]
@@ -517,8 +521,10 @@ def resample_variance_arrays(self, output_model, input_models):
                 if flat_var is not None:
                     mask = (flat_var >= 0) & np.isfinite(flat_var) & (weight > 0)
                     weighted_flat_var[mask] = np.nansum(
-                        [weighted_flat_var[mask],
-                        flat_var[mask] * weight[mask] * weight[mask]],
+                        [
+                            weighted_flat_var[mask],
+                            flat_var[mask] * weight[mask] * weight[mask]
+                        ],
                         axis=0
                     )
                     total_weight_flat_var[mask] += weight[mask]
@@ -549,7 +555,6 @@ def resample_variance_arrays(self, output_model, input_models):
             del weighted_rn_var, weighted_pn_var, weighted_flat_var
             del total_weight_rn_var, total_weight_pn_var, total_weight_flat_var
 
-
     def _resample_one_variance_array(self, name, input_model, output_model):
         """Resample one variance image from an input model.
 
diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py
index 0ff50902de..3775cb5377 100644
--- a/jwst/resample/resample_utils.py
+++ b/jwst/resample/resample_utils.py
@@ -4,7 +4,6 @@
 
 import numpy as np
 from astropy import units as u
-import gwcs
 from drizzle.utils import calc_pixmap
 
 from stdatamodels.dqflags import interpret_bit_flags
@@ -125,7 +124,7 @@ def shape_from_bounding_box(bounding_box):
 def calc_gwcs_pixmap(in_wcs, out_wcs, shape=None):
     """ Return a pixel grid map from input frame to output frame.
     """
-    if shape and not np.array_equiv(shape, in_wcs.array_shape):
+    if shape is not None and not np.array_equiv(shape, in_wcs.array_shape):
         in_wcs = deepcopy(in_wcs)
         in_wcs.array_shape = shape