diff --git a/drizzlepac/haputils/astrometric_utils.py b/drizzlepac/haputils/astrometric_utils.py
index c8bb530ec..6b98b39ad 100644
--- a/drizzlepac/haputils/astrometric_utils.py
+++ b/drizzlepac/haputils/astrometric_utils.py
@@ -46,16 +46,16 @@
 from astropy.utils.decorators import deprecated
 
 import photutils  # needed to check version
-if Version(photutils.__version__) < Version('1.1.0'):
+from photutils.segmentation import (detect_sources, SourceCatalog,
+                                    deblend_sources, SegmentationImage)
+if Version(photutils.__version__) < Version('1.5.0'):
     OLD_PHOTUTILS = True
-    from photutils.segmentation import (detect_sources,
-                                        deblend_sources, make_source_mask)
-    from photutils.segmentation import source_properties as SourceCatalog
+    from photutils.segmentation import make_source_mask
 else:
     OLD_PHOTUTILS = False
-    from photutils.segmentation import (detect_sources, SourceCatalog,
-                                        deblend_sources, make_source_mask,
-                                        SegmentationImage)
+    from photutils.segmentation import detect_threshold
+    from photutils.utils import circular_footprint
+    from astropy.stats import SigmaClip
 
 from photutils.detection import DAOStarFinder, find_peaks
 
@@ -694,8 +694,18 @@ def compute_2d_background(imgarr, box_size, win_size,
     # the background to be used in source identification
     if bkg is None:
         log.info("Background2D failure detected. Using alternative background calculation instead....")
-        mask = make_source_mask(imgarr, nsigma=2, npixels=5, dilate_size=11)
-        sigcl_mean, sigcl_median, sigcl_std = sigma_clipped_stats(imgarr, sigma=3.0, mask=mask, maxiters=9)
+
+        if OLD_PHOTUTILS:
+            mask = make_source_mask(imgarr, nsigma=2, npixels=5, dilate_size=11)
+            sigcl_mean, sigcl_median, sigcl_std = sigma_clipped_stats(imgarr, sigma=3.0, mask=mask, maxiters=9)
+        else:
+            sigma_clip = SigmaClip(sigma=3.0, maxiters=9)
+            threshold = detect_threshold(imgarr, nsigma=2.0, sigma_clip=sigma_clip)
+            segment_img = detect_sources(imgarr, threshold, npixels=5)
+            footprint = circular_footprint(radius=11)
+            mask = segment_img.make_source_mask(footprint=footprint)
+            sigcl_mean, sigcl_median, sigcl_std = sigma_clipped_stats(imgarr, sigma=3.0, mask=mask)
+
         bkg_median = max(0.0, sigcl_median)
         bkg_rms_median = sigcl_std
         # create background frame shaped like imgarr populated with sigma-clipped median value
@@ -1033,12 +1043,8 @@ def extract_sources(img, dqmask=None, fwhm=3.0, kernel=None, photmode=None,
 
         segm.remove_labels(bad_srcs)
 
-    if OLD_PHOTUTILS:
-        flux_colname = 'source_sum'
-        ferr_colname = 'source_sum_err'
-    else:
-        flux_colname = 'segment_flux'
-        ferr_colname = 'segment_fluxerr'
+    flux_colname = 'segment_flux'
+    ferr_colname = 'segment_fluxerr'
 
     # convert segm to mask for daofind
     if centering_mode == 'starfind':
@@ -1103,23 +1109,15 @@ def extract_sources(img, dqmask=None, fwhm=3.0, kernel=None, photmode=None,
                 if (detection_img > detection_img.max() * 0.9).sum() > 3:
 
                     # Revert to segmentation photometry for sat. source posns
-                    if OLD_PHOTUTILS:
-                        segment_properties = SourceCatalog(detection_img, segment.data)
-                    else:
-                        segimg = SegmentationImage(segment.data)
-                        segment_properties = SourceCatalog(detection_img, segimg)
+                    segimg = SegmentationImage(segment.data)
+                    segment_properties = SourceCatalog(detection_img, segimg)
 
                     sat_table = segment_properties.to_table()
                     seg_table['flux'][max_row] = sat_table[flux_colname][0]
                     seg_table['peak'][max_row] = sat_table['max_value'][0]
-                    if OLD_PHOTUTILS:
-                        xcentroid = sat_table['xcentroid'][0].value
-                        ycentroid = sat_table['ycentroid'][0].value
-                        sky = sat_table['background_mean'][0].value
-                    else:
-                        xcentroid = sat_table['xcentroid'][0]
-                        ycentroid = sat_table['ycentroid'][0]
-                        sky = sat_table['local_background'][0]
+                    xcentroid = sat_table['xcentroid'][0]
+                    ycentroid = sat_table['ycentroid'][0]
+                    sky = sat_table['local_background'][0]
                     seg_table['xcentroid'][max_row] = xcentroid
                     seg_table['ycentroid'][max_row] = ycentroid
                     seg_table['npix'][max_row] = sat_table['area'][0].value
diff --git a/pyproject.toml b/pyproject.toml
index e7343b81a..c474c5d70 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,7 @@ dependencies = [
     'pandas',
     'spherical_geometry>=1.2.22',
     'astroquery>=0.4',
-    'photutils>1.2.0',
+    'photutils>1.5.0',
     'lxml',
     'PyPDF2',
     'scikit-image>=0.14.2',
@@ -158,6 +158,7 @@ extend-exclude = [
     '.tox',
 ]
 ignore = [
+    'E401', # Multiple imports not allowed on single line
     'E402', # Module level import not at top of file
     'E501', # Line too long
     'E711', # Comparison to `None` should be `cond is None`
diff --git a/tests/hap/test_pipeline.py b/tests/hap/test_pipeline.py
index 272cc0a32..33e4b4d04 100644
--- a/tests/hap/test_pipeline.py
+++ b/tests/hap/test_pipeline.py
@@ -119,7 +119,7 @@ class BaseWFC3Pipeline(BasePipeline):
 class TestSingleton(BaseWFC3Pipeline):
 
     @pytest.mark.parametrize(
-        'dataset_names', ['iaaua1n4q', 'iacs01t4q']
+        'dataset_names', ['iaaua1n4q']
     )
 
     def test_astrometric_singleton(self, dataset_names):