diff --git a/piff/basis_interp.py b/piff/basis_interp.py
index 9d2f6973..29cfdc98 100644
--- a/piff/basis_interp.py
+++ b/piff/basis_interp.py
@@ -24,6 +24,7 @@
 from .interp import Interp
 from .star import Star
 from . import _piff
+from .config import LoggerWrapper
 
 try: 
     import jax
@@ -154,7 +155,7 @@ def solve(self, stars, logger=None):
         :param stars:       A list of Star instances to interpolate between
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         if self.q is None:
             raise RuntimeError("Attempt to solve() before initialize() of BasisInterp")
 
@@ -210,7 +211,7 @@ def solve(self, stars, logger=None):
         # QR method often has fewer numerical problems than the direct method for large matrices.
         # It just requires a lot of memory.
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         if self.solver == "qr":
             self._solve_qr(stars, logger)
         else:
@@ -260,15 +261,15 @@ def _solve_qr(self, stars, logger):
         if cond < 1.e-12:
             # Note: this calculation is much slower, but it is safe to use even for
             # singular inputs, so it will always produce a valid answer.
-            logger.info('Nominal condition is %s (min, max = %s, %s)', cond,
-                        np.min(abs_Rdiag), np.max(abs_Rdiag))
-            logger.info('Switching to QRP solution')
+            logger.verbose('Nominal condition is %s (min, max = %s, %s)', cond,
+                           np.min(abs_Rdiag), np.max(abs_Rdiag))
+            logger.verbose('Switching to QRP solution')
             QR, P, tau, work, info = scipy.linalg.lapack.dgeqp3(A, overwrite_a=True)
             P[:] -= 1  # Switch to python 0-based indexing.
             abs_Rdiag = np.abs(np.diag(QR))
             cond = np.min(abs_Rdiag) / np.max(abs_Rdiag)
-            logger.info('Condition for QRP is %s (min, max = %s, %s)', cond,
-                        np.min(abs_Rdiag), np.max(abs_Rdiag))
+            logger.verbose('Condition for QRP is %s (min, max = %s, %s)', cond,
+                           np.min(abs_Rdiag), np.max(abs_Rdiag))
             # Skip any rows of R that have essentially 0 on the diagonal.
             k = np.sum(abs_Rdiag > 1.e-15 * np.max(abs_Rdiag))
             logger.debug('k = %d, m = %d',k,m)
@@ -363,7 +364,7 @@ def _solve_direct_python(self, stars, logger):
                     ATb += A1.T.dot(s.fit.b)
                     ATA += A1.T.dot(A1)
 
-        logger.info('Beginning solution of matrix size %s',ATA.shape)
+        logger.verbose('Beginning solution of matrix size %s',ATA.shape)
         try:
             with warnings.catch_warnings(record=True) as w:
                 warnings.simplefilter("always")
@@ -383,25 +384,25 @@ def _solve_direct_python(self, stars, logger):
                 # size rediduals.  We don't have a unit test that would catch this, so be careful
                 # about changing the behavior of this part of the code!  For now, we just go to
                 # the svd solution when ATA is fully singular.
-                logger.info(w[0].message)
+                logger.verbose(w[0].message)
                 logger.debug('norm(ATA dq - ATb) = %s',scipy.linalg.norm(ATA.dot(dq) - ATb))
                 logger.debug('norm(dq) = %s',scipy.linalg.norm(dq))
 
         except (np.linalg.LinAlgError, scipy.linalg.LinAlgError) as e:
-            logger.info('Caught %s',str(e))
-            logger.info('Switching to svd solution')
+            logger.verbose('Caught %s',str(e))
+            logger.verbose('Switching to svd solution')
             Sd,U = scipy.linalg.eigh(ATA)
             nsvd = np.sum(np.abs(Sd) > 1.e-15 * np.abs(Sd[-1]))
-            logger.info('2-condition is %e',np.abs(Sd[-1]/Sd[0]))
-            logger.info('nsvd = %d of %d',nsvd,len(Sd))
+            logger.verbose('2-condition is %e',np.abs(Sd[-1]/Sd[0]))
+            logger.verbose('nsvd = %d of %d',nsvd,len(Sd))
             # Note: unlike scipy.linalg.svd, the Sd here is in *ascending* order, not descending.
             Sd[-nsvd:] = 1./Sd[-nsvd:]
             Sd[:-nsvd] = 0.
             S = np.diag(Sd)
             dq = U.dot(S.dot(U.T.dot(ATb)))
-            logger.info('norm(ATA dq - ATb) = %s',scipy.linalg.norm(ATA.dot(dq) - ATb))
-            logger.info('norm(dq) = %s',scipy.linalg.norm(dq))
-            logger.info('norm(q) = %s',scipy.linalg.norm(self.q))
+            logger.verbose('norm(ATA dq - ATb) = %s',scipy.linalg.norm(ATA.dot(dq) - ATb))
+            logger.verbose('norm(dq) = %s',scipy.linalg.norm(dq))
+            logger.verbose('norm(q) = %s',scipy.linalg.norm(self.q))
 
         logger.debug('...finished solution')
         # Reshape dq back into a 2d array and add it to the current solution.
@@ -483,7 +484,7 @@ def __init__(
         ):
         super(BasisPolynomial, self).__init__()
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
 
         self._keys = keys
         if hasattr(order,'__len__'):
