diff --git a/pfb/deconv/clean.py b/pfb/deconv/clean.py
deleted file mode 100644
index 99880157..00000000
--- a/pfb/deconv/clean.py
+++ /dev/null
@@ -1,79 +0,0 @@
-import numpy as np
-from pfb.operators.psf import PSF
-from pfb.opt.hogbom import hogbom
-import pyscilog
-log = pyscilog.get_logger('CLEAN')
-
-
-def resid_func(x, dirty, psfo):
-    residual = dirty - psfo.convolve(x)
-    return residual
-
-
-def clean(psf, model, residual, mask=None, beam=None,
-          nthreads=0, maxit=10, gamma=1.0, peak_factor=0.01, threshold=0.0,
-          hbgamma=0.1, hbpf=0.1, hbmaxit=5000, hbverbose=1):  # Hogbom options
-
-    if len(residual.shape) > 3:
-        raise ValueError("Residual must have shape (nband, nx, ny)")
-
-    nband, nx, ny = residual.shape
-
-    if beam is not None:
-        raise NotImplementedError("Beam not supported in clean minor cycle")
-
-    if mask is None:
-        def mask(x): return x
-    else:
-        try:
-            if mask.ndim == 2:
-                assert mask.shape == (nx, ny)
-                def mask(x): return mask[None] * x
-            elif mask.ndim == 3:
-                assert mask.shape == (1, nx, ny)
-                def mask(x): return mask * x
-            else:
-                raise ValueError
-        except BaseException:
-            raise ValueError("Mask has incorrect shape")
-
-    # PSF operator
-    psfo = PSF(psf, nthreads, imsize=residual.shape)
-
-    # init dirty
-    if model.any():
-        dirty = residual + psfo.convolve(model)
-    else:
-        dirty = residual
-
-    # residual
-    residual_mfs = np.sum(residual, axis=0)
-    rms = np.std(residual_mfs)
-    rmax = np.abs(residual_mfs).max()
-    rmax_orig = rmax
-    threshold = np.maximum(peak_factor * rmax, threshold)
-    # deconvolve
-    for i in range(0, maxit):
-        modelu = hogbom(mask(residual), psf,
-                        gamma=hbgamma, pf=hbpf,
-                        maxit=hbmaxit, verbosity=hbverbose)
-
-        model += gamma * modelu
-
-        residual = resid_func(model, dirty, psfo)
-
-        # check stopping criteria
-        residual_mfs = np.sum(residual, axis=0)
-        rmax = np.abs(residual_mfs).max()
-        rms = np.std(residual_mfs)
-
-        if rmax < threshold * rmax_orig:
-            print("Success, convergence after %i iterations."
-                  "Peak of residual is %f, rms is %f" % (i + 1, rmax, rms),
-                  file=log)
-            break
-        else:
-            print("At iteration %i peak of residual is %f, rms is %f" %
-                  (i + 1, rmax, rms), file=log)
-
-    return model
diff --git a/pfb/deconv/spotless.py b/pfb/deconv/spotless.py
deleted file mode 100644
index 83ae5639..00000000
--- a/pfb/deconv/spotless.py
+++ /dev/null
@@ -1,238 +0,0 @@
-import numpy as np
-import scipy
-from pfb.operators.psf import PSF
-from pfb.operators.dirac import Dirac
-from pfb.opt.primal_dual import primal_dual
-from pfb.opt.pcg import pcg
-from pfb.opt.power_method import power_method
-from pfb.opt.hogbom import hogbom
-from pfb.prox.prox_21m import prox_21m
-from pfb.utils.fits import save_fits
-from skimage.filters import threshold_mean
-import pyscilog
-log = pyscilog.get_logger('SPOTLESS')
-
-def make_noise_map(restored_image, boxsize):
-    # Modified version of Cyril's magic minimum filter
-    # Plundered from the depths of
-    # https://github.com/cyriltasse/DDFacet/blob/master/SkyModel/MakeMask.py
-    box = (boxsize, boxsize)
-    n = boxsize**2.0
-    x = np.linspace(-10, 10, 1000)
-    f = 0.5 * (1.0 + scipy.special.erf(x / np.sqrt(2.0)))
-    F = 1.0 - (1.0 - f)**n
-    ratio = np.abs(np.interp(0.5, F, x))
-    noise = -scipy.ndimage.filters.minimum_filter(restored_image, box) / ratio
-    negative_mask = noise < 0.0
-    noise[negative_mask] = 1.0e-10
-    median_noise = np.median(noise)
-    median_mask = noise < median_noise
-    noise[median_mask] = median_noise
-    return noise
-
-
-def resid_func(x, dirty, hessian, mask, beam, wsum):
-    """
-    Returns the unattenuated residual
-    """
-    residual = dirty - hessian(mask(beam(x)))/wsum
-    residual_mfs = np.sum(residual, axis=0)
-    residual = residual
-    return residual, residual_mfs
-
-
-def spotless(psf, model, residual, mask=None, beam_image=None, hessian=None,
-             wsum=1, adapt_sig21=False, cpsf=None, hdr=None, hdr_mfs=None, outfile=None,
-             nthreads=1, sig_21=1e-3, sigma_frac=100, maxit=10, tol=1e-4,
-             peak_factor=0.01, threshold=0.0, positivity=True, gamma=0.9999,
-             hbgamma=0.1, hbpf=0.1, hbmaxit=5000, hbverbose=1,
-             pdtol=1e-4, pdmaxit=250, pdverbose=1,  # primal dual options
-             cgtol=1e-4, cgminit=15, cgmaxit=150, cgverbose=1,  # pcg options
-             pmtol=1e-4, pmmaxit=50, pmverbose=1):  # power method options
-    """
-    Modified clean algorithm:
-
-    psf      - PSF image i.e. R.H W where W contains the weights.
-               Shape must be >= residual.shape
-    model    - current intrinsic model
-    residual - apparent residual image i.e. R.H W (V - R A x)
-
-    Note that peak finding happens in apparent residual because that
-    is where it is easiest to accommodate convolution by the PSF.
-    However, the beam and the mask have to be applied to the residual
-    before we solve for the pre-conditioned updates.
-
-    """
-
-    if len(residual.shape) > 3:
-        raise ValueError("Residual must have shape (nband, nx, ny)")
-
-    nband, nx, ny = residual.shape
-
-    if beam_image is None:
-        def beam(x): return x
-        def beaminv(x): return x
-    else:
-        try:
-            assert beam.shape == (nband, nx, ny)
-            def beam(x): return beam_image * x
-            def beaminv(x): return np.where(beam_image > 0.01,  x / beam_image, x)
-        except BaseException:
-            raise ValueError("Beam has incorrect shape")
-
-    if mask is None:
-        def mask(x): return x
-    else:
-        try:
-            if mask.ndim == 2:
-                assert mask.shape == (nx, ny)
-                def mask(x): return mask[None] * x
-            elif mask.ndim == 3:
-                assert mask.shape == (1, nx, ny)
-                def mask(x): return mask * x
-            else:
-                raise ValueError
-        except BaseException:
-            raise ValueError("Mask has incorrect shape")
-
-    # PSF operator
-    psfo = PSF(psf, residual.shape, nthreads=nthreads, backward_undersize=1.2)
-
-    # set up point sources
-    phi = Dirac(nband, nx, ny, mask=np.any(model, axis=0))
-    dual = np.zeros((nband, nx, ny), dtype=np.float64)
-
-    # clean beam
-    if cpsf is not None:
-        try:
-            assert cpsf.shape == (1,) + psf.shape[1::]
-        except Exception as e:
-            cpsf = cpsf[None, :, :]
-        cpsfo = PSF(cpsf, residual.shape, nthreads=nthreads)
-
-    residual_mfs = np.sum(residual, axis=0)
-    rmax = np.abs(residual_mfs).max()
-    rms = np.std(residual_mfs)
-
-    #  preconditioning operator
-    varmap = np.ones(model.shape) * (sigma_frac * rmax)
-    def hessb(x):
-        return phi.hdot(mask(beam(psfo.convolveb(mask(beam(phi.dot(x))))))) +\
-            x / varmap
-
-    def hessf(x):
-        return phi.hdot(mask(beam(psfo.convolve(mask(beam(phi.dot(x))))))) +\
-                    x / varmap
-
-    beta, betavec = power_method(hessb, residual.shape, tol=pmtol,
-                                 maxit=pmmaxit, verbosity=pmverbose)
-
-    if hessian is None:
-        hessian = psf.convolve
-        wsum = 1.0
-
-    if model.any():
-        dirty = residual + hessian(mask(beam(model)))/wsum
-    else:
-        dirty = residual
-
-    # deconvolve
-    threshold = np.maximum(peak_factor * rmax, threshold)
-    alpha = sig_21
-    for i in range(0, maxit):
-        # find point source candidates
-        modelu = hogbom(mask(residual), psf, gamma=hbgamma,
-                        pf=hbpf, maxit=hbmaxit, verbosity=hbverbose)
-
-        phi.update_locs(modelu)
-
-        # solve for beta updates
-        x = pcg(hessf,
-                phi.hdot(mask(beam(residual))),
-                phi.hdot(beaminv(modelu)),
-                M=lambda x: x * (sigma_frac * rmax), tol=cgtol,
-                maxit=cgmaxit, minit=cgminit, verbosity=cgverbose)
-
-        modelp = model.copy()
-        model += gamma * x
-
-        weights_21 = np.where(phi.mask,
-                              alpha/(alpha + np.abs(np.mean(modelp, axis=0))),
-                              1e10)  # 1e10 for effective infinity
-        beta, betavec = power_method(hessb, model.shape, b0=betavec,
-                                     tol=pmtol, maxit=pmmaxit,
-                                     verbosity=pmverbose)
-
-        model, dual = primal_dual(hessb, model, modelp, dual, sig_21,
-                                  phi, weights_21, beta, prox_21m,
-                                  tol=pdtol, maxit=pdmaxit, axis=0,
-                                  positivity=positivity, report_freq=50,
-                                  verbosity=pdverbose)
-
-        # update Dirac dictionary (remove zero components)
-        phi.trim_fat(model)
-        residual, residual_mfs = resid_func(model, dirty, hessian, mask, beam, wsum)
-
-        model_mfs = np.mean(model, axis=0)
-
-        # check stopping criteria
-        rmax = np.abs(mask(residual_mfs)).max()
-        rms = np.std(mask(residual_mfs))
-        eps = np.linalg.norm(model - modelp) / np.linalg.norm(model)
-
-        # update variance map (positivity constraint optional)
-        varmap = np.maximum(rmax*sigma_frac, sigma_frac*(rmax + model))
-
-        print("Iter %i: peak residual = %f, rms = %f, eps = %f" % (
-                i+1, rmax, rms, eps), file=log)
-
-
-        # save current iteration
-        if outfile is not None:
-            assert hdr is not None
-            assert hdr_mfs is not None
-
-            save_fits(outfile + str(i + 1) + '_model_mfs.fits',
-                      model_mfs, hdr_mfs)
-
-            save_fits(outfile + str(i + 1) + '_model.fits',
-                      model, hdr)
-
-            save_fits(outfile + str(i + 1) + '_update.fits',
-                      x, hdr)
-
-            save_fits(outfile + str(i + 1) + '_residual_mfs.fits',
-                      residual_mfs, hdr_mfs)
-
-            save_fits(outfile + str(i + 1) + '_residual.fits',
-                      residual*wsum, hdr)
-
-
-        if rmax < threshold or eps < tol:
-            print("Success, convergence after %i iterations", file=log)
-            break
-
-        if adapt_sig21:
-            # sig_21 should be set to the std of the image noise
-            from scipy.stats import skew, kurtosis
-            alpha = rms
-            tmp = residual_mfs
-            z = tmp/alpha
-            k = 0
-            while (np.abs(skew(z.ravel(), nan_policy='omit')) > 0.05 or
-                   np.abs(kurtosis(z.ravel(), fisher=True, nan_policy='omit')) > 0.5) and k < 10:
-                # eliminate outliers
-                tmp = np.where(np.abs(z) < 3, residual_mfs, np.nan)
-                alpha = np.nanstd(tmp)
-                z = tmp/alpha
-                print(alpha, skew(z.ravel(), nan_policy='omit'), kurtosis(z.ravel(), fisher=True, nan_policy='omit'))
-                k += 1
-
-            sig_21 = alpha
-            print("alpha set to %f"%(alpha), file=log)
-
-
-
-
-
-    return model
diff --git a/pfb/operators/gridder.py b/pfb/operators/gridder.py
index 107566d2..864d4a7c 100644
--- a/pfb/operators/gridder.py
+++ b/pfb/operators/gridder.py
@@ -9,12 +9,10 @@
 import numba
 import concurrent.futures as cf
 import xarray as xr
-import dask
 import dask.array as da
 from ducc0.wgridder.experimental import vis2dirty, dirty2vis
-from ducc0.fft import c2r, r2c, c2c
+from ducc0.fft import r2c
 from ducc0.misc import resize_thread_pool
-from africanus.constants import c as lightspeed
 from pfb.utils.weighting import counts_to_weights, _compute_counts
 from pfb.utils.beam import eval_beam
 from pfb.utils.naming import xds_from_list
diff --git a/pfb/operators/hessian.py b/pfb/operators/hessian.py
index fd16afc8..30bc27ee 100644
--- a/pfb/operators/hessian.py
+++ b/pfb/operators/hessian.py
@@ -1,12 +1,9 @@
 from functools import partial
 import numpy as np
 import numexpr as ne
-import dask
-import dask.array as da
 from ducc0.wgridder.experimental import vis2dirty, dirty2vis
 from ducc0.fft import r2c, c2r
 from ducc0.misc import empty_noncritical
-from uuid import uuid4
 from pfb.operators.psf import (psf_convolve_slice,
                                psf_convolve_cube)
 from pfb.opt.pcg import pcg
diff --git a/pfb/operators/psf.py b/pfb/operators/psf.py
index ad8a12db..4cd184af 100644
--- a/pfb/operators/psf.py
+++ b/pfb/operators/psf.py
@@ -1,11 +1,7 @@
 import numpy as np