@@ -515,7 +516,7 @@ def __init__(
             self.solver = "qr"
 
         if not CAN_USE_JAX and self.solver == "jax":
-            logger.warning("JAX not installed. Reverting to numpy/scipy.")
+            logger.info("JAX not installed. Reverting to numpy/scipy.")
             self.solver = "scipy"
 
         if self._max_order<0 or np.any(np.array(self._orders) < 0):
diff --git a/piff/config.py b/piff/config.py
index 781be187..6eb30c6e 100644
--- a/piff/config.py
+++ b/piff/config.py
@@ -19,6 +19,30 @@
 import yaml
 import os
 import galsim
+import logging
+
+
+# Add a logging level between INFO and DEBUG
+VERBOSE = (logging.INFO + logging.DEBUG) // 2
+logging.addLevelName(VERBOSE, "VERBOSE")
+
+class PiffLogger(logging.Logger):
+
+    def verbose(self, fmt, *args, **kwargs):
+        """Issue a VERBOSE level log message.
+
+        Arguments are as for `logging.info`.
+        ``VERBOSE`` is between ``DEBUG`` and ``INFO``.
+        """
+        # Based on the LSST logger at
+        # https://github.com/lsst/utils/blob/main/python/lsst/utils/logging.py
+        self.log(VERBOSE, fmt, *args, **kwargs)
+
+class LoggerWrapper(galsim.config.LoggerWrapper):
+    # Extend GalSim's wrapper to include the verbose level
+    def verbose(self, *args, **kwargs):
+        if self.logger and self.isEnabledFor(VERBOSE):
+            self.logger.verbose(*args, **kwargs)
 
 def setup_logger(verbose=1, log_file=None):
     """Build a logger object to use for logging progress
@@ -32,16 +56,18 @@ def setup_logger(verbose=1, log_file=None):
 
     :returns: a logging.Logger instance
     """
-    import logging
+    # Verbose logging is midway between INFO and DEBUG.
+
     # Parse the integer verbosity level from the command line args into a logging_level string
-    logging_levels = { 0: logging.CRITICAL,
-                       1: logging.WARNING,
-                       2: logging.INFO,
+    logging_levels = { 0: logging.ERROR,
+                       0.5: logging.WARNING,
+                       1: logging.INFO,
+                       2: VERBOSE,
                        3: logging.DEBUG }
     logging_level = logging_levels[verbose]
 
     # Setup logging to go to sys.stdout or (if requested) to an output file
-    logger = logging.getLogger('piff')
+    logger = PiffLogger(logging.getLogger('piff'))
     logger.handlers = []  # Remove any existing handlers
     if log_file is None:
         handle = logging.StreamHandler()
@@ -54,7 +80,6 @@ def setup_logger(verbose=1, log_file=None):
 
     return logger
 
-
 def parse_variables(config, variables, logger):
     """Parse configuration variables and add them to the config dict
 
@@ -78,8 +103,8 @@ def parse_variables(config, variables, logger):
             # Use YAML parser to evaluate the string in case it is a list for instance.
             value = yaml.safe_load(value)
         except yaml.YAMLError as e:  # pragma: no cover
-            logger.warning('Caught %r',e)
-            logger.warning('Unable to parse %s.  Treating it as a string.',value)
+            logger.info('Caught %r',e)
+            logger.info('Unable to parse %s.  Treating it as a string.',value)
         new_params[key] = value
     galsim.config.UpdateConfig(config, new_params)
 
@@ -201,7 +226,7 @@ def plotify(config, logger=None):
     file_name = config['output']['file_name']
     if 'dir' in config['output']:
         file_name = os.path.join(config['output']['dir'], file_name)
-    logger.info("Looking for PSF at %s", file_name)
+    logger.verbose("Looking for PSF at %s", file_name)
     psf = PSF.read(file_name, logger=logger)
 
     # The stars we care about for plotify are psf.stars, not what would be made from scratch by
@@ -296,7 +321,7 @@ def meanify(config, logger=None):
     config['output']['file_name'] = psf_list
 
     file_name_in = config['output']['file_name']
-    logger.info("Looking for PSF at %s", file_name_in)
+    logger.verbose("Looking for PSF at %s", file_name_in)
 
     file_name_out = config['hyper']['file_name']
     if 'dir' in config['hyper']:
@@ -316,7 +341,7 @@ def meanify(config, logger=None):
 
     params = np.concatenate(params, axis=0)
     coords = np.concatenate(coords, axis=0)
-    logger.info('Computing average for {0} params with {1} stars'.format(len(params[0]), len(coords)))
+    logger.verbose('Computing average for {0} params with {1} stars'.format(len(params[0]), len(coords)))
 
     if params_fitted is None:
         params_fitted = range(len(params[0]))
@@ -361,6 +386,6 @@ def meanify(config, logger=None):
     data['COORDS0'] = coords0
     data['PARAMS0'] = params0
 
-    logger.info('Writing average solution to {0}'.format(file_name_out))
+    logger.verbose('Writing average solution to {0}'.format(file_name_out))
     with fitsio.FITS(file_name_out,'rw',clobber=True) as f:
         f.write_table(data, extname='average_solution')
diff --git a/piff/convolvepsf.py b/piff/convolvepsf.py
index bc01f8fd..f1a38895 100644
--- a/piff/convolvepsf.py
+++ b/piff/convolvepsf.py
@@ -173,7 +173,7 @@ def single_iteration(self, stars, logger, convert_funcs, draw_method):
 
         # Fit each component in order
         for k, comp in enumerate(self.components):
-            logger.info("Starting work on component %d (%s)", k, comp._type_name)
+            logger.verbose("Starting work on component %d (%s)", k, comp._type_name)
 
             # Update the convert_funcs to add a convolution by the other components.
             new_convert_funcs = []
@@ -276,7 +276,8 @@ def _finish_write(self, writer, logger):
         :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         chisq_dict = {
             'chisq' : self.chisq,
             'last_delta_chisq' : self.last_delta_chisq,
diff --git a/piff/gp_interp.py b/piff/gp_interp.py
index b5ce2f87..bf670e9a 100644
--- a/piff/gp_interp.py
+++ b/piff/gp_interp.py
@@ -371,6 +371,7 @@ class GPInterpDepr(GPInterp):
     _type_name = 'GPInterp'
 
     def __init__(self, *args, logger=None, **kwargs):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.error("WARNING: The name GPInterp is deprecated. Use GP or GaussianProcess instead.")
         super().__init__(*args, logger=logger, **kwargs)
diff --git a/piff/gsobject_model.py b/piff/gsobject_model.py
index 0e8b1be4..5c57ddf2 100644
--- a/piff/gsobject_model.py
+++ b/piff/gsobject_model.py
@@ -260,7 +260,9 @@ def least_squares_fit(self, star, logger=None, convert_func=None, draw_method=No
 
         :returns: (flux, dx, dy, scale, g1, g2, flag)
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        import warnings
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.debug("Start least_squares")
         params = self._get_params(star)
 
@@ -270,9 +272,11 @@ def least_squares_fit(self, star, logger=None, convert_func=None, draw_method=No
         if params[3]**2 < star.data.local_wcs.pixelArea():
             params[3] = star.data.local_wcs.pixelArea()**0.5
 
-        results = scipy.optimize.least_squares(self._resid, params,
-                                               args=(star,convert_func,draw_method),
-                                               **self._scipy_kwargs)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            results = scipy.optimize.least_squares(self._resid, params,
+                                                   args=(star,convert_func,draw_method),
+                                                   **self._scipy_kwargs)
         #logger.debug(results)
         logger.debug("results = %s",results.x)
         if not results.success:
@@ -306,7 +310,8 @@ def fit(self, star, fastfit=None, logger=None, convert_func=None, draw_method=No
 
         :returns: a new Star with the fitted parameters in star.fit
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         if fastfit is None:
             fastfit = self._fastfit
         if convert_func is not None:
@@ -370,7 +375,8 @@ def initialize(self, star, logger=None, default_init=None):
 
         :returns: a new initialized Star.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         init = self._init if self._init is not None else default_init
         if init is None: init = 'hsm'
         logger.debug("initializing GSObject with method %s",init)
@@ -385,7 +391,7 @@ def initialize(self, star, logger=None, default_init=None):
             flux_scaling = 1.e-6
             # Also, make sure fit_flux=True.  Otherwise, this won't work properly.
             if not self._fit_flux:
-                logger.info('Setting fit_flux=True, since init=zero')
+                logger.verbose('Setting fit_flux=True, since init=zero')
                 self._fit_flux = True
                 self.kwargs['fit_flux'] = True
         elif init == 'delta':
@@ -394,7 +400,7 @@ def initialize(self, star, logger=None, default_init=None):
             flux_scaling, size_scaling = init
             size *= size_scaling
             if not self._fit_flux:
-                logger.info(f'Setting fit_flux=True, since init={init}')
+                logger.verbose(f'Setting fit_flux=True, since init={init}')
                 self._fit_flux = True
                 self.kwargs['fit_flux'] = True
         elif init.startswith('(') and init.endswith(')'):
@@ -520,6 +526,7 @@ class GSObjectModelDepr(GSObjectModel):
     _type_name = 'GSObjectModel'
 
     def __init__(self, *args, logger=None, **kwargs):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.error("WARNING: The name GSObjectModel is deprecated. Use GSObject instead.")
         super().__init__(*args, logger=logger, **kwargs)
diff --git a/piff/input.py b/piff/input.py
index 444bea52..267bd044 100644
--- a/piff/input.py
+++ b/piff/input.py
@@ -102,7 +102,8 @@ def makeStars(self, logger=None):
 
         :returns: a list of Star instances
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         if self.nimages == 1:
             logger.debug("Making star list")
@@ -122,8 +123,8 @@ def makeStars(self, logger=None):
         # Concatenate the star lists into a single list
         stars = [s for slist in all_stars if slist is not None for s in slist if slist]
 
-        logger.warning("Read a total of %d objects from %d image%s",len(stars),self.nimages,
-                       "s" if self.nimages > 1 else "")
+        logger.info("Read a total of %d objects from %d image%s",len(stars),self.nimages,
+                    "s" if self.nimages > 1 else "")
 
         return stars
 
@@ -268,7 +269,8 @@ class InputFiles(Input):
 
     def __init__(self, config, logger=None):
         import copy
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         req = { 'image_file_name': str,
                 'cat_file_name': str,
@@ -412,7 +414,7 @@ def __init__(self, config, logger=None):
         if cat_list is not None:
             logger.debug('cat_list = %s',cat_list)
             if len(cat_list) == 1 and nimages > 1:
-                logger.info("Using the same catlist for all image")
+                logger.verbose("Using the same catlist for all image")
                 cat_list = cat_list * nimages
             elif nimages != len(cat_list):
                 raise ValueError("nimages = %s doesn't match length of cat_file_name list (%d)"%(
@@ -433,7 +435,7 @@ def __init__(self, config, logger=None):
         self.remove_signal_from_weight = config.get('remove_signal_from_weight', False)
         self.invert_weight = config.get('invert_weight', False)
 
-        logger.info("Reading in %d images",nimages)
+        logger.verbose("Reading in %d images",nimages)
         for image_num in range(nimages):
 
             # This changes for each input image.
@@ -540,7 +542,8 @@ def load_images(self, stars, logger=None):
 
         :returns: a new list of Stars with the images information loaded.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         images = {}
         weight_images = {}
@@ -570,7 +573,8 @@ def getRawImageData(self, image_num, logger=None):
     def _getRawImageData(image_kwargs, cat_kwargs, wcs,
                          invert_weight, remove_signal_from_weight,
                          logger=None):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         image, weight = InputFiles.readImage(logger=logger, **image_kwargs)
 
         if invert_weight:
@@ -599,8 +603,8 @@ def _getRawImageData(image_kwargs, cat_kwargs, wcs,
                 # If no gain provided, then we want to estimate the gain from the weight image.
                 weight, g = InputFiles._removeSignalFromWeight(signal, weight)
                 extra_props['gain'] = np.array([g] * len(image_pos), dtype=float)
-                logger.warning("Empirically determined gain = %f",g)
-            logger.info("Removed signal from weight image.")
+                logger.info("Empirically determined gain = %f",g)
+            logger.verbose("Removed signal from weight image.")
 
         return image, weight, image_pos, extra_props
 
@@ -613,7 +617,7 @@ def _makeStarsFromImage(image_kwargs, cat_kwargs, wcs, chipnum,
         """
         image, wt, image_pos, extra_props = InputFiles._getRawImageData(
                 image_kwargs, cat_kwargs, wcs, invert_weight, remove_signal_from_weight, logger)
-        logger.info("Processing catalog %s with %d objects",chipnum,len(image_pos))
+        logger.verbose("Processing catalog %s with %d objects",chipnum,len(image_pos))
 
         objects = []
 
@@ -636,11 +640,11 @@ def _makeStarsFromImage(image_kwargs, cat_kwargs, wcs, chipnum,
                     logger.warning("Skipping this object.")
                     continue
                 if use_partial:
-                    logger.info("Object at position %f,%f overlaps the edge of the image.  "
+                    logger.verbose("Object at position %f,%f overlaps the edge of the image.  "
                                 "Using smaller than the full stamp size: %s", x, y, bounds)
                 else:
-                    logger.warning("Object at position %f,%f overlaps the edge of the image.", x, y)
-                    logger.warning("Skipping this object.")
+                    logger.info("Object at position %f,%f overlaps the edge of the image.", x, y)
+                    logger.info("Skipping this object.")
                     continue
             stamp = image[bounds].copy()
             wt_stamp = wt[bounds].copy()
@@ -663,9 +667,9 @@ def _makeStarsFromImage(image_kwargs, cat_kwargs, wcs, chipnum,
             if 'satur' in props:
                 max_val = np.max(stamp.array)
                 if max_val > props['satur']:
-                    logger.warning("Object at position %f,%f has saturated pixels.", x, y)
-                    logger.warning("Maximum value is %f.", max_val)
-                    logger.warning("Skipping this object.")
+                    logger.info("Object at position %f,%f has saturated pixels.", x, y)
+                    logger.info("Maximum value is %f.", max_val)
+                    logger.info("Skipping this object.")
                     continue
 
             # Subtract the sky
@@ -693,11 +697,11 @@ def setWCS(self, config, logger):
             image_hdu = kwargs['image_hdu']
             image = galsim.fits.read(image_file_name, hdu=image_hdu)
             if 'wcs' in config:
-                logger.warning("Using custom wcs from config for %s",image_file_name)
+                logger.info("Using custom wcs from config for %s",image_file_name)
                 base = { 'input' : config, 'index_key' : 'image_num', 'image_num' : image_num }
                 wcs = galsim.config.BuildWCS(config, 'wcs', base, logger)
             else:
-                logger.warning("Getting wcs from image file %s",image_file_name)
+                logger.info("Getting wcs from image file %s",image_file_name)
                 wcs = image.wcs
             self.wcs_list.append(wcs)
             self.center_list.append(image.true_center)
@@ -747,26 +751,26 @@ def readImage(image_file_name, image_hdu, weight_hdu, badpix_hdu,
         :returns: image, weight
         """
         # Read in the image
-        logger.warning("Reading image file %s",image_file_name)
+        logger.info("Reading image file %s",image_file_name)
         image = galsim.fits.read(image_file_name, hdu=image_hdu)
 
         # If requested, subtract a sky image
         if sky_file_name is not None:
             hdu_str = " from hdu %s"%(sky_hdu) if sky_hdu is not None else ""
-            logger.info("Reading sky image %s.", sky_file_name + hdu_str)
+            logger.verbose("Reading sky image %s.", sky_file_name + hdu_str)
             sky = galsim.fits.read(sky_file_name, hdu=sky_hdu)
             image -= sky
 
         # Either read in the weight image, or build a dummy one
         if weight_hdu is not None:
-            logger.info("Reading weight image from hdu %d.", weight_hdu)
+            logger.verbose("Reading weight image from hdu %d.", weight_hdu)
             weight = galsim.fits.read(image_file_name, hdu=weight_hdu)
             if np.all(weight.array == 0):
                 logger.error("According to the weight mask in %s, all pixels have zero weight!",
                              image_file_name)
             if np.any(weight.array < 0):
-                logger.error("Warning: weight map has invalid negative-valued pixels. "+
-                             "Taking them to be 0.0")
+                logger.warning("Warning: weight map has invalid negative-valued pixels. "+
+                               "Taking them to be 0.0")
                 weight.array[weight.array < 0] = 0.
         elif noise is not None:
             logger.debug("Making uniform weight image based on noise variance = %f", noise)
@@ -777,7 +781,7 @@ def readImage(image_file_name, image_hdu, weight_hdu, badpix_hdu,
 
         # If requested, set wt=0 for any bad pixels
         if badpix_hdu is not None:
-            logger.info("Reading badpix image from hdu %d.", badpix_hdu)
+            logger.verbose("Reading badpix image from hdu %d.", badpix_hdu)
             badpix = galsim.fits.read(image_file_name, hdu=badpix_hdu)
             # The badpix image may be offset by 32768 from the true value.
             # If so, subtract it off.
@@ -855,11 +859,11 @@ def readStarCatalog(cat_file_name, cat_hdu, x_col, y_col,
 
         if cat_file_name is None:
             # This is possible e.g. when loading images into an existing list of star instances.
-            logger.warning("Not reading input catalog.")
+            logger.info("Not reading input catalog.")
             return None, None
 
         # Read in the star catalog
-        logger.warning("Reading input catalog %s.",cat_file_name)
+        logger.info("Reading input catalog %s.",cat_file_name)
         cat = fitsio.read(cat_file_name, cat_hdu)
 
         if flag_col is not None:
@@ -867,27 +871,27 @@ def readStarCatalog(cat_file_name, cat_hdu, x_col, y_col,
                 raise ValueError("flag_col = %s is not a column in %s"%(flag_col,cat_file_name))
             col = cat[flag_col]
             if len(col.shape) == 2:
-                logger.warning("Flag col (%s) is multidimensional.  Treating as an array of bool",
-                               flag_col)
+                logger.info("Flag col (%s) is multidimensional.  Treating as an array of bool",
+                            flag_col)
             if use_flag is not None:
                 # Remove any objects with flag & use_flag == 0
                 mask = InputFiles._flag_select(col, use_flag) == 0
-                logger.info("Removing objects with flag (col %s) & %d == 0",flag_col,use_flag)
+                logger.verbose("Removing objects with flag (col %s) & %d == 0",flag_col,use_flag)
                 if skip_flag != -1:
                     mask |= InputFiles._flag_select(col, skip_flag) != 0
-                    logger.info("Removing objects with flag (col %s) & %d != 0",flag_col,skip_flag)
+                    logger.verbose("Removing objects with flag (col %s) & %d != 0",flag_col,skip_flag)
             else:
                 # Remove any objects with flag & skip_flag != 0
                 mask = InputFiles._flag_select(col, skip_flag) != 0
                 if skip_flag == -1:
-                    logger.info("Removing objects with flag (col %s) != 0",flag_col)
+                    logger.verbose("Removing objects with flag (col %s) != 0",flag_col)
                 else:
-                    logger.info("Removing objects with flag (col %s) & %d != 0",flag_col,skip_flag)
+                    logger.verbose("Removing objects with flag (col %s) & %d != 0",flag_col,skip_flag)
             cat = cat[mask == 0]
 
         # Limit to nstars objects
         if nstars is not None and nstars < len(cat):
-            logger.info("Limiting to %d objects for %s",nstars,cat_file_name)
+            logger.verbose("Limiting to %d objects for %s",nstars,cat_file_name)
             cat = cat[:nstars]
 
         # Make the list of positions:
@@ -1044,7 +1048,8 @@ def setPointing(self, ra, dec, logger=None):
            individual WCS functions. [Not implemented currently.]
         """
         import fitsio
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         if (ra is None) != (dec is None):
             raise ValueError("Only one of ra, dec was specified")
@@ -1056,9 +1061,9 @@ def setPointing(self, ra, dec, logger=None):
                     wcs = self.wcs_list[0]
                     center = self.center_list[0]
                     self.pointing = wcs.toWorld(center)
-                    logger.info("Setting pointing to image center: %.3f h, %.3f d",
-                                self.pointing.ra / galsim.hours,
-                                self.pointing.dec / galsim.degrees)
+                    logger.verbose("Setting pointing to image center: %.3f h, %.3f d",
+                                   self.pointing.ra / galsim.hours,
+                                   self.pointing.dec / galsim.degrees)
                 else:
                     # Use the mean of all the image centers
                     plist = [wcs.toWorld(center)
@@ -1069,18 +1074,18 @@ def setPointing(self, ra, dec, logger=None):
                     y = np.mean(ylist)
                     z = np.mean(zlist)
                     self.pointing = galsim.CelestialCoord.from_xyz(x,y,z)
-                    logger.info("Setting pointing to mean of image centers: %.3f h, %.3f d",
-                                self.pointing.ra / galsim.hours,
-                                self.pointing.dec / galsim.degrees)
+                    logger.verbose("Setting pointing to mean of image centers: %.3f h, %.3f d",
+                                   self.pointing.ra / galsim.hours,
+                                   self.pointing.dec / galsim.degrees)
             else:
                 self.pointing = None
         elif type(ra) in [float, int]:
             ra = float(ra) * galsim.hours
             dec = float(dec) * galsim.degrees
             self.pointing = galsim.CelestialCoord(ra,dec)
-            logger.info("Setting pointing to: %.3f h, %.3f d",
-                        self.pointing.ra / galsim.hours,
-                        self.pointing.dec / galsim.degrees)
+            logger.verbose("Setting pointing to: %.3f h, %.3f d",
+                           self.pointing.ra / galsim.hours,
+                           self.pointing.dec / galsim.degrees)
         elif str(ra) != ra:
             raise ValueError("Unable to parse input ra: %s"%ra)
         elif str(dec) != dec:
@@ -1089,15 +1094,15 @@ def setPointing(self, ra, dec, logger=None):
             ra = galsim.Angle.from_hms(ra)
             dec = galsim.Angle.from_dms(dec)
             self.pointing = galsim.CelestialCoord(ra,dec)
-            logger.info("Setting pointing to: %.3f h, %.3f d",
-                        self.pointing.ra / galsim.hours,
-                        self.pointing.dec / galsim.degrees)
+            logger.verbose("Setting pointing to: %.3f h, %.3f d",
+                           self.pointing.ra / galsim.hours,
+                           self.pointing.dec / galsim.degrees)
         else:
             file_name = self.image_file_name[0]
             if len(self.chipnums) == 1:
-                logger.info("Setting pointing from keywords %s, %s", ra, dec)
+                logger.verbose("Setting pointing from keywords %s, %s", ra, dec)
             else:
-                logger.info("Setting pointing from keywords %s, %s in %s", ra, dec, file_name)
+                logger.verbose("Setting pointing from keywords %s, %s in %s", ra, dec, file_name)
             fits = fitsio.FITS(file_name)
             hdu = 1 if file_name.endswith('.fz') else 0
             header = fits[hdu].read_header()
diff --git a/piff/knn_interp.py b/piff/knn_interp.py
index 7e1966c2..a83a2788 100644
--- a/piff/knn_interp.py
+++ b/piff/knn_interp.py
@@ -77,7 +77,8 @@ def _fit(self, locations, targets, logger=None):
         :param targets:     The target values. (n_samples, n_targets).
                             (In sklearn parlance, this is 'y'.)
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         self.knn.fit(locations, targets)
         self.locations = locations
         logger.debug('locations updated to shape: %s', self.locations.shape)
@@ -92,7 +93,8 @@ def _predict(self, locations, logger=None):
 
         :returns:   Regressed parameters y (n_samples, n_targets)
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         regression = self.knn.predict(locations)
         logger.debug('Regression shape: %s', regression.shape)
         return regression
@@ -147,7 +149,8 @@ def interpolateList(self, stars, logger=None):
 
         :returns: a list of new Star instances with interpolated parameters
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         locations = np.array([self.getProperties(star) for star in stars])
         targets = self._predict(locations)
         stars_fitted = []
@@ -192,7 +195,8 @@ class kNNInterp(KNNInterp):
     _type_name = 'kNNInterp'
 
     def __init__(self, *args, logger=None, **kwargs):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.error("WARNING: The name kNNInterp is deprecated. "
                      "Use KNN or KNearestNeighbors instead.")
         super().__init__(*args, logger=logger, **kwargs)
diff --git a/piff/optical_model.py b/piff/optical_model.py
index 2db1ad39..2e745061 100644
--- a/piff/optical_model.py
+++ b/piff/optical_model.py
@@ -169,7 +169,8 @@ class Optical(Model):
     _model_can_be_offset = False
 
     def __init__(self, template=None, gsparams=None, atmo_type='VonKarman', logger=None, **kwargs):
-        self.logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        self.logger = LoggerWrapper(logger)
         self.set_num(None)
 
         # If pupil_angle and strut angle are provided as strings, eval them.
diff --git a/piff/outliers.py b/piff/outliers.py
index 9f0b10fe..afb2001e 100644
--- a/piff/outliers.py
+++ b/piff/outliers.py
@@ -22,7 +22,6 @@
 import galsim
 from scipy.stats import chi2
 
-
 class Outliers(object):
     """The base class for handling outliers.
 
@@ -247,6 +246,8 @@ class ChisqOutliers(Outliers):
 
     def __init__(self, thresh=None, ndof=None, prob=None, nsigma=None, max_remove=None,
                  include_reserve=None, logger=None):
+        from .config import LoggerWrapper
+
         if all( (thresh is None, ndof is None, prob is None, nsigma is None) ):
             raise TypeError("One of thresh, ndof, prob, or nsigma is required.")
         if thresh is not None and any( (ndof is not None, prob is not None, nsigma is not None) ):
@@ -256,7 +257,7 @@ def __init__(self, thresh=None, ndof=None, prob=None, nsigma=None, max_remove=No
         if prob is not None and nsigma is not None:
             raise TypeError("Only one of thresh, ndof, prob, or nsigma may be given.")
         if include_reserve is not None:
-            logger = galsim.config.LoggerWrapper(logger)
+            logger = LoggerWrapper(logger)
             logger.error("WARNING: include_reserve is no longer used.")
 
         # The only one of these we can convert now is nsigma, which we can convert into prob.
@@ -293,7 +294,8 @@ def removeOutliers(self, stars, logger=None):
         :returns: stars, nremoved   A new list of stars with outliers flagged, and how many outliers
                                     were flagged.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         # First figure out the threshold we actually want to use given max_remove.
 
@@ -301,7 +303,7 @@ def removeOutliers(self, stars, logger=None):
         # But note that we apply the threshold to everything eventually.
         use_stars = [star for star in stars if not star.is_flagged and not star.is_reserve]
         nstars = len(use_stars)
-        logger.info("Checking %d stars for outliers", nstars)
+        logger.verbose("Checking %d stars for outliers", nstars)
 
         chisq = np.array([s.fit.chisq for s in use_stars])
         dof = np.array([s.fit.dof for s in use_stars])
@@ -313,17 +315,17 @@ def removeOutliers(self, stars, logger=None):
         thresh = np.array([self._get_thresh(d) for d in dof]) * factor
 
         if np.all(dof == dof[0]):
-            logger.info("dof = %f, thresh = %f * %f = %f",
+            logger.verbose("dof = %f, thresh = %f * %f = %f",
                          dof[0], self._get_thresh(dof[0]), factor, thresh[0])
         else:
             min_dof = np.min(dof)
             max_dof = np.max(dof)
             min_thresh = self._get_thresh(min_dof)
             max_thresh = self._get_thresh(max_dof)
-            logger.info("Minimum dof = %d with thresh = %f * %f = %f",
-                         min_dof, min_thresh, factor, min_thresh*factor)
-            logger.info("Maximum dof = %d with thresh = %f * %f = %f",
-                         max_dof, max_thresh, factor, max_thresh*factor)
+            logger.verbose("Minimum dof = %d with thresh = %f * %f = %f",
+                           min_dof, min_thresh, factor, min_thresh*factor)
+            logger.verbose("Maximum dof = %d with thresh = %f * %f = %f",
+                           max_dof, max_thresh, factor, max_thresh*factor)
 
         nremoved = np.sum(~(chisq <= thresh))  # Write it as not chisq <= thresh in case of nans.
 
@@ -335,7 +337,7 @@ def removeOutliers(self, stars, logger=None):
                      for s in stars]
             return stars, 0
 
-        logger.info("Found %d stars with chisq > thresh", nremoved)
+        logger.verbose("Found %d stars with chisq > thresh", nremoved)
 
         # Update max_remove if necessary
         max_remove = self.max_remove
diff --git a/piff/output.py b/piff/output.py
index 2784e97e..0cedc96e 100644
--- a/piff/output.py
+++ b/piff/output.py
@@ -160,7 +160,8 @@ def write(self, psf, logger=None):
         :param psf:         A PSF instance
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         ensure_dir(self.file_name)
         psf.write(self.file_name, logger=logger)
 
diff --git a/piff/pixelgrid.py b/piff/pixelgrid.py
index 02e21b4d..2bfa52a9 100644
--- a/piff/pixelgrid.py
+++ b/piff/pixelgrid.py
@@ -70,8 +70,8 @@ class specifies a different default.]
 
     def __init__(self, scale, size, interp=None, centered=True, init=None, fit_flux=False,
                  logger=None):
-
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.debug("Building Pixel model with the following parameters:")
         logger.debug("scale = %s",scale)
         logger.debug("size = %s",size)
@@ -124,7 +124,8 @@ def initialize(self, star, logger=None, default_init=None):
 
         :returns: a star instance with the appropriate initial fit values
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         init = self._init if self._init is not None else default_init
         if init is None: init = 'hsm'
         logger.debug("initializing PixelGrid with method %s",init)
@@ -161,7 +162,7 @@ def initialize(self, star, logger=None, default_init=None):
 
                 # Also, make sure fit_flux=True.  Otherwise, this won't work properly.
                 if not self._fit_flux:
-                    logger.info('Setting fit_flux=True, since init=zero')
+                    logger.verbose('Setting fit_flux=True, since init=zero')
                     self._fit_flux = True
                     self.kwargs['fit_flux'] = True
 
@@ -192,8 +193,9 @@ def fit(self, star, logger=None, convert_func=None, draw_method=None):
 
         :returns: a new Star instance with updated fit information
         """
+        from .config import LoggerWrapper
         assert draw_method in (None, 'no_pixel')
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         # Get chisq Taylor expansion for linearized model
         star1 = self.chisq(star, logger=logger, convert_func=convert_func)
 
@@ -229,7 +231,7 @@ def fit(self, star, logger=None, convert_func=None, draw_method=None):
             params_var = np.diagonal(scipy.linalg.inv(star1.fit.A.T.dot(star1.fit.A)))
         except np.linalg.LinAlgError as e:
             # If we get an error, set the variance to "infinity".
-            logger.info("Caught error %s making params_var.  Setting all to 1.e100",e)
+            logger.verbose("Caught error %s making params_var.  Setting all to 1.e100",e)
             params_var = np.ones_like(dparam) * 1.e100
 
         params = star.fit.get_params(self._num)
@@ -258,8 +260,9 @@ def chisq(self, star, logger=None, convert_func=None, draw_method=None):
 
         :returns: a new Star instance with updated fit parameters. (esp. A,b)
         """
+        from .config import LoggerWrapper
         assert draw_method in (None, 'no_pixel')
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         logger.debug('Start chisq function')
         logger.debug('initial params = %s',star.fit.get_params(self._num))
 
diff --git a/piff/polynomial_interp.py b/piff/polynomial_interp.py
index 11b83880..daf8e190 100644
--- a/piff/polynomial_interp.py
+++ b/piff/polynomial_interp.py
@@ -260,8 +260,9 @@ def solve(self, stars, logger=None):
         :param stars:       A list of Star instances to use for the interpolation.
         :param logger:      A logger object for logging debug info. [default: None]
         """
+        from .config import LoggerWrapper
         import scipy.optimize
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
 
         # We will want to index things later, so useful
         # to convert these to numpy arrays and transpose
@@ -275,9 +276,9 @@ def solve(self, stars, logger=None):
         npos = len(positions)
         self._setup_indices(nparam)
 
-        logger.info("Fitting %d parameter vectors using "\
-                    "polynomial type %s with %d positions",
-                    nparam,self.poly_type,npos)
+        logger.verbose("Fitting %d parameter vectors using "\
+                       "polynomial type %s with %d positions",
+                       nparam,self.poly_type,npos)
 
         # This model function adapts our _interpolationModel method
         # into the form that the scipy curve_fit function is expecting.
diff --git a/piff/psf.py b/piff/psf.py
index 303cbfd8..c668ecd2 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -64,7 +64,8 @@ def process(cls, config_psf, logger=None):
 
         :returns: a PSF instance of the appropriate type.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.debug("Parsing PSF based on config dict:")
 
         # Get the class to use for the PSF
@@ -80,7 +81,7 @@ def process(cls, config_psf, logger=None):
         kwargs = psf_cls.parseKwargs(config_psf, logger)
 
         # Build PSF object
-        logger.info("Building %s",psf_type)
+        logger.verbose("Building %s",psf_type)
         psf = psf_cls(**kwargs)
         logger.debug("Done building PSF")
 
@@ -221,7 +222,8 @@ def reflux(self, star, logger=None):
 
         :returns:           New Star instance, with updated flux, center, chisq, dof
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.debug("Reflux for star:")
         logger.debug("    flux = %s",star.fit.flux)
         logger.debug("    center = %s",star.fit.center)
@@ -337,7 +339,7 @@ def reflux_stars(self, stars, logger):
                 if not star.is_flagged:
                     logger.warning("Failed trying to reflux star at %s.  Excluding it.",
                                    star.image_pos)
-                    logger.warning("  -- Caught exception: %s", e)
+                    logger.info("  -- Caught exception: %s", e)
                     nremoved += 1
                     star = star.flag_if(True)
             new_stars.append(star)
@@ -359,7 +361,7 @@ def remove_outliers(self, stars, iteration, logger):
             if nremoved == 0:
                 logger.debug("             No outliers found")
             else:
-                logger.info("             Removed %d outliers", nremoved)
+                logger.verbose("             Removed %d outliers", nremoved)
         else:
             nremoved = 0
         return stars, nremoved
@@ -379,7 +381,8 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
         :param draw_method:     The method to use with the GalSim drawImage command. If not given,
                                 use the default method for the PSF model being fit. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         self.wcs = wcs
         self.pointing = pointing
@@ -393,9 +396,9 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
         for iteration in range(self.max_iter):
 
             nstars = np.sum([not star.is_reserve and not star.is_flagged for star in stars])
-            logger.warning("Iteration %d: Fitting %d stars", iteration+1, nstars)
+            logger.info("Iteration %d: Fitting %d stars", iteration+1, nstars)
             if nreserve != 0:
-                logger.warning("             (%d stars are reserved)", nreserve)
+                logger.info("             (%d stars are reserved)", nreserve)
 
             # Run a single iteration of the fitter.
             # Different PSF types do different things here.
@@ -416,7 +419,7 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
 
             chisq = np.sum([s.fit.chisq for s in stars if not s.is_reserve and not s.is_flagged])
             dof   = np.sum([s.fit.dof for s in stars if not s.is_reserve and not s.is_flagged])
-            logger.warning("             Total chisq = %.2f / %d dof", chisq, dof)
+            logger.info("             Total chisq = %.2f / %d dof", chisq, dof)
 
             # Save these so we can write them to the output file.
             self.chisq = chisq
@@ -506,7 +509,8 @@ def draw(self, x, y, chipnum=None, flux=1.0, center=None, offset=None, stamp_siz
 
         :returns:           A GalSim Image of the PSF
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         chipnum = self._check_chipnum(chipnum)
 
@@ -580,7 +584,8 @@ def get_profile(self, x, y, chipnum=None, flux=1.0, logger=None, **kwargs):
                             method = either 'no_pixel' or 'auto' indicating which method to use
                             when drawing the profile on an image.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         chipnum = self._check_chipnum(chipnum)
 
@@ -721,9 +726,10 @@ def write(self, file_name, logger=None):
         :param logger:      A logger object for logging debug info. [default: None]
         """
         from .writers import FitsWriter
+        from .config import LoggerWrapper
 
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.warning("Writing PSF to file %s",file_name)
+        logger = LoggerWrapper(logger)
+        logger.info("Writing PSF to file %s",file_name)
 
         with FitsWriter.open(file_name) as w:
             self._write(w, 'psf', logger=logger)
@@ -739,14 +745,14 @@ def _write(self, writer, name, logger):
         from . import __version__ as piff_version
         psf_type = self._type_name
         writer.write_struct(name, dict(self.kwargs, type=psf_type, piff_version=piff_version))
-        logger.info("Wrote the basic PSF information to name %s", writer.get_full_name(name))
+        logger.verbose("Wrote the basic PSF information to name %s", writer.get_full_name(name))
         with writer.nested(name) as w:
             if hasattr(self, 'stars'):
                 Star.write(self.stars, w, 'stars')
-                logger.info("Wrote the PSF stars to name %s", w.get_full_name('stars'))
+                logger.verbose("Wrote the PSF stars to name %s", w.get_full_name('stars'))
             if hasattr(self, 'wcs'):
                 w.write_wcs_map('wcs', self.wcs, self.pointing)
-                logger.info("Wrote the PSF WCS to name %s", w.get_full_name('wcs'))
+                logger.verbose("Wrote the PSF WCS to name %s", w.get_full_name('wcs'))
             self._finish_write(w, logger=logger)
 
     @classmethod
@@ -759,9 +765,10 @@ def read(cls, file_name, logger=None):
         :returns: a PSF instance
         """
         from .readers import FitsReader
+        from .config import LoggerWrapper
 
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.warning("Reading PSF from file %s",file_name)
+        logger = LoggerWrapper(logger)
+        logger.info("Reading PSF from file %s",file_name)
 
         with FitsReader.open(file_name) as r:
             logger.debug('opened FITS file')
diff --git a/piff/readers.py b/piff/readers.py
index c44a574b..581c06a8 100644
--- a/piff/readers.py
+++ b/piff/readers.py
@@ -134,8 +134,8 @@ def read_wcs_map(self, name, logger):
             # If this doesn't work, then the file was probably written by py2, not py3
             repr(wcs)
         except Exception:
-            logger.info('Failed to decode wcs with bytes encoding.')
-            logger.info('Retry with encoding="latin1" in case file written with python 2.')
+            logger.verbose('Failed to decode wcs with bytes encoding.')
+            logger.verbose('Retry with encoding="latin1" in case file written with python 2.')
             wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
             wcs = dict(zip(chipnums, wcs_list))
             repr(wcs)
diff --git a/piff/select.py b/piff/select.py
index 487df33f..cd7f5c60 100644
--- a/piff/select.py
+++ b/piff/select.py
@@ -196,10 +196,12 @@ def reserveStars(self, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
+        from .config import LoggerWrapper
+
         if self.reserve_frac == 0:
             return
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
 
         # Apply the reserve separately on each ccd, so they each reserve 20% of their stars
         # (or whatever fraction).  We wouldn't want to accidentally reserve all the stars on
@@ -211,12 +213,12 @@ def reserveStars(self, stars, logger=None):
             # Mark a fraction of the stars as reserve stars
             nreserve = int(self.reserve_frac * len(chip_stars))  # round down
             nreserve_all += nreserve
-            logger.info("Reserve %s of %s (reserve_frac=%s) input stars on chip %s",
-                        nreserve, len(stars), self.reserve_frac, chip_stars[0]['chipnum'])
+            logger.verbose("Reserve %s of %s (reserve_frac=%s) input stars on chip %s",
+                           nreserve, len(stars), self.reserve_frac, chip_stars[0]['chipnum'])
             reserve_list = self.rng.choice(len(chip_stars), nreserve, replace=False)
             for i, star in enumerate(chip_stars):
                 star.data.properties['is_reserve'] = i in reserve_list
-        logger.warning("Reserved %s of %s total stars", nreserve_all, len(stars))
+        logger.info("Reserved %s of %s total stars", nreserve_all, len(stars))
 
     def rejectStars(self, stars, logger=None):
         """Reject some nominal stars that may not be good exemplars of the PSF.
@@ -226,8 +228,9 @@ def rejectStars(self, stars, logger=None):
 
         :returns: The subset of the input list that passed the rejection cuts.
         """
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.info('start rejectStars: %s',len(stars))
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
+        logger.verbose('start rejectStars: %s',len(stars))
 
         if self.max_edge_frac is not None and len(stars) > 0:
             stamp_size = stars[0].image.array.shape[0]
@@ -245,17 +248,17 @@ def rejectStars(self, stars, logger=None):
             if self.max_mask_pixels is not None:
                 n_masked = np.prod(star.weight.array.shape) - np.count_nonzero(star.weight.array)
                 if n_masked >= self.max_mask_pixels:
-                    logger.info("Star at position %f,%f has %i masked pixels, ",
-                                star.image_pos.x, star.image_pos.y, n_masked)
-                    logger.info("Skipping this star.")
+                    logger.verbose("Star at position %f,%f has %i masked pixels, ",
+                                   star.image_pos.x, star.image_pos.y, n_masked)
+                    logger.verbose("Skipping this star.")
                     continue
 
             # Check the snr and limit it if appropriate
             snr = calculateSNR(star.image, star.weight)
             logger.debug("SNR = %f",snr)
             if self.min_snr is not None and snr < self.min_snr:
-                logger.info("Skipping star at position %f,%f with snr=%f.",
-                            star.image_pos.x, star.image_pos.y, snr)
+                logger.verbose("Skipping star at position %f,%f with snr=%f.",
+                               star.image_pos.x, star.image_pos.y, snr)
                 continue
             if self.max_snr > 0 and snr > self.max_snr:
                 factor = (self.max_snr / snr)**2
@@ -273,25 +276,25 @@ def rejectStars(self, stars, logger=None):
                     flux_extra = np.sum(star.image.array[edge_mask])
                     flux_frac = flux_extra / flux
                 except IndexError:
-                    logger.info("Star at position %f,%f overlaps the edge of the image and "+
-                                "max_edge_frac cut is set.",
-                                star.image_pos.x, star.image_pos.y)
-                    logger.info("Skipping this star.")
+                    logger.verbose("Star at position %f,%f overlaps the edge of the image and "+
+                                   "max_edge_frac cut is set.",
+                                   star.image_pos.x, star.image_pos.y)
+                    logger.verbose("Skipping this star.")
                     continue
                 if flux_frac > self.max_edge_frac:
-                    logger.info("Star at position %f,%f fraction of flux near edge of stamp "+
-                                "exceeds cut: %f > %f",
-                                star.image_pos.x, star.image_pos.y,
-                                flux_frac, self.max_edge_frac)
-                    logger.info("Skipping this star.")
+                    logger.verbose("Star at position %f,%f fraction of flux near edge of stamp "+
+                                   "exceeds cut: %f > %f",
+                                   star.image_pos.x, star.image_pos.y,
+                                   flux_frac, self.max_edge_frac)
+                    logger.verbose("Skipping this star.")
                     continue
 
             if self.reject_where is not None:
                 # Use the eval_where function of PropertiesSelect
                 reject = PropertiesSelect.eval_where([star], self.reject_where, logger=logger)
                 if reject:
-                    logger.info("Skipping star at position %f,%f due to reject_where",
-                                star.image_pos.x, star.image_pos.y)
+                    logger.verbose("Skipping star at position %f,%f due to reject_where",
+                                   star.image_pos.x, star.image_pos.y)
                     logger.debug("reject_where string: %s",self.reject_where)
                     logger.debug("star properties = %s",star.data.properties)
                     continue
@@ -338,13 +341,13 @@ def rejectStars(self, stars, logger=None):
                 ratio = np.median(max_pixel[bright] / flux[bright])
                 flux_cut = self.max_pixel_cut / ratio
                 logger.debug("max_pixel/flux ratio = %.4f, flux_cut is %.1f", ratio, flux_cut)
-                logger.info("Rejected %d stars for having a flux > %.1f, "
-                            "which implies max_pixel > ~%s",
-                            np.sum(flux>=flux_cut), flux_cut, self.max_pixel_cut)
+                logger.verbose("Rejected %d stars for having a flux > %.1f, "
+                               "which implies max_pixel > ~%s",
+                               np.sum(flux>=flux_cut), flux_cut, self.max_pixel_cut)
                 good_stars = [s for f,s in zip(flux,good_stars) if f < flux_cut]
 
-        logger.warning("Rejected a total of %d stars out of %s total candidates",
-                       len(stars) - len(good_stars), len(stars))
+        logger.info("Rejected a total of %d stars out of %s total candidates",
+                    len(stars) - len(good_stars), len(stars))
         return good_stars
 
 
@@ -391,13 +394,14 @@ def selectStars(self, objects, logger=None):
 
         :returns: a list of Star instances
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         if self.flag_name is None:
-            logger.info("Using all input objects as stars")
+            logger.verbose("Using all input objects as stars")
             return objects
 
-        logger.info("Selecting stars according to flag %r", self.flag_name)
+        logger.verbose("Selecting stars according to flag %r", self.flag_name)
         try:
             flag_array = np.array([obj[self.flag_name] for obj in objects])
         except KeyError:
@@ -412,7 +416,7 @@ def selectStars(self, objects, logger=None):
             select = InputFiles._flag_select(flag_array, self.skip_flag) == 0
 
         stars = [obj for use, obj in zip(select, objects) if use]
-        logger.info("Seleced %d stars from %d total candidates.", len(stars), len(objects))
+        logger.verbose("Seleced %d stars from %d total candidates.", len(stars), len(objects))
         return stars
 
 class PropertiesSelect(Select):
@@ -446,7 +450,8 @@ def eval_where(cls, objects, where, logger=None):
 
         Used by both PropertiesSelect and the reject_where option.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         # Build appropriate locals and globals for the eval statement.
         gdict = globals().copy()
         # Import some likely packages in case needed.
@@ -461,9 +466,9 @@ def eval_where(cls, objects, where, logger=None):
         try:
             select = eval(where, gdict, ldict)
         except Exception as e:
-            logger.info("Caught exception trying to evaluate where string")
-            logger.info("%r",e)
-            logger.info("Trying slower non-numpy array method")
+            logger.verbose("Caught exception trying to evaluate where string")
+            logger.verbose("%r",e)
+            logger.verbose("Trying slower non-numpy array method")
             select = []
             for obj in objects:
                 ldict = {}
@@ -480,12 +485,13 @@ def selectStars(self, objects, logger=None):
 
         :returns: a list of Star instances
         """
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.info("Selecting stars according to %r", self.where)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
+        logger.verbose("Selecting stars according to %r", self.where)
 
         select = self.eval_where(objects, self.where)
         logger.debug("select = %s",select)
 
         stars = [obj for use, obj in zip(select, objects) if use]
-        logger.info("Seleced %d stars from %d total candidates.", len(stars), len(objects))
+        logger.verbose("Seleced %d stars from %d total candidates.", len(stars), len(objects))
         return stars
diff --git a/piff/simplepsf.py b/piff/simplepsf.py
index 0992aaea..31c083dd 100644
--- a/piff/simplepsf.py
+++ b/piff/simplepsf.py
@@ -141,7 +141,7 @@ def initialize_params(self, stars, logger=None, default_init=None):
         if nremoved == 0:
             logger.debug("No stars removed in initialize step")
         else:
-            logger.info("Removed %d stars in initialize", nremoved)
+            logger.verbose("Removed %d stars in initialize", nremoved)
 
         logger.debug("Initializing interpolator")
         stars = self.interp.initialize(new_stars, logger=logger)
@@ -247,7 +247,8 @@ def _finish_write(self, writer, logger):
         :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         chisq_dict = {
             'chisq' : self.chisq,
             'last_delta_chisq' : self.last_delta_chisq,
diff --git a/piff/singlechip.py b/piff/singlechip.py
index 6bc088fc..1c4c43cc 100644
--- a/piff/singlechip.py
+++ b/piff/singlechip.py
@@ -33,7 +33,7 @@ def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_funcs, dr
     wcs_chip = { chipnum : wcs[chipnum] }
 
     # Run the psf_chip fit function using this stars and wcs (and the same pointing)
-    logger.warning("Building solution for chip %s with %d stars", chipnum, len(stars_chip))
+    logger.info("Building solution for chip %s with %d stars", chipnum, len(stars_chip))
     psf_chip.fit(stars_chip, wcs_chip, pointing, logger=logger, convert_funcs=convert_funcs,
                  draw_method=draw_method)
 
@@ -110,7 +110,8 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
         :param draw_method:     The method to use with the GalSim drawImage command. If not given,
                                 use the default method for the PSF model being fit. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         self.wcs = wcs
         self.pointing = pointing
         self.psf_by_chip = {}
diff --git a/piff/size_mag.py b/piff/size_mag.py
index cf84f886..53af2108 100644
--- a/piff/size_mag.py
+++ b/piff/size_mag.py
@@ -41,10 +41,11 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         # measure everything with hsm
-        logger.info("Measuring Star and Model Shapes")
+        logger.verbose("Measuring Star and Model Shapes")
         _, star_shapes, psf_shapes = self.measureShapes(psf, stars, logger=logger)
 
         if hasattr(psf, 'initial_stars'):
@@ -94,7 +95,8 @@ def plot(self, logger=None, **kwargs):
         :returns: fig, ax
         """
         from matplotlib.figure import Figure
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         fig = Figure(figsize=(8,6))
         ax = fig.add_subplot(1,1,1)
 
@@ -203,9 +205,10 @@ def selectStars(self, objects, logger=None):
 
         :returns: a list of Star instances
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
-        logger.warning("Selecting small/bright objects as stars")
+        logger.info("Selecting small/bright objects as stars")
 
         logger.debug("Initial count = %s", len(objects))
 
@@ -273,7 +276,7 @@ def selectStars(self, objects, logger=None):
         imin = np.argmin(delta_T)
         logger.debug("imin = %s", imin)
         star_logT = bright_small_logT[imin:imin+half_len+1]
-        logger.info("Initial bright/small selection includes %d objects",half_len+1)
+        logger.verbose("Initial bright/small selection includes %d objects",half_len+1)
 
         # Expand this to include all stars that are within twice the interquarile range of
         # these candidate stars.  Keep doing so until we converge on a good set of stars.
@@ -286,7 +289,7 @@ def selectStars(self, objects, logger=None):
             logger.debug("Iteration %d",it)
             logger.debug("Sizes of candidate stars = %s", np.exp(star_logT))
             med = np.median(star_logT)
-            logger.info("Median size = %s", np.exp(med))
+            logger.verbose("Median size = %s", np.exp(med))
             q25, q75 = np.percentile(star_logT, [25,75])
             iqr = q75 - q25
             logger.debug("Range of star logT size = %s, %s", np.min(star_logT), np.max(star_logT))
@@ -300,10 +303,10 @@ def selectStars(self, objects, logger=None):
             if np.array_equal(select, old_select):
                 break
             old_select = select
-            logger.info("Expand this to include %d selected stars",new_count)
+            logger.verbose("Expand this to include %d selected stars",new_count)
             star_logT = logT[select]
         else:
-            logger.info("Max iter = 10 reached.  Stop updating based on median/IQR.")
+            logger.verbose("Max iter = 10 reached.  Stop updating based on median/IQR.")
 
         # Get the initial indexes of these objects
         select_index = orig_index[select]
@@ -311,7 +314,7 @@ def selectStars(self, objects, logger=None):
         stars = [objects[i] for i in select_index]
         logger.debug("sizes of stars = %s",[2*s.hsm[3]**2 for s in stars])
         logger.debug("fluxs of stars = %s",[s.hsm[0] for s in stars])
-        logger.warning("Bright/small selection found %d likely stars",len(stars))
+        logger.info("Bright/small selection found %d likely stars",len(stars))
 
         return stars
 
@@ -371,6 +374,8 @@ class SizeMagSelect(Select):
     _type_name = 'SizeMag'
 
     def __init__(self, config, logger=None):
+        from .config import LoggerWrapper
+
         super(SizeMagSelect, self).__init__(config, logger)
 
         opt = {
@@ -382,7 +387,7 @@ def __init__(self, config, logger=None):
 
         # This used to be (badly) named purity.
         # For backwards compatibility, switch it if someone is using the old name.
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         if 'purity' in config:
             logger.error("WARNING: The parameter name 'purity' should now be called 'impurity'.")
             config['impurity'] = config.pop('purity')
@@ -405,9 +410,10 @@ def selectStars(self, objects, logger=None):
 
         :returns: a list of Star instances
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
-        logger.warning("Selecting stars according to locus in size-magnitude diagram")
+        logger.info("Selecting stars according to locus in size-magnitude diagram")
 
         stars = Select.process(self.initial_select, objects, logger=logger, select_only=True)
 
@@ -542,15 +548,15 @@ def selectStars(self, objects, logger=None):
             logT_star = logT_obj[select]
             u_star = u_obj[select]
             v_star = v_obj[select]
-            logger.info("SizeMag iteration %d => N stars = %d", i_iter, len(logf_star))
-            logger.info("Mean logT of stars = %.3f, std = %.3f", np.mean(logT_star), np.std(logT_star))
+            logger.verbose("SizeMag iteration %d => N stars = %d", i_iter, len(logf_star))
+            logger.verbose("Mean logT of stars = %.3f, std = %.3f", np.mean(logT_star), np.std(logT_star))
 
         select_index = orig_index[select]
         logger.debug("select_index = %s",select_index)
         stars = [objects[i] for i in select_index]
         logger.debug("sizes of stars = %s",[2*s.hsm[3]**2 for s in stars])
         logger.debug("fluxs of stars = %s",[s.hsm[0] for s in stars])
-        logger.warning("SizeMag selection found %d likely stars",len(stars))
+        logger.info("SizeMag selection found %d likely stars",len(stars))
 
         return stars
 
diff --git a/piff/star.py b/piff/star.py
index 8ae5b2b5..2c6cc715 100644
--- a/piff/star.py
+++ b/piff/star.py
@@ -560,13 +560,14 @@ def load_images(stars, file_name, pointing=None,
         :returns: a new list of Stars with the images information loaded.
         """
         from .input import InputFiles
+        from .config import LoggerWrapper
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
 
         logger.error("WARNING: The Star.load_images function is deprecated."
                      "Use InputFiles.load_images instead.")
 
-        logger.info("Loading image information from file %s",file_name)
+        logger.verbose("Loading image information from file %s",file_name)
         config = {
             'image_file_name': file_name,
             'cat_file_name': None,  # We don't need this, but it needs to be present.
diff --git a/piff/star_stats.py b/piff/star_stats.py
index 3a4f2f5b..12dd5be4 100644
--- a/piff/star_stats.py
+++ b/piff/star_stats.py
@@ -67,7 +67,8 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         # Determine which stars to plot
         possible_indices = []
         if self.include_reserve:
@@ -91,7 +92,7 @@ def compute(self, psf, stars, logger=None):
         else:
             calculate_indices = self.indices
 
-        logger.info("Making {0} model stars".format(len(calculate_indices)))
+        logger.verbose("Making {0} model stars".format(len(calculate_indices)))
         calculated_stars = []
         calculated_models = []
         for i, star in enumerate(stars):
@@ -107,7 +108,7 @@ def compute(self, psf, stars, logger=None):
                 calculated_models.append(None)
 
         # if including the average image, put that first.
-        logger.info("Making average star and model")
+        logger.verbose("Making average star and model")
         if self.include_ave:
             ave_star_image = np.mean([s.image.array for s in calculated_stars if s is not None],
                                      axis=0)
@@ -136,10 +137,12 @@ def plot(self, logger=None, **kwargs):
         :returns: fig, ax
         """
         from matplotlib.figure import Figure
+        from .config import LoggerWrapper
+
         if not hasattr(self, 'indices'):
             raise RuntimeError("Must call compute before calling plot or write")
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
 
         # 6 x nplot/2 images, with each image (3.5 x 3)
         nplot = len(self.indices)
@@ -147,7 +150,7 @@ def plot(self, logger=None, **kwargs):
         fig = Figure(figsize = (21,3*nrows))
         axs = fig.subplots(ncols=6, nrows=nrows, squeeze=False)
 
-        logger.info("Creating %d Star plots", self.nplot)
+        logger.verbose("Creating %d Star plots", self.nplot)
 
         for i in range(nplot):
             star = self.stars[i]
@@ -194,6 +197,7 @@ class StarStatsDepr(StarStats):
     _type_name = 'Star'
 
     def __init__(self, *args, logger=None, **kwargs):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         logger.error("WARNING: The name Star is deprecated. Use StarImages instead.")
         super().__init__(*args, logger=logger, **kwargs)
diff --git a/piff/stats.py b/piff/stats.py
index 82dd9f5c..88c5948d 100644
--- a/piff/stats.py
+++ b/piff/stats.py
@@ -143,9 +143,10 @@ def write(self, file_name=None, logger=None, **kwargs):
         # cf. http://www.dalkescientific.com/writings/diary/archive/2005/04/23/matplotlib_without_gui.html
         import matplotlib
         from matplotlib.backends.backend_agg import FigureCanvasAgg
+        from .config import LoggerWrapper
 
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.info("Creating plot for %s", self._type_name)
+        logger = LoggerWrapper(logger)
+        logger.verbose("Creating plot for %s", self._type_name)
         fig, ax = self.plot(logger=logger, **kwargs)
 
         if file_name is None:
@@ -153,7 +154,7 @@ def write(self, file_name=None, logger=None, **kwargs):
         if file_name is None:
             raise ValueError("No file_name specified for %s"%self._type_name)
 
-        logger.warning("Writing %s plot to file %s",self._type_name, file_name)
+        logger.info("Writing %s plot to file %s",self._type_name, file_name)
 
         canvas = FigureCanvasAgg(fig)
         # Do this after we've set the canvas to use Agg to avoid warning.
@@ -183,7 +184,9 @@ def measureShapes(self, psf, stars, model_properties=None, fourth_order=False,
                             models of stars (T, g1, g2)
         """
         from .util import calculate_moments
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+
+        logger = LoggerWrapper(logger)
         # measure moments with Gaussian on image
         logger.debug("Measuring shapes of real stars")
 
@@ -348,9 +351,10 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         # get the shapes
-        logger.warning("Calculating shape histograms for %d stars",len(stars))
+        logger.info("Calculating shape histograms for %d stars",len(stars))
         positions, shapes_data, shapes_model = self.measureShapes(
                 psf, stars, model_properties=self.model_properties, logger=logger)
 
@@ -383,7 +387,8 @@ def plot(self, logger=None, **kwargs):
 
         :returns: fig, ax
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         from matplotlib.figure import Figure
         fig = Figure(figsize = (15,10))
         # In matplotlib 2.0, this will be
@@ -415,7 +420,7 @@ def plot(self, logger=None, **kwargs):
         nbins = self.nbins
         if nbins is None:
             nbins = int(np.sqrt(len(self.T))+1)
-            logger.info("nstars = %d, using %d bins for Shape Histograms",len(self.T),nbins)
+            logger.verbose("nstars = %d, using %d bins for Shape Histograms",len(self.T),nbins)
 
         lower = dict(method='lower')
         higher = dict(method='higher')
@@ -423,21 +428,21 @@ def plot(self, logger=None, **kwargs):
         # axs[0,0] = size distributions
         ax = axs[0, 0]
         all_T = np.concatenate([self.T_model, self.T])
-        logger.info("nbins = %s",nbins)
-        logger.info("cut_frac = %s",self.cut_frac)
+        logger.verbose("nbins = %s",nbins)
+        logger.verbose("cut_frac = %s",self.cut_frac)
         rng = (np.quantile(all_T, self.cut_frac, **lower),
                np.quantile(all_T, 1-self.cut_frac, **higher))
-        logger.info("T_d: Full range = (%f, %f)",np.min(self.T),np.max(self.T))
-        logger.info("T_m: Full range = (%f, %f)",np.min(self.T_model),np.max(self.T_model))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("T_d: Full range = (%f, %f)",np.min(self.T),np.max(self.T))
+        logger.verbose("T_m: Full range = (%f, %f)",np.min(self.T_model),np.max(self.T_model))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist([self.T, self.T_model], bins=nbins, range=rng, label=['data', 'model'], **kwargs)
         ax.legend(loc='upper right')
         # axs[0,1] = size difference
         ax = axs[1, 0]
         rng = (np.quantile(self.dT, self.cut_frac, **lower),
                np.quantile(self.dT, 1-self.cut_frac, **higher))
-        logger.info("dT: Full range = (%f, %f)",np.min(self.dT),np.max(self.dT))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("dT: Full range = (%f, %f)",np.min(self.dT),np.max(self.dT))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist(self.dT, bins=nbins, range=rng, **kwargs)
 
         # axs[1,0] = g1 distribution
@@ -445,17 +450,17 @@ def plot(self, logger=None, **kwargs):
         all_g1 = np.concatenate([self.g1_model, self.g1])
         rng = (np.quantile(all_g1, self.cut_frac, **lower),
                np.quantile(all_g1, 1-self.cut_frac, **higher))
-        logger.info("g1_d: Full range = (%f, %f)",np.min(self.g1),np.max(self.g1))
-        logger.info("g1_m: Full range = (%f, %f)",np.min(self.g1_model),np.max(self.g1_model))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("g1_d: Full range = (%f, %f)",np.min(self.g1),np.max(self.g1))
+        logger.verbose("g1_m: Full range = (%f, %f)",np.min(self.g1_model),np.max(self.g1_model))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist([self.g1, self.g1_model], bins=nbins, range=rng, label=['data', 'model'], **kwargs)
         ax.legend(loc='upper right')
         # axs[1,0] = g1 difference
         ax = axs[1, 1]
         rng = (np.quantile(self.dg1, self.cut_frac, **lower),
                np.quantile(self.dg1, 1-self.cut_frac, **higher))
-        logger.info("dg1: Full range = (%f, %f)",np.min(self.dg1),np.max(self.dg1))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("dg1: Full range = (%f, %f)",np.min(self.dg1),np.max(self.dg1))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist(self.dg1, bins=nbins, range=rng, **kwargs)
 
         # axs[2,0] = g2 distribution
@@ -463,17 +468,17 @@ def plot(self, logger=None, **kwargs):
         all_g2 = np.concatenate([self.g2_model, self.g2])
         rng = (np.quantile(all_g2, self.cut_frac, **lower),
                np.quantile(all_g2, 1-self.cut_frac, **higher))
-        logger.info("g2_d: Full range = (%f, %f)",np.min(self.g2),np.max(self.g2))
-        logger.info("g2_m: Full range = (%f, %f)",np.min(self.g2_model),np.max(self.g2_model))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("g2_d: Full range = (%f, %f)",np.min(self.g2),np.max(self.g2))
+        logger.verbose("g2_m: Full range = (%f, %f)",np.min(self.g2_model),np.max(self.g2_model))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist([self.g2, self.g2_model], bins=nbins, range=rng, label=['data', 'model'], **kwargs)
         ax.legend(loc='upper right')
         # axs[2,0] = g2 difference
         ax = axs[1, 2]
         rng = (np.quantile(self.dg2, self.cut_frac, **lower),
                np.quantile(self.dg2, 1-self.cut_frac, **higher))
-        logger.info("dg2: Full range = (%f, %f)",np.min(self.dg2),np.max(self.dg2))
-        logger.info("Display range = (%f, %f)",rng[0],rng[1])
+        logger.verbose("dg2: Full range = (%f, %f)",np.min(self.dg2),np.max(self.dg2))
+        logger.verbose("Display range = (%f, %f)",rng[0],rng[1])
         ax.hist(self.dg2, bins=nbins, range=rng, **kwargs)
 
         return fig, ax
@@ -539,12 +544,14 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
+        from .config import LoggerWrapper
         import treecorr
+
         treecorr.set_max_omp_threads(1)
 
-        logger = galsim.config.LoggerWrapper(logger)
+        logger = LoggerWrapper(logger)
         # get the shapes
-        logger.warning("Calculating rho statistics for %d stars",len(stars))
+        logger.info("Calculating rho statistics for %d stars",len(stars))
         positions, shapes_data, shapes_model = self.measureShapes(
                 psf, stars, model_properties=self.model_properties, logger=logger)
 
@@ -568,7 +575,7 @@ def compute(self, psf, stars, logger=None):
         dg2 = g2 - shapes_model[mask, 5]
 
         # make the treecorr catalogs
-        logger.info("Creating Treecorr Catalogs")
+        logger.verbose("Creating Treecorr Catalogs")
 
         cat_g = treecorr.Catalog(x=u, y=v, x_units='arcsec', y_units='arcsec',
                                  g1=g1, g2=g2)
@@ -578,7 +585,7 @@ def compute(self, psf, stars, logger=None):
                                     g1=g1 * dT / T, g2=g2 * dT / T)
 
         # setup and run the correlations
-        logger.info("Processing rho PSF statistics")
+        logger.verbose("Processing rho PSF statistics")
 
         # save the rho objects
         self.rho1 = treecorr.GGCorrelation(self.tckwargs)
@@ -814,9 +821,10 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         # get the shapes
-        logger.warning("Calculating shape histograms for %d stars",len(stars))
+        logger.info("Calculating shape histograms for %d stars",len(stars))
         positions, shapes_data, shapes_model = self.measureShapes(
                 psf, stars, model_properties=self.model_properties,
                 fourth_order=self.fourth_order, raw_moments=self.raw_moments,
@@ -911,8 +919,10 @@ def write(self, file_name=None, logger=None):
         :param logger:      A logger object for logging debug info. [default: None]
         """
         from . import __version__ as piff_version
-        logger = galsim.config.LoggerWrapper(logger)
         import fitsio
+        from .config import LoggerWrapper
+
+        logger = LoggerWrapper(logger)
         if file_name is None:
             file_name = self.file_name
         if file_name is None:
@@ -920,7 +930,7 @@ def write(self, file_name=None, logger=None):
         if not hasattr(self, 'cols'):
             raise RuntimeError("Must call compute before calling write")
 
-        logger.warning("Writing HSM catalog to file %s",file_name)
+        logger.info("Writing HSM catalog to file %s",file_name)
 
         data = np.array(list(zip(*self.cols)), dtype=self.dtypes)
         header = {'piff_version': piff_version}
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 59eb4a5f..b64b0d00 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -176,7 +176,7 @@ def single_iteration(self, stars, logger, convert_funcs, draw_method):
 
         # Fit each component in order
         for k, comp in enumerate(self.components):
-            logger.info("Starting work on component %d (%s)", k, comp._type_name)
+            logger.verbose("Starting work on component %d (%s)", k, comp._type_name)
             # Update stars to subtract current fit from other components.
             modified_stars = []
             for i, star in enumerate(stars):
@@ -293,7 +293,8 @@ def _finish_write(self, writer, logger):
         :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         chisq_dict = {
             'chisq' : self.chisq,
             'last_delta_chisq' : self.last_delta_chisq,
diff --git a/piff/twod_stats.py b/piff/twod_stats.py
index 27845841..28854e0a 100644
--- a/piff/twod_stats.py
+++ b/piff/twod_stats.py
@@ -67,11 +67,12 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         self.twodhists = {}
 
         # get the shapes
-        logger.info("Measuring Star and Model Shapes")
+        logger.verbose("Measuring Star and Model Shapes")
         positions, shapes_truth, shapes_model = self.measureShapes(
             psf, stars, model_properties=self.model_properties, logger=logger)
 
@@ -79,7 +80,7 @@ def compute(self, psf, stars, logger=None):
         flag_truth = shapes_truth[:, 6]
         flag_model = shapes_model[:, 6]
         mask = (flag_truth == 0) & (flag_model == 0)
-        logger.info("%d/%d measurements were successful",np.sum(mask),len(mask))
+        logger.verbose("%d/%d measurements were successful",np.sum(mask),len(mask))
         if np.sum(mask) == 0:
             logger.warning("All stars had hsm errors.  TwoDHist plot will be empty.")
             self.skip = True
@@ -99,7 +100,7 @@ def compute(self, psf, stars, logger=None):
         dg2 = g2 - g2_model
 
         # compute the indices
-        logger.info("Computing TwoDHist indices")
+        logger.verbose("Computing TwoDHist indices")
 
         self.bins_u = np.linspace(np.min(u), np.max(u), num=self.nbins_u + 1)
         self.bins_v = np.linspace(np.min(v), np.max(v), num=self.nbins_v + 1)
@@ -116,7 +117,7 @@ def compute(self, psf, stars, logger=None):
         unique_indx = np.vstack([tuple(row) for row in np.vstack((indx_u, indx_v)).T])
 
         # compute the arrays
-        logger.info("Computing TwoDHist arrays")
+        logger.verbose("Computing TwoDHist arrays")
 
         # throw in coordinates for good measure
         self.twodhists['u'] = self._array_to_2dhist(u, indx_u, indx_v, unique_indx)
@@ -158,7 +159,9 @@ def plot(self, logger=None, **kwargs):
         :returns: fig, ax
         """
         from matplotlib.figure import Figure
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+
+        logger = LoggerWrapper(logger)
         fig = Figure(figsize=(12,9))
         # In matplotlib 2.0, this will be
         # axs = fig.subplots(ncols=3, nrows=3)
@@ -199,7 +202,7 @@ def plot(self, logger=None, **kwargs):
             raise RuntimeError("Must call compute before calling plot or write")
 
         # make the colormaps
-        logger.info("Creating TwoDHist colormaps")
+        logger.verbose("Creating TwoDHist colormaps")
         # T and T_model share colorbar
         vmin__T = np.min([self.twodhists['T'], self.twodhists['T_model']])
         vmax__T = np.max([self.twodhists['T'], self.twodhists['T_model']])
@@ -221,7 +224,7 @@ def plot(self, logger=None, **kwargs):
         cmap__dg = self._shift_cmap(vmin__dg, vmax__dg)
 
         # make the plots
-        logger.info("Creating TwoDHist plots")
+        logger.verbose("Creating TwoDHist plots")
         ax = axs[0, 0]
         IM = ax.pcolormesh(self.bins_u, self.bins_v, self.twodhists['T'], cmap=cmap__T,
                        vmin=vmin__T, vmax=vmax__T)
@@ -442,11 +445,12 @@ def compute(self, psf, stars, logger=None):
         :param stars:       A list of Star instances.
         :param logger:      A logger object for logging debug info. [default: None]
         """
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         self.twodhists = {}
 
         # get the shapes
-        logger.info("Measuring Star and Model Shapes")
+        logger.verbose("Measuring Star and Model Shapes")
         positions, shapes_truth, shapes_model = self.measureShapes(
             psf, stars, model_properties=self.model_properties, logger=logger)
 
@@ -454,7 +458,7 @@ def compute(self, psf, stars, logger=None):
         flag_truth = shapes_truth[:, 6]
         flag_model = shapes_model[:, 6]
         mask = (flag_truth == 0) & (flag_model == 0)
-        logger.info("%d/%d measurements were successful",np.sum(mask),len(mask))
+        logger.verbose("%d/%d measurements were successful",np.sum(mask),len(mask))
         if np.sum(mask) == 0:
             logger.warning("All stars had hsm errors.  Whisker plot will be empty.")
             self.skip = True
@@ -487,7 +491,7 @@ def compute(self, psf, stars, logger=None):
         dw2 = dmag_w * np.sin(dphi)
 
         # compute the indices
-        logger.info("Computing Whisker indices")
+        logger.verbose("Computing Whisker indices")
 
         self.bins_u = np.linspace(np.min(u), np.max(u), num=self.nbins_u + 1)
         self.bins_v = np.linspace(np.min(v), np.max(v), num=self.nbins_v + 1)
@@ -504,7 +508,7 @@ def compute(self, psf, stars, logger=None):
         unique_indx = np.vstack([tuple(row) for row in np.vstack((indx_u, indx_v)).T])
 
         # compute the arrays
-        logger.info("Computing Whisker arrays")
+        logger.verbose("Computing Whisker arrays")
 
         self.twodhists['u'] = self._array_to_2dhist(u, indx_u, indx_v, unique_indx)
         self.twodhists['v'] = self._array_to_2dhist(v, indx_u, indx_v, unique_indx)
@@ -537,7 +541,9 @@ def plot(self, logger=None, **kwargs):
         :returns: fig, ax
         """
         from matplotlib.figure import Figure
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+
+        logger = LoggerWrapper(logger)
         fig = Figure(figsize=(7.5,4))
         ax0, ax1 = fig.subplots(1, 2, sharey=True, subplot_kw={'aspect' : 'equal'})
 
@@ -556,7 +562,7 @@ def plot(self, logger=None, **kwargs):
             raise RuntimeError("Must call compute before calling plot or write")
 
         # make the plots
-        logger.info("Creating Whisker plots")
+        logger.verbose("Creating Whisker plots")
 
         # configure to taste
         # bigger scale = smaller whiskers
diff --git a/piff/util.py b/piff/util.py
index c2c9ad62..84421efd 100644
--- a/piff/util.py
+++ b/piff/util.py
@@ -22,6 +22,8 @@
 import traceback
 import galsim
 
+from .config import PiffLogger
+
 def ensure_dir(target):
     """Ensure that the directory for a target output file exists.
 
@@ -142,7 +144,7 @@ def _run_multi_helper(func, i, args, kwargs, log_level): # pragma: no cover
 
     # In multiprocessing, we cannot pass in the logger, so log to a string and then
     # return that back at the end to be logged by the parent process.
-    logger = logging.getLogger('logtostring_%d'%i)
+    logger = PiffLogger(logging.getLogger('logtostring_%d'%i))
     buf = StringIO()
     handler = logging.StreamHandler(buf)
     logger.addHandler(handler)
@@ -155,7 +157,7 @@ def _run_multi_helper(func, i, args, kwargs, log_level): # pragma: no cover
         # is to catch it and return it.  We can deal with it somehow on the other end.
         # Also add more details here with verbose>=2 to help with debugging.
         tr = traceback.format_exc()
-        logger.info("Caught exception:\n%s",tr)
+        logger.verbose("Caught exception:\n%s",tr)
         out = e
 
     handler.flush()
@@ -192,7 +194,7 @@ def run_multi(func, nproc, raise_except, args, logger, kwargs=None):
 
     def log_output(result): # pragma: no cover (It is covered, but in an async process.)
         i, out, log = result
-        logger.info(log)
+        logger.verbose(log)
         if isinstance(out, Exception):
             logger.warning("Caught exception in multiprocessing job: %r",out)
             err_list[i] = out
diff --git a/piff/wavefront.py b/piff/wavefront.py
index c9a51164..ec306074 100644
--- a/piff/wavefront.py
+++ b/piff/wavefront.py
@@ -1,8 +1,23 @@
+# Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at
+# https://github.com/rmjarvis/Piff  All rights reserved.
+#
+# Piff is free software: Redistribution and use in source and binary forms
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+
 import numpy as np
 import galsim
-from .star import Star
 from galsim import LookupTable2D
 
+from .star import Star
+
 def convert_zernikes_des(a_fp):
     """ This method converts an array of Noll Zernike coefficients from
     Focal Plane coordinates to a set suitable for u,v Sky coordinates.
@@ -76,7 +91,8 @@ class Wavefront(object):
     :param logger:             A logger object for logging debug info. [default: None]
     """
     def __init__(self,wavefront_kwargs,logger=None):
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
         self.maxnZ = 37            # hardcoded maximum Zernike index, Noll parameterization
         zformatstr = "z%d"         # hardcoded format string for Zernike coefficients
         self.interp_objects = {}   # store all interpolation objects
@@ -177,8 +193,8 @@ def fillWavefront(self,star_list,wavelength=-1.0,logger=None,addtostars=True):
         :param addtostars:      A boolean flag to return stars with wavefront added to star.data or
                                 return array of wavefronts. [default: True]
         """
-
-        logger = galsim.config.LoggerWrapper(logger)
+        from .config import LoggerWrapper
+        logger = LoggerWrapper(logger)
 
         nstars = len(star_list)
         wf_arr = np.zeros((nstars,self.maxnZ+1))
@@ -196,7 +212,7 @@ def fillWavefront(self,star_list,wavelength=-1.0,logger=None,addtostars=True):
         # loop over wavefront sources
         for isource in range(self.nsources):
 
-            logger.info("Filling Zernike coefficients for source %d" % (isource))
+            logger.verbose("Filling Zernike coefficients for source %d" % (isource))
             xkey,ykey = self.starxykeys[isource]
 
             # get x,y from Star.data, or use x_fp,y_fp above
@@ -231,7 +247,7 @@ def fillWavefront(self,star_list,wavelength=-1.0,logger=None,addtostars=True):
 
         # convert Zernike coeff for DES
         if self.survey == 'des':
-            logger.info("Converting Zernike coefficients for Survey %s" % (self.survey))
+            logger.verbose("Converting Zernike coefficients for Survey %s" % (self.survey))
             wf_arr_final = np.zeros_like(wf_arr)
             for i in range(wf_arr_final.shape[0]):
                 wf_arr_final[i,:] = convert_zernikes_des(wf_arr[i,:])
diff --git a/tests/piff_test_helper.py b/tests/piff_test_helper.py
index 555032a8..612356eb 100644
--- a/tests/piff_test_helper.py
+++ b/tests/piff_test_helper.py
@@ -14,6 +14,7 @@
 
 
 import logging
+from piff.config import VERBOSE, PiffLogger
 
 # Some helper functions that mutliple test files might want to use
 
@@ -98,11 +99,11 @@ class CaptureLog(object):
 
     """
     def __init__(self, level=3):
-        logging_levels = { 0: logging.CRITICAL,
-                           1: logging.WARNING,
-                           2: logging.INFO,
+        logging_levels = { 0: logging.ERROR,
+                           1: logging.INFO,
+                           2: VERBOSE,
                            3: logging.DEBUG }
-        self.logger = logging.getLogger('CaptureLog')
+        self.logger = PiffLogger(logging.getLogger('CaptureLog'))
         self.logger.setLevel(logging_levels[level])
         try:
             from StringIO import StringIO