-import numexpr as ne
 import dask.array as da
 from uuid import uuid4
-from ducc0.fft import r2c, c2r, c2c, good_size
-from ducc0.misc import roll_resize_roll as rrr
-from uuid import uuid4
-import gc
+from ducc0.fft import r2c, c2r
 
 
 def psf_convolve_slice(
diff --git a/pfb/operators/psi.py b/pfb/operators/psi.py
index af1362d9..fc48ba0f 100644
--- a/pfb/operators/psi.py
+++ b/pfb/operators/psi.py
@@ -5,9 +5,7 @@
 from numba import types, typed
 from numba.experimental import jitclass
 import pywt
-from scipy.datasets import ascent
 from pfb.wavelets import coeff_size, signal_size, dwt2d, idwt2d, copyT
-from time import time
 
 
 @numba.njit
diff --git a/pfb/opt/pcg.py b/pfb/opt/pcg.py
index 4e63c896..0e2fc422 100644
--- a/pfb/opt/pcg.py
+++ b/pfb/opt/pcg.py
@@ -5,12 +5,10 @@
 import numexpr as ne
 from functools import partial
 import dask.array as da
-from distributed import wait
 from uuid import uuid4
-from pfb.utils.misc import norm_diff, fitcleanbeam, Gaussian2D, taperf, JIT_OPTIONS
+from pfb.utils.misc import norm_diff, JIT_OPTIONS
 from pfb.utils.naming import xds_from_list
 from ducc0.misc import empty_noncritical
-from ducc0.fft import c2c
 iFs = np.fft.ifftshift
 Fs = np.fft.fftshift
 
diff --git a/pfb/opt/power_method.py b/pfb/opt/power_method.py
index f491fc21..a3acae71 100644
--- a/pfb/opt/power_method.py
+++ b/pfb/opt/power_method.py
@@ -1,9 +1,6 @@
 import numpy as np
-import dask.array as da
-from operator import getitem
-from distributed import wait, get_client, as_completed
+from distributed import get_client
 from scipy.linalg import norm
-from copy import deepcopy
 import pyscilog
 log = pyscilog.get_logger('PM')
 
@@ -68,8 +65,6 @@ def bnormf(bsumsq):
 def betaf(beta_num, beta_den):
     return np.sum(beta_num)/np.sum(beta_den)
 
-
-from time import time
 def power_method_dist(actors,
                       nx,
                       ny,
diff --git a/pfb/opt/primal_dual.py b/pfb/opt/primal_dual.py
index cbd82b3b..f7bbce03 100644
--- a/pfb/opt/primal_dual.py
+++ b/pfb/opt/primal_dual.py
@@ -1,16 +1,11 @@
 import numpy as np
 import numexpr as ne
-import dask.array as da
 from pfb.utils.dist import l1reweight_func
-from operator import getitem
-from ducc0.misc import make_noncritical
 from pfb.utils.misc import norm_diff
 from numba import njit, prange
 from uuid import uuid4
 import pyscilog
-import gc
 from time import time
-# gc.set_debug(gc.DEBUG_LEAK)
 log = pyscilog.get_logger('PD')
 
 
diff --git a/pfb/utils/beam.py b/pfb/utils/beam.py
index 21786542..93a8dc0c 100644
--- a/pfb/utils/beam.py
+++ b/pfb/utils/beam.py
@@ -1,9 +1,6 @@
 import numpy as np
 from scipy.interpolate import RegularGridInterpolator as RGI
-from functools import partial
 from katbeam import JimBeam
-import dask.array as da
-from numba.core.errors import NumbaDeprecationWarning
 from africanus.rime.fast_beam_cubes import beam_cube_dde
 from africanus.rime import parallactic_angles
 
diff --git a/pfb/utils/correlations.py b/pfb/utils/correlations.py
index a29f53f9..b46a9ae4 100644
--- a/pfb/utils/correlations.py
+++ b/pfb/utils/correlations.py
@@ -1,7 +1,5 @@
 import numpy as np
 from numba import njit
-from numba.types import literal
-from dask.graph_manipulation import clone
 import dask.array as da
 from xarray import Dataset
 from operator import getitem
diff --git a/pfb/utils/fits.py b/pfb/utils/fits.py
index 8c1c1982..c1b12856 100644
--- a/pfb/utils/fits.py
+++ b/pfb/utils/fits.py
@@ -1,8 +1,6 @@
 import numpy as np
 from astropy.io import fits
 from astropy.wcs import WCS
-import dask.array as da
-from dask import delayed
 from datetime import datetime
 from casacore.quanta import quantity
 from astropy.time import Time
diff --git a/pfb/utils/misc.py b/pfb/utils/misc.py
index 87fc5548..c196168e 100644
--- a/pfb/utils/misc.py
+++ b/pfb/utils/misc.py
@@ -2,8 +2,7 @@
 from contextlib import nullcontext
 import numpy as np
 import numexpr as ne
-import numba
-from numba import jit, njit, prange
+from numba import njit, prange
 from numba.extending import overload
 import dask
 import dask.array as da
@@ -17,9 +16,7 @@
 from daskms.experimental.zarr import xds_from_zarr
 from omegaconf import ListConfig
 from skimage.morphology import label
-from scipy.optimize import curve_fit, fmin_l_bfgs_b
-from collections import namedtuple
-from africanus.coordinates.coordinates import radec_to_lmn
+from scipy.optimize import fmin_l_bfgs_b
 import xarray as xr
 from scipy.interpolate import RegularGridInterpolator
 from scipy.linalg import solve_triangular
diff --git a/pfb/utils/naming.py b/pfb/utils/naming.py
index cd1a7887..abe2376d 100644
--- a/pfb/utils/naming.py
+++ b/pfb/utils/naming.py
@@ -4,7 +4,6 @@
 import pickle
 import xarray as xr
 import concurrent.futures as cf
-import time
 
 def set_output_names(opts):
     '''
diff --git a/pfb/utils/restoration.py b/pfb/utils/restoration.py
index ccf4cf7e..32aa827b 100644
--- a/pfb/utils/restoration.py
+++ b/pfb/utils/restoration.py
@@ -1,9 +1,7 @@
-from functools import partial
 import numpy as np
-from pfb.utils.misc import Gaussian2D, convolve2gaussres, fitcleanbeam
+from pfb.utils.misc import convolve2gaussres
 from pfb.utils.fits import set_wcs, add_beampars, save_fits
 from pfb.utils.naming import xds_from_list
-from pfb.opt.pcg import pcg
 
 
 def restore_cube(ds_name,
diff --git a/pfb/utils/stokes.py b/pfb/utils/stokes.py
index 1d5c3082..8f43a157 100644
--- a/pfb/utils/stokes.py
+++ b/pfb/utils/stokes.py
@@ -218,10 +218,4 @@ def vfunc(gp, gq, W, V):
     else:
         raise ValueError(f"Jones term has incorrect number of dimensions")
 
-    # import inspect
-    # print(inspect.getsource(Djfn))
-    # print(inspect.getsource(Wjfn))
-
-    # quit()
-
     return vfunc, wfunc
diff --git a/pfb/utils/stokes2im.py b/pfb/utils/stokes2im.py
index 8dbaa20e..5356367b 100644
--- a/pfb/utils/stokes2im.py
+++ b/pfb/utils/stokes2im.py
@@ -4,7 +4,6 @@
 from distributed import worker_client
 import xarray as xr
 from uuid import uuid4
-from pfb.utils.stokes import stokes_funcs
 from pfb.utils.weighting import (_compute_counts, counts_to_weights,
                                  weight_data, filter_extreme_counts)
 from pfb.utils.fits import set_wcs, save_fits
@@ -12,9 +11,7 @@
 from ducc0.wgridder.experimental import vis2dirty
 from casacore.quanta import quantity
 from datetime import datetime
-from ducc0.fft import c2r, r2c, good_size
 from ducc0.misc import resize_thread_pool
-from africanus.constants import c as lightspeed
 from pfb.utils.astrometry import get_coordinates
 import gc
 iFs = np.fft.ifftshift
diff --git a/pfb/utils/stokes2vis.py b/pfb/utils/stokes2vis.py
index f2941fb9..cafb411a 100644
--- a/pfb/utils/stokes2vis.py
+++ b/pfb/utils/stokes2vis.py
@@ -1,17 +1,11 @@
 import numpy as np
 import numexpr as ne
 import xarray as xr
-import numba
-from numba import njit, prange, literally
+from numba import literally
 from dask.graph_manipulation import clone
 from distributed import worker_client
-import dask.array as da
-from xarray import Dataset
 from operator import getitem
 from pfb.utils.beam import interp_beam
-from pfb.utils.misc import weight_from_sigma, combine_columns
-import dask
-from pfb.utils.stokes import stokes_funcs
 from pfb.utils.weighting import weight_data
 from uuid import uuid4
 import gc
@@ -105,7 +99,6 @@ def single_stokes(
 
     nrow, nchan, ncorr = data.shape
 
-
     if opts.sigma_column is not None:
         weight = ne.evaluate('1.0/sigma**2',
                              local_dict={'sigma': getattr(ds, opts.sigma_column).values})
@@ -300,7 +293,6 @@ def single_stokes(
 
 
     coords = {'chan': (('chan',), freq),
-            #   'time': (('time',), utime),
               'l_beam': (('l_beam',), l_beam),
               'm_beam': (('m_beam',), m_beam)
     }
@@ -308,7 +300,6 @@ def single_stokes(
     unix_time = quantity(f'{time_out}s').to_unix_time()
     utc = datetime.utcfromtimestamp(unix_time).strftime('%Y-%m-%d %H:%M:%S')
 
-    # TODO - provide time and freq centroids
     attrs = {
         'ra' : radec[0],
         'dec': radec[1],
diff --git a/pfb/utils/weighting.py b/pfb/utils/weighting.py
index 0fbb968d..0c19ddde 100644
--- a/pfb/utils/weighting.py
+++ b/pfb/utils/weighting.py
@@ -2,8 +2,6 @@
 import numpy as np
 from numba import njit, prange, literally
 from numba.extending import overload
-import dask.array as da
-from ducc0.fft import c2c
 from africanus.constants import c as lightspeed
 from pfb.utils.misc import JIT_OPTIONS
 from pfb.utils.stokes import stokes_funcs
diff --git a/pfb/workers/degrid.py b/pfb/workers/degrid.py
index 24d026c8..50323718 100644
--- a/pfb/workers/degrid.py
+++ b/pfb/workers/degrid.py
@@ -1,7 +1,5 @@
 # flake8: noqa
-from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 import time
 from omegaconf import OmegaConf
 import pyscilog
@@ -104,20 +102,15 @@ def _degrid(**kw):
     import numpy as np
     from pfb.utils.misc import construct_mappings
     import dask
-    from dask.distributed import performance_report, wait, get_client
+    from dask.distributed import wait, get_client
     from dask.graph_manipulation import clone
-    from daskms.experimental.zarr import xds_from_zarr
     from daskms import xds_from_storage_ms as xds_from_ms
-    from daskms import xds_from_storage_table as xds_from_table
     from daskms import xds_to_storage_table as xds_to_table
     from daskms.fsspec_store import DaskMSStore
     import dask.array as da
-    from africanus.constants import c as lightspeed
-    from africanus.gridding.wgridder.dask import model as im2vis
-    from pfb.operators.gridder import comps2vis, _comps2vis_impl
-    from pfb.utils.fits import load_fits, data_from_header, set_wcs
+    from pfb.operators.gridder import comps2vis
+    from pfb.utils.fits import set_wcs
     from regions import Regions
-    from astropy.io import fits
     from pfb.utils.naming import xds_from_url
     import xarray as xr
     import sympy as sm
diff --git a/pfb/workers/fluxtractor.py b/pfb/workers/fluxtractor.py
index 7698b9c4..85dc6d6d 100644
--- a/pfb/workers/fluxtractor.py
+++ b/pfb/workers/fluxtractor.py
@@ -1,7 +1,6 @@
 # flake8: noqa
 from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
@@ -38,7 +37,7 @@ def fluxtractor(**kw):
     OmegaConf.set_struct(opts, True)
 
     from pfb import set_envs
-    from ducc0.misc import resize_thread_pool, thread_pool_size
+    from ducc0.misc import resize_thread_pool
     resize_thread_pool(opts.nthreads)
     set_envs(opts.nthreads, ncpu)
 
@@ -53,7 +52,6 @@ def fluxtractor(**kw):
     for key in opts.keys():
         print('     %25s = %s' % (key, opts[key]), file=log)
 
-    from daskms.fsspec_store import DaskMSStore
     from pfb.utils.naming import xds_from_url
 
     basename = opts.output_filename
@@ -138,10 +136,9 @@ def _fluxtractor(**kw):
     from itertools import cycle
     import numpy as np
     import xarray as xr
-    from pfb.utils.fits import load_fits, set_wcs, save_fits
+    from pfb.utils.fits import load_fits, set_wcs
     from daskms.fsspec_store import DaskMSStore
-    from pfb.utils.naming import xds_from_url, xds_from_list
-    from pfb.utils.misc import init_mask, dds2cubes
+    from pfb.utils.naming import xds_from_url
     from pfb.opt.pcg import pcg_dds
     from ducc0.misc import resize_thread_pool, thread_pool_size
     from ducc0.fft import c2c
@@ -149,10 +146,6 @@ def _fluxtractor(**kw):
     Fs = np.fft.fftshift
 
     basename = opts.output_filename
-    if opts.fits_output_folder is not None:
-        fits_oname = opts.fits_output_folder + '/' + basename.split('/')[-1]
-    else:
-        fits_oname = basename
 
     dds_name = f'{basename}_{opts.suffix}.dds'
     dds_store = DaskMSStore(dds_name)
diff --git a/pfb/workers/grid.py b/pfb/workers/grid.py
index e9652ff1..dbbc30ba 100644
--- a/pfb/workers/grid.py
+++ b/pfb/workers/grid.py
@@ -1,8 +1,5 @@
 # flake8: noqa
-import os
-from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
diff --git a/pfb/workers/hci.py b/pfb/workers/hci.py
index 2b8da7b7..9cab8e4c 100644
--- a/pfb/workers/hci.py
+++ b/pfb/workers/hci.py
@@ -1,6 +1,4 @@
 # flake8: noqa
-import os
-import sys
 from pfb.workers.main import cli
 from omegaconf import OmegaConf
 import pyscilog
diff --git a/pfb/workers/init.py b/pfb/workers/init.py
index 503053dd..f7000fcc 100644
--- a/pfb/workers/init.py
+++ b/pfb/workers/init.py
@@ -1,6 +1,4 @@
 # flake8: noqa
-import os
-import sys
 from pfb.workers.main import cli
 from omegaconf import OmegaConf
 import pyscilog
diff --git a/pfb/workers/kclean.py b/pfb/workers/kclean.py
index 99c00bf8..34a432d6 100644
--- a/pfb/workers/kclean.py
+++ b/pfb/workers/kclean.py
@@ -1,9 +1,5 @@
 # flake8: noqa
-import os
-import sys
-from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
@@ -34,7 +30,7 @@ def kclean(**kw):
     OmegaConf.set_struct(opts, True)
 
     from pfb import set_envs
-    from ducc0.misc import resize_thread_pool, thread_pool_size
+    from ducc0.misc import resize_thread_pool
     resize_thread_pool(opts.nthreads)
     set_envs(opts.nthreads, ncpu)
 
diff --git a/pfb/workers/model2comps.py b/pfb/workers/model2comps.py
index a007110e..e3634144 100644
--- a/pfb/workers/model2comps.py
+++ b/pfb/workers/model2comps.py
@@ -1,7 +1,6 @@
 # flake8: noqa
 from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
diff --git a/pfb/workers/restore.py b/pfb/workers/restore.py
index 44637267..b120b30c 100644
--- a/pfb/workers/restore.py
+++ b/pfb/workers/restore.py
@@ -1,7 +1,6 @@
 # flake8: noqa
 from contextlib import ExitStack
 from pfb.workers.main import cli
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
diff --git a/pfb/workers/sara.py b/pfb/workers/sara.py
index d3cd784a..cce03e5a 100644
--- a/pfb/workers/sara.py
+++ b/pfb/workers/sara.py
@@ -1,6 +1,4 @@
 # flake8: noqa
-import os
-import sys
 import concurrent.futures as cf
 from pfb.workers.main import cli
 from omegaconf import OmegaConf
diff --git a/pfb/workers/smoovie.py b/pfb/workers/smoovie.py
index df5826f8..9e854a40 100644
--- a/pfb/workers/smoovie.py
+++ b/pfb/workers/smoovie.py
@@ -1,6 +1,3 @@
-import os
-from pathlib import Path
-from contextlib import ExitStack
 from pfb.workers.main import cli
 from omegaconf import OmegaConf
 import pyscilog
diff --git a/pfb/workers/spotless.py b/pfb/workers/spotless.py
index 8df5e747..c8aa299c 100644
--- a/pfb/workers/spotless.py
+++ b/pfb/workers/spotless.py
@@ -1,10 +1,6 @@
 # flake8: noqa
-import os
-import sys
 from contextlib import ExitStack
 from pfb.workers.main import cli
-from functools import partial
-import click
 from omegaconf import OmegaConf
 import pyscilog
 pyscilog.init('pfb')
@@ -33,7 +29,7 @@ def spotless(**kw):
         ncpu = ncpu//2
 
     from pfb import set_envs
-    from ducc0.misc import resize_thread_pool, thread_pool_size
+    from ducc0.misc import resize_thread_pool
     resize_thread_pool(opts.nthreads)
     set_envs(opts.nthreads, ncpu)