Repository: rmjarvis/Piff (language: python)

Issue:
Serialization for embedding in larger (Rubin) files
Rubin's been storing its Piff-generated PSFs via pickle thus far, and we all know that's a terrible idea.  We've had a few vague thoughts about what we might do instead, but some of those are finally firming up and it's past time we did something.

We'd like to be able to store arbitrary Piff PSFs in a backwards-compatible way (future Piff versions can read files written by older versions) inside larger files.  Most of those files will be FITS, but we don't want to assume they all will be.  And while we could have Rubin-side code that introspects Piff objects of the types we actually use and handles the serialization, it be better to do this directly in Piff and more naturally get support for all types of models, interpolants, etc.

Here are some possibilities:

1. One option is Piff's existing FITS serialization, which we probably should have tried to use instead of pickle from the beginning.  In many cases we could give Piff a `fitsio.FITS` object and let it append some additional HDUs.  We'd need to make sure it didn't assume it was writing a full file, but I assume that's either already the case or not hard. And we'd probably need it to report something about what it did write, so that we can have pointers to those HDUs elsewhere in the file.  The big downside of this approach is that it does require FITS, and that means in other contexts (and maybe in some FITS contexts where appending HDUs might be problematic) we'd have to write to FITS and then consider that a binary blob to be stuffed into some other file.

2. If we're going to write a binary blob sometimes (and consider Piff-appended HDUs to be totally opaque in terms of how we document our files for science users in other places), maybe we should just always store Piff PSFs as binary blobs, but make it a Piff method that does that.  The first version could even delegate to the FITS serialization, then run that through gzip or something and prepend a serialization version number for backwards compatibility, and if we wanted to rewrite that without delegating to FITS later it'd be a relatively transparent change (as long as we could still read the FITS version in the future).

3. I've been putting together a new pure-Python library [ShoeFits](https://github.com/lsst-dm/shoefits) for serialization on top of [pydantic](https://docs.pydantic.dev/latest/).  On its own Pydantic lets you define dataclass-like structs and get them mapped to and from JSON automatically (with nice support for JSON schemas and hence format documentation), and what ShoeFits adds is the ability to pull numpy arrays out of those structs and save them elsewhere, with pointers replacing them in the JSON tree that are automatically deferenced on read.  It's not on PyPI yet, but it will be (pending formal Rubin DM adoption; this is still something we're discussing on an [RFC](https://rubinobs.atlassian.net/browse/RFC-1030)), so I'm hoping adding it as an (optional) dependency of Piff might be possible.  I think this is the best approach to take if you want the Piff representation in Rubin files to be human-readable or self-describing (this is not something DM requires, though it might be nice), and it's probably easier to implement and maintain than a bespoke non-FITS binary blob format.

Rubin DM can definitely contribute effort to any of these (probably mostly @PFLeget, maybe some from @jmeyers314 and me), and we're open to other options too if you have ideas.


Discussion / hints:
Case 1 is already enabled.  The UI is
```
psf._write(fits, extname, logger)
```
Indeed the normal front end writer, which creates the FITS file, is just implemented as
```
        with fitsio.FITS(file_name,'rw',clobber=True) as f:
            self._write(f, 'psf', logger)
```
If 'psf' is too generic, you can change that to something longer, like 'piff_psf' or whatever when you call the underscore version.
Relevant links:
`_write`:
https://github.com/rmjarvis/Piff/blob/releases/1.4/piff/psf.py#L731

The corresponding `_read` on an existing open FITS file:
https://github.com/rmjarvis/Piff/blob/releases/1.4/piff/psf.py#L771
Case 3 sounds reasonable.  I'd definitely consider any PRs towards that end.  

We can also do something like what I do in TreeCorr where I have several writers and readers, so the first step is to create one of these according to the desired file format.  And then the actual I/O all works through methods of these classes.

cf. https://github.com/rmjarvis/TreeCorr/blob/releases/5.0/treecorr/writer.py
https://github.com/rmjarvis/TreeCorr/blob/releases/5.0/treecorr/reader.py

The front end for this looks like (https://github.com/rmjarvis/TreeCorr/blob/releases/5.0/treecorr/zzcorrelation.py#L191)
```
        with make_writer(file_name, precision, file_type, self.logger) as writer:
            self._write(writer, None, write_patch_results, write_cov=write_cov)
```
Then the underscore write functions work with methods of the `writer` object.  cf. https://github.com/rmjarvis/TreeCorr/blob/releases/5.0/treecorr/corr2base.py#L1700

So we could maintain a FITS writer/reader, and something else that Rubin DM could maintain.  We could even get fancy and let you maintain that external to the Piff repo if that makes more sense and use a registration step to let you register your writer class to be used by Piff.
It'd be really easy to map case 3 to `reader` / `writer` interfaces like the TreeCorr ones, and it does look like that would help a new serialization format share code with the existing FITS one.  And while I don't think Rubin DM needs to have something external to the Piff repo, that might help sidestep dependency ordering problems in the future.

So I think I'll start by trying to put together a PR that refactors the FITS I/O into `reader` / `writer` classes like the TreeCorr ones, while making sure the methods on those are compatible with what I'd want to implement them with ShoeFits.

If that all goes well, I'll put together a ShoeFits implementation in a Rubin DM repo, and we can discuss whether to upstream it to Piff later.

Fix patch:
diff --git a/piff/__init__.py b/piff/__init__.py
index 685b063e..1c2c1f90 100644
--- a/piff/__init__.py
+++ b/piff/__init__.py
@@ -120,3 +120,5 @@
 from . import des
 from . import wavefront
 from . import meta_data
+from . import readers
+from . import writers
diff --git a/piff/basis_interp.py b/piff/basis_interp.py
index bd02f8c1..81581e94 100644
--- a/piff/basis_interp.py
+++ b/piff/basis_interp.py
@@ -434,11 +434,10 @@ def constant(self, value=1.):
         out[0] = value  # The constant term is always first.
         return out
 
-    def _finish_write(self, fits, extname):
-        """Write the solution to a FITS binary table.
+    def _finish_write(self, writer):
+        """Write the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension.
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         if self.q is None:
             raise RuntimeError("Solution not set yet.  Cannot write this BasisPolynomial.")
@@ -446,14 +445,14 @@ def _finish_write(self, fits, extname):
         dtypes = [ ('q', float, self.q.shape) ]
         data = np.zeros(1, dtype=dtypes)
         data['q'] = self.q
-        fits.write_table(data, extname=extname + '_solution')
+        writer.write_table('solution', data)
 
-    def _finish_read(self, fits, extname):
-        """Read the solution from a FITS binary table.
+    def _finish_read(self, reader):
+        """Read the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The name of the extension with the interpolator information.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
-        data = fits[extname + '_solution'].read()
+        data = reader.read_table('solution')
+        assert data is not None
         self.q = data['q'][0]
 
diff --git a/piff/convolvepsf.py b/piff/convolvepsf.py
index cabf1503..bc01f8fd 100644
--- a/piff/convolvepsf.py
+++ b/piff/convolvepsf.py
@@ -20,7 +20,6 @@
 import galsim
 
 from .psf import PSF
-from .util import write_kwargs, read_kwargs
 from .star import Star, StarFit
 from .outliers import Outliers
 
@@ -271,11 +270,10 @@ def _getRawProfile(self, star, skip=None):
         else:
             return galsim.Convolve(profiles), method
 
-    def _finish_write(self, fits, extname, logger):
+    def _finish_write(self, writer, logger):
         """Finish the writing process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
         logger = galsim.config.LoggerWrapper(logger)
@@ -285,32 +283,28 @@ def _finish_write(self, fits, extname, logger):
             'dof' : self.dof,
             'nremoved' : self.nremoved,
         }
-        write_kwargs(fits, extname + '_chisq', chisq_dict)
-        logger.debug("Wrote the chisq info to extension %s",extname + '_chisq')
+        writer.write_struct('chisq', chisq_dict)
+        logger.debug("Wrote the chisq info to extension %s", writer.get_full_name('chisq'))
         for k, comp in enumerate(self.components):
-            comp._write(fits, extname + '_' + str(k), logger=logger)
+            comp._write(writer, str(k), logger=logger)
         if self.outliers:
-            self.outliers.write(fits, extname + '_outliers')
-            logger.debug("Wrote the PSF outliers to extension %s",extname + '_outliers')
+            self.outliers.write(writer, 'outliers')
+            logger.debug("Wrote the PSF outliers to extension %s", writer.get_full_name('outliers'))
 
-    def _finish_read(self, fits, extname, logger):
+    def _finish_read(self, reader, logger):
         """Finish the reading process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param reader:      A reader object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        chisq_dict = read_kwargs(fits, extname + '_chisq')
+        chisq_dict = reader.read_struct('chisq')
         for key in chisq_dict:
             setattr(self, key, chisq_dict[key])
 
         ncomponents = self.components
         self.components = []
         for k in range(ncomponents):
-            self.components.append(PSF._read(fits, extname + '_' + str(k), logger=logger))
-        if extname + '_outliers' in fits:
-            self.outliers = Outliers.read(fits, extname + '_outliers')
-        else:
-            self.outliers = None
+            self.components.append(PSF._read(reader, str(k), logger=logger))
+        self.outliers = Outliers.read(reader, 'outliers')
         # Set up all the num's properly now that everything is constructed.
         self.set_num(None)
diff --git a/piff/des/decam_wavefront.py b/piff/des/decam_wavefront.py
index 843ce9db..ef34e07b 100644
--- a/piff/des/decam_wavefront.py
+++ b/piff/des/decam_wavefront.py
@@ -154,13 +154,12 @@ def _predict(self, locations, logger=None):
 
         return targets
 
-    def _finish_write(self, fits, extname):
-        """Write the solution to a FITS binary table.
+    def _finish_write(self, writer):
+        """Write the solution.
 
         Save the knn params and the locations and targets arrays
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension with the interp information.
+        :param writer:      A writer object that encapsulates the serialization format.
         """
 
         dtypes = [('LOCATIONS', self.locations.dtype, self.locations.shape),
@@ -175,15 +174,15 @@ def _finish_write(self, fits, extname):
         data['MISALIGNMENT'] = self.misalignment
 
         # write to fits
-        fits.write_table(data, extname=extname + '_solution')
+        writer.write_table('solution', data)
 
-    def _finish_read(self, fits, extname):
-        """Read the solution from a FITS binary table.
+    def _finish_read(self, reader):
+        """Read the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The name of the extension with the interp information.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
-        data = fits[extname + '_solution'].read()
+        data = reader.read_table('solution')
+        assert data is not None
 
         # self.locations and self.targets assigned in _fit
         self._fit(data['LOCATIONS'][0], data['TARGETS'][0])
diff --git a/piff/gp_interp.py b/piff/gp_interp.py
index 8b33c86c..b5ce2f87 100644
--- a/piff/gp_interp.py
+++ b/piff/gp_interp.py
@@ -292,11 +292,10 @@ def interpolateList(self, stars, logger=None):
             fitted_stars.append(Star(star.data, fit))
         return fitted_stars
 
-    def _finish_write(self, fits, extname):
+    def _finish_write(self, writer):
         """Finish the writing process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         # Note, we're only storing the training data and hyperparameters here, which means the
         # Cholesky decomposition will have to be re-computed when this object is read back from
@@ -320,15 +319,15 @@ def _finish_write(self, fits, extname):
         data['ROWS'] = self.rows
         data['OPTIMIZER'] = self.optimizer
 
-        fits.write_table(data, extname=extname+'_kernel')
+        writer.write_table('kernel', data)
 
-    def _finish_read(self, fits, extname):
+    def _finish_read(self, reader):
         """Finish the reading process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
-        data = fits[extname+'_kernel'].read()
+        data = reader.read_table('kernel')
+        assert data is not None
         # Run fit to set up GP, but don't actually do any hyperparameter optimization. Just
         # set the GP up using the current hyperparameters.
         # Need to give back average fits files if needed.
diff --git a/piff/interp.py b/piff/interp.py
index 22cde835..add8cbff 100644
--- a/piff/interp.py
+++ b/piff/interp.py
@@ -17,7 +17,6 @@
 """
 
 import numpy as np
-from .util import write_kwargs, read_kwargs
 
 class Interp(object):
     """The base class for interpolating a set of data vectors across the field of view.
@@ -174,77 +173,68 @@ def interpolateList(self, stars, logger=None):
         """
         return [ self.interpolate(star) for star in stars ]
 
-    def write(self, fits, extname):
-        """Write an Interp to a FITS file.
+    def write(self, writer, name):
+        """Write an Interp via a writer object.
 
         Note: this only writes the initialization kwargs to the fits extension, not the parameters.
 
-        The base class implemenation works if the class has a self.kwargs attribute and these
+        The base class implementation works if the class has a self.kwargs attribute and these
         are all simple values (str, float, or int).
 
         However, the derived class will need to implement _finish_write to write the solution
         parameters to a binary table.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to write the interpolator information.
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with this interpolator in the serialized output.
         """
         # First write the basic kwargs that works for all Interp classes
         interp_type = self.__class__._type_name
-        write_kwargs(fits, extname, dict(self.kwargs, type=interp_type))
+        writer.write_struct(name, dict(self.kwargs, type=interp_type))
 
         # Now do the class-specific steps.  Typically, this will write out the solution parameters.
-        self._finish_write(fits, extname)
+        with writer.nested(name) as w:
+            self._finish_write(w)
 
-    def _finish_write(self, fits, extname):
+    def _finish_write(self, writer):
         """Finish the writing process with any class-specific steps.
 
         The base class implementation doesn't do anything, but this will probably always be
         overridden by the derived class.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         raise NotImplementedError("Derived classes must define the _finish_write method.")
 
     @classmethod
-    def read(cls, fits, extname):
-        """Read an Interp from a FITS file.
+    def read(cls, reader, name):
+        """Read an Interp via a reader object.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension with the interpolator information.
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:        Name associated with this interpolator in the serialized output.
 
-        :returns: an interpolator built with a information in the FITS file.
+        :returns: an interpolator built from serialized information.
         """
-        assert extname in fits
-        assert 'type' in fits[extname].get_colnames()
-        assert 'type' in fits[extname].read().dtype.names
-        # interp_type = fits[extname].read_column('type')
-        interp_type = fits[extname].read()['type']
-        assert len(interp_type) == 1
-        try:
-            interp_type = str(interp_type[0].decode())
-        except AttributeError:
-            # fitsio 1.0 returns strings
-            interp_type = interp_type[0]
+        kwargs = reader.read_struct(name)
+        assert kwargs is not None
+        assert 'type' in kwargs
+        interp_type = kwargs.pop('type')
 
         # Check that interp_type is a valid Interp type.
         if interp_type not in Interp.valid_interp_types:
             raise ValueError("interp type %s is not a valid Piff Interpolation"%interp_type)
         interp_cls = Interp.valid_interp_types[interp_type]
 
-        kwargs = read_kwargs(fits, extname)
-        kwargs.pop('type',None)
         interp = interp_cls(**kwargs)
-        interp._finish_read(fits, extname)
+        with reader.nested(name) as r:
+            interp._finish_read(r)
         return interp
 
-    def _finish_read(self, fits, extname):
+    def _finish_read(self, reader):
         """Finish the reading process with any class-specific steps.
 
         The base class implementation doesn't do anything, but this will probably always be
         overridden by the derived class.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
         raise NotImplementedError("Derived classes must define the _finish_read method.")
diff --git a/piff/knn_interp.py b/piff/knn_interp.py
index fbbf5b5e..7e1966c2 100644
--- a/piff/knn_interp.py
+++ b/piff/knn_interp.py
@@ -156,13 +156,12 @@ def interpolateList(self, stars, logger=None):
             stars_fitted.append(Star(star.data, fit))
         return stars_fitted
 
-    def _finish_write(self, fits, extname):
-        """Write the solution to a FITS binary table.
+    def _finish_write(self, writer):
+        """Write the solution.
 
         Save the knn params and the locations and targets arrays
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension with the interp information.
+        :param writer:      A writer object that encapsulates the serialization format.
         """
 
         dtypes = [('LOCATIONS', self.locations.dtype, self.locations.shape),
@@ -173,16 +172,14 @@ def _finish_write(self, fits, extname):
         data['LOCATIONS'] = self.locations
         data['TARGETS'] = self.targets
 
-        # write to fits
-        fits.write_table(data, extname=extname + '_solution')
+        writer.write_table('solution', data)
 
-    def _finish_read(self, fits, extname):
-        """Read the solution from a FITS binary table.
+    def _finish_read(self, reader):
+        """Read the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension with the interp information.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
-        data = fits[extname + '_solution'].read()
+        data = reader.read_table('solution')
 
         # self.locations and self.targets assigned in _fit
         self._fit(data['LOCATIONS'][0], data['TARGETS'][0])
diff --git a/piff/mean_interp.py b/piff/mean_interp.py
index 1ba1ad0a..156468fb 100644
--- a/piff/mean_interp.py
+++ b/piff/mean_interp.py
@@ -60,22 +60,21 @@ def interpolate(self, star, logger=None):
             fit = star.fit.newParams(self.mean, num=self._num)
         return Star(star.data, fit)
 
-    def _finish_write(self, fits, extname):
-        """Write the solution to a FITS binary table.
+    def _finish_write(self, writer):
+        """Write the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         cols = [ self.mean ]
         dtypes = [ ('mean', float) ]
         data = np.array(list(zip(*cols)), dtype=dtypes)
-        fits.write_table(data, extname=extname + '_solution')
+        writer.write_table('solution', data)
 
-    def _finish_read(self, fits, extname):
-        """Read the solution from a FITS binary table.
+    def _finish_read(self, reader):
+        """Read the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension
+        :param reader:      A reader object that encapsulates the serialization format.
         """
-        data = fits[extname + '_solution'].read()
+        data = reader.read_table('solution')
+        assert data is not None
         self.mean = data['mean']
diff --git a/piff/model.py b/piff/model.py
index 393a4057..eef4a954 100644
--- a/piff/model.py
+++ b/piff/model.py
@@ -16,10 +16,6 @@
 .. module:: model
 """
 
-import numpy as np
-import galsim
-
-from .util import write_kwargs, read_kwargs
 from .star import Star
 
 
@@ -154,71 +150,63 @@ def draw(self, star, copy_image=True):
         prof.drawImage(image, method=self._method, center=star.image_pos)
         return Star(star.data.withNew(image=image), star.fit)
 
-    def write(self, fits, extname):
-        """Write a Model to a FITS file.
+    def write(self, writer, name):
+        """Write a Model via a Writer object.
 
         Note: this only writes the initialization kwargs to the fits extension, not the parameters.
 
         The base class implemenation works if the class has a self.kwargs attribute and these
         are all simple values (str, float, or int)
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to write the model information.
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with this model in the serialized output.
         """
         # First write the basic kwargs that works for all Model classes
         model_type = self._type_name
-        write_kwargs(fits, extname, dict(self.kwargs, type=model_type))
-
+        writer.write_struct(name, dict(self.kwargs, type=model_type))
         # Now do any class-specific steps.
-        self._finish_write(fits, extname)
+        with writer.nested(name) as w:
+            self._finish_write(w)
 
-    def _finish_write(self, fits, extname):
+    def _finish_write(self, writer):
         """Finish the writing process with any class-specific steps.
 
         The base class implementation doesn't do anything, which is often appropriate, but
         this hook exists in case any Model classes need to write extra information to the
         fits file.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         pass
 
     @classmethod
-    def read(cls, fits, extname):
+    def read(cls, reader, name):
         """Read a Model from a FITS file.
 
         Note: the returned Model will not have its parameters set.  This just initializes a fresh
         model that can be used to interpret interpolated vectors.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension with the model information.
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:        Name associated with this model in the serialized output.
 
-        :returns: a model built with a information in the FITS file.
+        :returns: a model built with a information in the FITS file
         """
-        assert extname in fits
-        assert 'type' in fits[extname].get_colnames()
-        model_type = fits[extname].read()['type']
-        assert len(model_type) == 1
-        try:
-            model_type = str(model_type[0].decode())
-        except AttributeError:
-            # fitsio 1.0 returns strings
-            model_type = model_type[0]
+        kwargs = reader.read_struct(name)
+        assert kwargs is not None
+        assert 'type' in kwargs
+        model_type = kwargs.pop('type')
 
         # Check that model_type is a valid Model type.
         if model_type not in Model.valid_model_types:
             raise ValueError("model type %s is not a valid Piff Model"%model_type)
         model_cls = Model.valid_model_types[model_type]
 
-        kwargs = read_kwargs(fits, extname)
-        kwargs.pop('type',None)
         if 'force_model_center' in kwargs: # pragma: no cover
             # old version of this parameter name.
             kwargs['centered'] = kwargs.pop('force_model_center')
         model_cls._fix_kwargs(kwargs)
         model = model_cls(**kwargs)
-        model._finish_read(fits, extname)
+        model._finish_read(reader)
         return model
 
     @classmethod
@@ -235,14 +223,13 @@ def _fix_kwargs(cls, kwargs):
         """
         pass
 
-    def _finish_read(self, fits, extname):
+    def _finish_read(self, reader):
         """Finish the reading process with any class-specific steps.
 
         The base class implementation doesn't do anything, which is often appropriate, but
         this hook exists in case any Model classes need to read extra information from the
         fits file.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
         pass
diff --git a/piff/outliers.py b/piff/outliers.py
index fbb6bf93..9f0b10fe 100644
--- a/piff/outliers.py
+++ b/piff/outliers.py
@@ -22,7 +22,6 @@
 import galsim
 from scipy.stats import chi2
 
-from .util import write_kwargs, read_kwargs
 
 class Outliers(object):
     """The base class for handling outliers.
@@ -91,49 +90,46 @@ def parseKwargs(cls, config_outliers, logger=None):
         kwargs['logger'] = logger
         return kwargs
 
-    def write(self, fits, extname):
-        """Write an Outliers to a FITS file.
+    def write(self, writer, name):
+        """Write an Outers via a writer object.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to write the outliers information.
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with the Ootliers in the serialized output.
         """
         # First write the basic kwargs that works for all Outliers classes
         outliers_type = self._type_name
-        write_kwargs(fits, extname, dict(self.kwargs, type=outliers_type))
+        writer.write_struct(name, dict(self.kwargs, type=outliers_type))
 
         # Now do any class-specific steps.
-        self._finish_write(fits, extname)
+        with writer.nested(name) as w:
+            self._finish_write(w)
 
-    def _finish_write(self, fits, extname):
+    def _finish_write(self, writer):
         """Finish the writing process with any class-specific steps.
 
         The base class implementation doesn't do anything, which is often appropriate, but
         this hook exists in case any Outliers classes need to write extra information to the
         fits file.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with the outliers in the serialized output.
         """
         pass
 
     @classmethod
-    def read(cls, fits, extname):
+    def read(cls, reader, name):
         """Read a Outliers from a FITS file.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension with the outliers information.
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:        Name associated with the outliers in the serialized output.
 
-        :returns: an Outliers handler
+        :returns: an Outliers handler, or None if there isn't one.
         """
-        assert extname in fits
-        assert 'type' in fits[extname].get_colnames()
-        outliers_type = fits[extname].read()['type']
-        assert len(outliers_type) == 1
-        try:
-            outliers_type = str(outliers_type[0].decode())
-        except AttributeError:
-            # fitsio 1.0 returns strings
-            outliers_type = outliers_type[0]
+        kwargs = reader.read_struct(name)
+        if kwargs is None:
+            return None
+        assert 'type' in kwargs
+        outliers_type = kwargs.pop('type')
 
         # Old output files had the full name.  Fix it if necessary.
         if (outliers_type.endswith('Outliers')
@@ -145,21 +141,19 @@ def read(cls, fits, extname):
             raise ValueError("outliers type %s is not a valid Piff Outliers"%outliers_type)
         outliers_cls = Outliers.valid_outliers_types[outliers_type]
 
-        kwargs = read_kwargs(fits, extname)
-        kwargs.pop('type',None)
         outliers = outliers_cls(**kwargs)
-        outliers._finish_read(fits, extname)
+        with reader.nested(name) as r:
+            outliers._finish_read(r)
         return outliers
 
-    def _finish_read(self, fits, extname):
+    def _finish_read(self, reader):
         """Finish the reading process with any class-specific steps.
 
         The base class implementation doesn't do anything, which is often appropriate, but
         this hook exists in case any Outliers classes need to read extra information from the
         fits file.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension.
+        :param reader:      A reader object that encapsulates the serialization format.
         """
         pass
 
diff --git a/piff/polynomial_interp.py b/piff/polynomial_interp.py
index 9e3c72ae..11b83880 100644
--- a/piff/polynomial_interp.py
+++ b/piff/polynomial_interp.py
@@ -320,11 +320,10 @@ def model(uv,*coeffs):
         # self._unpack_coefficients
         self.coeffs = coeffs
 
-    def _finish_write(self, fits, extname):
-        """Write the solution to a FITS binary table.
+    def _finish_write(self, writer):
+        """Write the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension
+        :param writer:      A writer object that encapsulates the serialization format.
         """
         if self.coeffs is None:
             raise RuntimeError("Coeffs not set yet.  Cannot write this Polynomial.")
@@ -368,20 +367,19 @@ def _finish_write(self, fits, extname):
 
         # Finally, write all of this to a FITS table.
         data = np.array(list(zip(*cols)), dtype=dtypes)
-        fits.write_table(data, extname=extname + '_solution', header=header)
+        writer.write_table('solution', data, metadata=header)
 
 
-    def _finish_read(self, fits, extname):
-        """Read the solution from a fits file.
+    def _finish_read(self, reader):
+        """Read the solution.
 
-        :param fits:        An open fitsio.FITS object.
-        :param extname:     The base name of the extension
+        :param reader:      A reader object that encapsulates the serialization format.
         """
         # Read the solution extension.
-        data = fits[extname + '_solution'].read()
-        header = fits[extname + '_solution'].read_header()
-
-        self.nparam = header['NPARAM']
+        metadata = {}
+        data = reader.read_table('solution', metadata=metadata)
+        assert data is not None
+        self.nparam = metadata['NPARAM']
 
         # Run setup functions to get these values right.
         self._set_function(self.poly_type)
diff --git a/piff/psf.py b/piff/psf.py
index fbb4ccd6..6d1c62b2 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -17,12 +17,9 @@
 """
 
 import numpy as np
-import fitsio
 import galsim
-import sys
 
 from .star import Star, StarData
-from .util import write_kwargs, read_kwargs
 
 class PSF(object):
     """The base class for describing a PSF model across a field of view.
@@ -722,34 +719,34 @@ def write(self, file_name, logger=None):
         :param file_name:   The name of the file to write to.
         :param logger:      A logger object for logging debug info. [default: None]
         """
+        from .writers import FitsWriter
+
         logger = galsim.config.LoggerWrapper(logger)
         logger.warning("Writing PSF to file %s",file_name)
 
-        with fitsio.FITS(file_name,'rw',clobber=True) as f:
-            self._write(f, 'psf', logger)
+        with FitsWriter.open(file_name) as w:
+            self._write(w, 'psf', logger=logger)
 
-    def _write(self, fits, extname, logger=None):
+    def _write(self, writer, name, logger):
         """This is the function that actually does the work for the write function.
         Composite PSF classes that need to iterate can call this multiple times as needed.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension with the psf information.
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with this PSF in the serialized output.
         :param logger:      A logger object for logging debug info.
         """
         from . import __version__ as piff_version
-        if len(fits) == 1:
-            header = {'piff_version': piff_version}
-            fits.write(data=None, header=header)
         psf_type = self._type_name
-        write_kwargs(fits, extname, dict(self.kwargs, type=psf_type, piff_version=piff_version))
-        logger.info("Wrote the basic PSF information to extname %s", extname)
-        if hasattr(self, 'stars'):
-            Star.write(self.stars, fits, extname=extname + '_stars')
-            logger.info("Wrote the PSF stars to extname %s", extname + '_stars')
-        if hasattr(self, 'wcs'):
-            self.writeWCS(fits, extname=extname + '_wcs', logger=logger)
-            logger.info("Wrote the PSF WCS to extname %s", extname + '_wcs')
-        self._finish_write(fits, extname=extname, logger=logger)
+        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))
+        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'))
+            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'))
+            self._finish_write(w, logger=logger)
 
     @classmethod
     def read(cls, file_name, logger=None):
@@ -760,26 +757,27 @@ def read(cls, file_name, logger=None):
 
         :returns: a PSF instance
         """
+        from .readers import FitsReader
+
         logger = galsim.config.LoggerWrapper(logger)
         logger.warning("Reading PSF from file %s",file_name)
 
-        with fitsio.FITS(file_name,'r') as f:
+        with FitsReader.open(file_name) as r:
             logger.debug('opened FITS file')
-            return cls._read(f, 'psf', logger)
+            return cls._read(r, 'psf', logger)
 
     @classmethod
-    def _read(cls, fits, extname, logger):
+    def _read(cls, reader, name, logger):
         """This is the function that actually does the work for the read function.
         Composite PSF classes that need to iterate can call this multiple times as needed.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension with the psf information.
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:        Name associated with this PSF in the serialized output.
         :param logger:      A logger object for logging debug info.
         """
         # Read the type and kwargs from the base extension
-        assert extname in fits
-        assert 'type' in fits[extname].get_colnames()
-        kwargs = read_kwargs(fits, extname)
+        kwargs = reader.read_struct(name)
+        assert kwargs is not None
         psf_type = kwargs.pop('type')
 
         # Old output files had the full class name.  Fix it if necessary.
@@ -798,144 +796,26 @@ def _read(cls, fits, extname, logger):
         # Make the PSF instance
         psf = psf_cls(**kwargs)
 
-        # Read the stars, wcs, pointing values
-        if extname + '_stars' in fits:
-            stars = Star.read(fits, extname + '_stars')
-            logger.debug("stars = %s",stars)
-            psf.stars = stars
-        if extname + '_wcs' in fits:
-            wcs, pointing = cls.readWCS(fits, extname + '_wcs', logger=logger)
-            logger.debug("wcs = %s, pointing = %s",wcs,pointing)
-            psf.wcs = wcs
-            psf.pointing = pointing
-
-        # Just in case the class needs to do something else at the end.
-        psf._finish_read(fits, extname, logger)
-
-        # Save the piff version as an attibute.
+        with reader.nested(name) as r:
+            # Read the stars, wcs, pointing values
+            stars = Star.read(r, 'stars')
+            if stars is not None:
+                logger.debug("stars = %s", stars)
+                psf.stars = stars
+            wcs, pointing = r.read_wcs_map('wcs', logger=logger)
+            if wcs is not None:
+                logger.debug("wcs = %s, pointing = %s",wcs,pointing)
+                psf.wcs = wcs
+                psf.pointing = pointing
+
+            # Just in case the class needs to do something else at the end.
+            psf._finish_read(r, logger)
+
+        # Save the piff version as an attribute.
         psf.piff_version = piff_version
 
         return psf
 
-    def writeWCS(self, fits, extname, logger):
-        """Write the WCS information to a FITS file.
-
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to write to
-        :param logger:      A logger object for logging debug info.
-        """
-        import base64
-        try:
-            import cPickle as pickle
-        except ImportError:
-            import pickle
-        logger = galsim.config.LoggerWrapper(logger)
-
-        # Start with the chipnums
-        chipnums = list(self.wcs.keys())
-        cols = [ chipnums ]
-        dtypes = [ ('chipnums', int) ]
-
-        # GalSim WCS objects can be serialized via pickle
-        wcs_str = [ base64.b64encode(pickle.dumps(w)) for w in self.wcs.values() ]
-        max_len = np.max([ len(s) for s in wcs_str ])
-        # Some GalSim WCS serializations are rather long.  In particular, the Pixmappy one
-        # is longer than the maximum length allowed for a column in a fits table (28799).
-        # So split it into chunks of size 2**14 (mildly less than this maximum).
-        chunk_size = 2**14
-        nchunks = max_len // chunk_size + 1
-        cols.append( [nchunks]*len(chipnums) )
-        dtypes.append( ('nchunks', int) )
-
-        # Update to size of chunk we actually need.
-        chunk_size = (max_len + nchunks - 1) // nchunks
-
-        chunks = [ [ s[i:i+chunk_size] for i in range(0, max_len, chunk_size) ] for s in wcs_str ]
-        cols.extend(zip(*chunks))
-        dtypes.extend( ('wcs_str_%04d'%i, bytes, chunk_size) for i in range(nchunks) )
-
-        if self.pointing is not None:
-            # Currently, there is only one pointing for all the chips, but write it out
-            # for each row anyway.
-            dtypes.extend( (('ra', float), ('dec', float)) )
-            ra = [self.pointing.ra / galsim.hours] * len(chipnums)
-            dec = [self.pointing.dec / galsim.degrees] * len(chipnums)
-            cols.extend( (ra, dec) )
-
-        data = np.array(list(zip(*cols)), dtype=dtypes)
-        fits.write_table(data, extname=extname)
-
-    @classmethod
-    def readWCS(cls, fits, extname, logger):
-        """Read the WCS information from a FITS file.
-
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to read from
-        :param logger:      A logger object for logging debug info.
-
-        :returns: wcs, pointing where wcs is a dict of galsim.BaseWCS instances and
-                                      pointing is a galsim.CelestialCoord instance
-        """
-        import base64
-        try:
-            import cPickle as pickle
-        except ImportError:
-            import pickle
-
-        assert extname in fits
-        assert 'chipnums' in fits[extname].get_colnames()
-        assert 'nchunks' in fits[extname].get_colnames()
-
-        data = fits[extname].read()
-
-        chipnums = data['chipnums']
-        nchunks = data['nchunks']
-        nchunks = nchunks[0]  # These are all equal, so just take first one.
-
-        wcs_keys = [ 'wcs_str_%04d'%i for i in range(nchunks) ]
-        wcs_str = [ data[key] for key in wcs_keys ] # Get all wcs_str columns
-        try:
-            wcs_str = [ b''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
-        except TypeError:  # pragma: no cover
-            # fitsio 1.0 returns strings
-            wcs_str = [ ''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
-
-        wcs_str = [ base64.b64decode(s) for s in wcs_str ] # Convert back from b64 encoding
-        # Convert back into wcs objects
-        try:
-            wcs_list = [ pickle.loads(s, encoding='bytes') for s in wcs_str ]
-        except Exception:
-            # If the file was written by py2, the bytes encoding might raise here,
-            # or it might not until we try to use it.
-            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
-
-        wcs = dict(zip(chipnums, wcs_list))
-
-        try:
-            # 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.')
-            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
-            wcs = dict(zip(chipnums, wcs_list))
-            repr(wcs)
-
-        # Work-around for a change in the GalSim API with 2.0
-        # If the piff file was written with pre-2.0 GalSim, this fixes it.
-        for key in wcs:
-            w = wcs[key]
-            if hasattr(w, '_origin') and  isinstance(w._origin, galsim._galsim.PositionD):
-                w._origin = galsim.PositionD(w._origin)
-
-        if 'ra' in fits[extname].get_colnames():
-            ra = data['ra']
-            dec = data['dec']
-            pointing = galsim.CelestialCoord(ra[0] * galsim.hours, dec[0] * galsim.degrees)
-        else:
-            pointing = None
-
-        return wcs, pointing
 
 # Make a global function, piff.read, as an alias for piff.PSF.read, since that's the main thing
 # users will want to do as their starting point for using a piff file.
diff --git a/piff/readers.py b/piff/readers.py
new file mode 100644
index 00000000..3a0c4cd0
--- /dev/null
+++ b/piff/readers.py
@@ -0,0 +1,177 @@
+# Copyright (c) 2024 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.
+
+"""
+.. module:: readers
+"""
+
+from contextlib import contextmanager
+
+import fitsio
+import galsim
+
+
+class FitsReader:
+    """A reader object that reads from to multiple FITS HDUs.
+
+    This reader is intended to read files written by `FitsWriter`.
+
+    :param fits:        An already-open `fitsio.FITS` object.
+    :param base_name:   Base name to prepend to all object names, or `None`.
+    """
+    def __init__(self, fits, base_name):
+        self._fits = fits
+        self._base_name = base_name
+
+    @classmethod
+    @contextmanager
+    def open(cls, file_name):
+        """Return a context manager that opens the given FITS file and yields
+        a `FitsReader` object.
+
+        :param file_name:   Name of the file to read.
+
+        :returns:  A context manager that yields a `FitsReader` instance.
+        """
+        with fitsio.FITS(file_name, "r") as f:
+            yield cls(f, base_name=None)
+
+    def read_struct(self, name):
+        """Load a simple flat dictionary.
+
+        :param name:   Name used to save this struct in `FitsWriter.write_struct`.
+
+        :returns:  A `dict` with `str` keys and `int`, `float`, `str`, `bool`,
+                   or `None` values.  If nothing was stored with the given
+                   name, `None` is returned.
+        """
+        extname = self.get_full_name(name)
+        if extname not in self._fits:
+            return None
+        cols = self._fits[extname].get_colnames()
+        data = self._fits[extname].read()
+        assert len(data) == 1
+        return dict([ (col, data[col][0]) for col in cols ])
+
+    def read_table(self, name, metadata=None):
+        """Load a table as a numpy array.
+
+        :param name:      Name used to save this table in `FitsWriter.write_table`.
+        :param metadata:  If not `None`, a `dict` to be filled with any
+                          metadata associated with the table on write.  Key
+                          case may not be preserved!
+
+        :returns:  A numpy array with a structured dtype, or `None` if no
+                   object with this name was saved.
+        """
+        extname = self.get_full_name(name)
+        if extname not in self._fits:
+            return None
+        if metadata is not None:
+            metadata.update(self._fits[extname].read_header())
+        return self._fits[extname].read()
+
+    def read_wcs_map(self, name, logger):
+        """Load a map of WCS objects and an optional pointing coord.
+
+        :param name:      Name used to save this map in `FitsWriter.write_wcs_map`.
+        :param logger:    A logger object for logging debug info.
+
+        :returns:  A 2-element tuple, where the first entry is a `dict` with
+                   `int` (chipnum) keys and `galsim.BaseWCS` values, and the
+                   second entry is a `galsim.CelestialCoord` object.  Both
+                   entries may be `None`, and if the WCS `dict` is `None` the
+                   pointing coord is always `None`.
+        """
+        import base64
+        import pickle
+
+        extname = self.get_full_name(name)
+        if extname not in self._fits:
+            return None, None
+
+        assert 'chipnums' in self._fits[extname].get_colnames()
+        assert 'nchunks' in self._fits[extname].get_colnames()
+
+        data = self._fits[extname].read()
+
+        chipnums = data['chipnums']
+        nchunks = data['nchunks']
+        nchunks = nchunks[0]  # These are all equal, so just take first one.
+
+        wcs_keys = [ 'wcs_str_%04d'%i for i in range(nchunks) ]
+        wcs_str = [ data[key] for key in wcs_keys ] # Get all wcs_str columns
+        try:
+            wcs_str = [ b''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
+        except TypeError:  # pragma: no cover
+            # fitsio 1.0 returns strings
+            wcs_str = [ ''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
+
+        wcs_str = [ base64.b64decode(s) for s in wcs_str ] # Convert back from b64 encoding
+        # Convert back into wcs objects
+        try:
+            wcs_list = [ pickle.loads(s, encoding='bytes') for s in wcs_str ]
+        except Exception:
+            # If the file was written by py2, the bytes encoding might raise here,
+            # or it might not until we try to use it.
+            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
+
+        wcs = dict(zip(chipnums, wcs_list))
+
+        try:
+            # 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.')
+            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
+            wcs = dict(zip(chipnums, wcs_list))
+            repr(wcs)
+
+        # Work-around for a change in the GalSim API with 2.0
+        # If the piff file was written with pre-2.0 GalSim, this fixes it.
+        for key in wcs:
+            w = wcs[key]
+            if hasattr(w, '_origin') and isinstance(w._origin, galsim._galsim.PositionD):
+                w._origin = galsim.PositionD(w._origin)
+
+        if 'ra' in self._fits[extname].get_colnames():
+            ra = data['ra']
+            dec = data['dec']
+            pointing = galsim.CelestialCoord(ra[0] * galsim.hours, dec[0] * galsim.degrees)
+        else:
+            pointing = None
+
+        return wcs, pointing
+
+    @contextmanager
+    def nested(self, name):
+        """Return a context manager that yields a new `FitsReader` that reads
+        content written by a similarly nested `FitsWriter`.
+
+        :param name:     Base name for all objects read with the returned object.
+
+        :returns:     A context manager that yields a nested reader object.
+        """
+        yield FitsReader(self._fits, base_name=self.get_full_name(name))
+
+    def get_full_name(self, name: str) -> str:
+        """Return the full name of a data structure saved with the given name,
+        combining it with any base names added if this `FitsReader` was created
+        by the `nested` method.
+
+        The FITS implementation concatenates with ``_`` as a delimiter, and
+        uses the full name as the EXTNAME header value.
+        """
+        return name if self._base_name is None else f"{self._base_name}_{name}"
diff --git a/piff/simplepsf.py b/piff/simplepsf.py
index b48a1de7..0992aaea 100644
--- a/piff/simplepsf.py
+++ b/piff/simplepsf.py
@@ -23,7 +23,6 @@
 from .interp import Interp
 from .outliers import Outliers
 from .psf import PSF
-from .util import write_kwargs, read_kwargs
 
 class SimplePSF(PSF):
     """A PSF class that uses a single model and interpolator.
@@ -242,11 +241,10 @@ def _drawStar(self, star):
     def _getRawProfile(self, star):
         return self.model.getProfile(star.fit.get_params(self._num)), self.model._method
 
-    def _finish_write(self, fits, extname, logger):
+    def _finish_write(self, writer, logger):
         """Finish the writing process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
         logger = galsim.config.LoggerWrapper(logger)
@@ -257,29 +255,26 @@ def _finish_write(self, fits, extname, logger):
             'nremoved' : self.nremoved,
             'niter' : self.niter,
         }
-        write_kwargs(fits, extname + '_chisq', chisq_dict)
-        logger.debug("Wrote the chisq info to extension %s",extname + '_chisq')
-        self.model.write(fits, extname + '_model')
-        logger.debug("Wrote the PSF model to extension %s",extname + '_model')
-        self.interp.write(fits, extname + '_interp')
-        logger.debug("Wrote the PSF interp to extension %s",extname + '_interp')
+        writer.write_struct('chisq', chisq_dict)
+        logger.debug("Wrote the chisq info to %s", writer.get_full_name('chisq'))
+        self.model.write(writer, 'model')
+        logger.debug("Wrote the PSF model to %s", writer.get_full_name('model'))
+        self.interp.write(writer, 'interp')
+        logger.debug("Wrote the PSF interp to %s", writer.get_full_name('interp'))
         if self.outliers:
-            self.outliers.write(fits, extname + '_outliers')
-            logger.debug("Wrote the PSF outliers to extension %s",extname + '_outliers')
+            self.outliers.write(writer, 'outliers')
+            logger.debug("Wrote the PSF outliers to %s", writer.get_full_name('outliers'))
 
-    def _finish_read(self, fits, extname, logger):
+    def _finish_read(self, reader, logger):
         """Finish the reading process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:        Name associated with this PSF in the serialized output.
         :param logger:      A logger object for logging debug info.
         """
-        chisq_dict = read_kwargs(fits, extname + '_chisq')
+        chisq_dict = reader.read_struct('chisq')
         for key in chisq_dict:
             setattr(self, key, chisq_dict[key])
-        self.model = Model.read(fits, extname + '_model')
-        self.interp = Interp.read(fits, extname + '_interp')
-        if extname + '_outliers' in fits:
-            self.outliers = Outliers.read(fits, extname + '_outliers')
-        else:
-            self.outliers = None
+        self.model = Model.read(reader, 'model')
+        self.interp = Interp.read(reader, 'interp')
+        self.outliers = Outliers.read(reader, 'outliers')
diff --git a/piff/singlechip.py b/piff/singlechip.py
index f9fdff79..6bc088fc 100644
--- a/piff/singlechip.py
+++ b/piff/singlechip.py
@@ -21,7 +21,7 @@
 import galsim
 
 from .psf import PSF
-from .util import write_kwargs, read_kwargs, make_dtype, adjust_value, run_multi
+from .util import make_dtype, adjust_value, run_multi
 
 # Used by SingleChipPSF.fit
 def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_funcs, draw_method, logger):
@@ -178,14 +178,13 @@ def _getRawProfile(self, star):
         chipnum = star['chipnum']
         return self.psf_by_chip[chipnum]._getRawProfile(star)
 
-    def _finish_write(self, fits, extname, logger):
+    def _finish_write(self, writer, logger):
         """Finish the writing process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        # Write the colnums to an extension.
+        # Write the colnums to a table.
         chipnums = list(self.psf_by_chip.keys())
         chipnums = [c for c in chipnums if self.psf_by_chip[c] is not None]
         dt = make_dtype('chipnums', chipnums[0])
@@ -193,21 +192,22 @@ def _finish_write(self, fits, extname, logger):
         cols = [ chipnums ]
         dtypes = [ dt ]
         data = np.array(list(zip(*cols)), dtype=dtypes)
-        fits.write_table(data, extname=extname + '_chipnums')
+        writer.write_table('chipnums', data)
 
-        # Add _1, _2, etc. to the extname for the psf model of each chip.
+        # Append 1, 2, etc. to the name for the psf model of each chip.
         for chipnum in chipnums:
-            self.psf_by_chip[chipnum]._write(fits, extname + '_%s'%chipnum, logger)
+            self.psf_by_chip[chipnum]._write(writer, str(chipnum), logger)
 
-    def _finish_read(self, fits, extname, logger):
+    def _finish_read(self, reader, logger):
         """Finish the reading process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param reader:      A reader object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        chipnums = fits[extname + '_chipnums'].read()['chipnums']
+        table = reader.read_table('chipnums')
+        assert table is not None
+        chipnums = table['chipnums']
         self.psf_by_chip = {}
         for chipnum in chipnums:
-            self.psf_by_chip[chipnum] = PSF._read(fits, extname + '_%s'%chipnum, logger)
+            self.psf_by_chip[chipnum] = PSF._read(reader, str(chipnum), logger)
         self.single_psf = self.psf_by_chip[chipnums[0]]
diff --git a/piff/star.py b/piff/star.py
index 0ceda3bd..8ae5b2b5 100644
--- a/piff/star.py
+++ b/piff/star.py
@@ -328,12 +328,12 @@ def makeTarget(cls, x=None, y=None, u=None, v=None, properties={}, wcs=None, sca
         return cls(data, fit)
 
     @classmethod
-    def write(self, stars, fits, extname):
-        """Write a list of stars to a FITS file.
+    def write(cls, stars, writer, name):
+        """Write a list of stars to a writer object.
 
         :param stars:       A list of stars to write
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to write to
+        :param writer:      A writer object that encapsulates the serialization format.
+        :param name:        A name to associate with these stars in the serialized output.
         """
         # TODO This doesn't write everything out.  Probably want image as an optional I/O.
 
@@ -406,7 +406,7 @@ def write(self, stars, fits, extname):
             cols.append( [s.data.pointing.dec / galsim.degrees for s in stars ] )
 
         data = np.array(list(zip(*cols)), dtype=dtypes)
-        fits.write_table(data, extname=extname, header=header)
+        writer.write_table(name, data, metadata=header)
 
     @classmethod
     def read_coords_params(cls, fits, extname):
@@ -432,19 +432,21 @@ def read_coords_params(cls, fits, extname):
         return coords, params
 
     @classmethod
-    def read(cls, fits, extname):
-        """Read stars from a FITS file.
-
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The name of the extension to read from
+    def read(cls, reader, name):
+        """Read a list of stars via an open reader object.
 
-        :returns: a list of Star instances
+        :param reader:      A reader object that encapsulates the serialization format.
+        :param name:         Name associated with the stars in the serialized output.
+        :returns: a list of Star instances, or None if there aren't any
         """
-        assert extname in fits
-        colnames = fits[extname].get_colnames()
-        header = fits[extname].read_header()
-        if 'PARAMS_LENS' in header:
-            params_lens = np.atleast_1d(eval(header['PARAMS_LENS']))
+        metadata = {}
+        data = reader.read_table(name, metadata)
+        if data is None:
+            return None
+
+        colnames = list(data.dtype.names)
+        if 'PARAMS_LENS' in metadata:
+            params_lens = np.atleast_1d(eval(metadata['PARAMS_LENS']))
         else:
             params_lens = None
 
@@ -457,7 +459,6 @@ def read(cls, fits, extname):
         # These two might not be there, but if they are, remove them.
         colnames = [key for key in colnames if key not in ['reserve', 'flag_psf']]
 
-        data = fits[extname].read()
         x_list = data['x']
         y_list = data['y']
         u_list = data['u']
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 75c9d865..59eb4a5f 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -20,7 +20,6 @@
 import galsim
 
 from .psf import PSF
-from .util import write_kwargs, read_kwargs
 from .star import Star, StarFit
 from .outliers import Outliers
 
@@ -288,11 +287,10 @@ def _getRawProfile(self, star):
         # Add them up.
         return galsim.Sum(profiles), method
 
-    def _finish_write(self, fits, extname, logger):
+    def _finish_write(self, writer, logger):
         """Finish the writing process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param writer:      A writer object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
         logger = galsim.config.LoggerWrapper(logger)
@@ -302,32 +300,28 @@ def _finish_write(self, fits, extname, logger):
             'dof' : self.dof,
             'nremoved' : self.nremoved,
         }
-        write_kwargs(fits, extname + '_chisq', chisq_dict)
-        logger.debug("Wrote the chisq info to extension %s",extname + '_chisq')
+        writer.write_struct('chisq', chisq_dict)
+        logger.debug("Wrote the chisq info to %s", writer.get_full_name('chisq'))
         for k, comp in enumerate(self.components):
-            comp._write(fits, extname + '_' + str(k), logger=logger)
+            comp._write(writer, str(k), logger=logger)
         if self.outliers:
-            self.outliers.write(fits, extname + '_outliers')
-            logger.debug("Wrote the PSF outliers to extension %s",extname + '_outliers')
+            self.outliers.write(writer, 'outliers')
+            logger.debug("Wrote the PSF outliers to %s", writer.get_full_name('outliers'))
 
-    def _finish_read(self, fits, extname, logger):
+    def _finish_read(self, reader, logger):
         """Finish the reading process with any class-specific steps.
 
-        :param fits:        An open fitsio.FITS object
-        :param extname:     The base name of the extension to write to.
+        :param reader:      A reader object that encapsulates the serialization format.
         :param logger:      A logger object for logging debug info.
         """
-        chisq_dict = read_kwargs(fits, extname + '_chisq')
+        chisq_dict = reader.read_struct('chisq')
         for key in chisq_dict:
             setattr(self, key, chisq_dict[key])
 
         ncomponents = self.components
         self.components = []
         for k in range(ncomponents):
-            self.components.append(PSF._read(fits, extname + '_' + str(k), logger=logger))
-        if extname + '_outliers' in fits:
-            self.outliers = Outliers.read(fits, extname + '_outliers')
-        else:
-            self.outliers = None
+            self.components.append(PSF._read(reader, str(k), logger=logger))
+        self.outliers = Outliers.read(reader, 'outliers')
         # Set up all the num's properly now that everything is constructed.
         self.set_num(None)
diff --git a/piff/util.py b/piff/util.py
index 7503e9f5..78646af5 100644
--- a/piff/util.py
+++ b/piff/util.py
@@ -101,44 +101,6 @@ def adjust_value(value, dtype):
             # For other numpy arrays, we can use astype instead.
             return np.array(value).astype(t)
 
-def write_kwargs(fits, extname, kwargs):
-    """A helper function for writing a single row table into a fits file with the values
-    and column names given by a kwargs dict.
-
-    :param fits:        An open fitsio.FITS instance
-    :param extname:     The extension to write to
-    :param kwargs:      A kwargs dict to be written as a FITS binary table.
-    """
-    from . import __version__ as piff_version
-    cols = []
-    dtypes = []
-    for key, value in kwargs.items():
-        # Don't add values that are None to the table.
-        if value is None:
-            continue
-        dt = make_dtype(key, value)
-        value = adjust_value(value,dt)
-        cols.append([value])
-        dtypes.append(dt)
-    data = np.array(list(zip(*cols)), dtype=dtypes)
-    header = {'piff_version': piff_version}
-    fits.write_table(data, extname=extname, header=header)
-
-def read_kwargs(fits, extname):
-    """A helper function for reading a single row table from a fits file returning the values
-    and column names as a kwargs dict.
-
-    :param fits:        An open fitsio.FITS instance
-    :param extname:     The extension to read.
-
-    :returns: A dict of the kwargs that were read from the file.
-    """
-    cols = fits[extname].get_colnames()
-    data = fits[extname].read()
-    assert len(data) == 1
-    kwargs = dict([ (col, data[col][0]) for col in cols ])
-    return kwargs
-
 def estimate_cov_from_jac(jac):
     """Estimate a covariance matrix from a jacobian as returned by scipy.optimize.least_squares
     .. math::
diff --git a/piff/writers.py b/piff/writers.py
new file mode 100644
index 00000000..1c6c2597
--- /dev/null
+++ b/piff/writers.py
@@ -0,0 +1,165 @@
+# Copyright (c) 2024 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.
+
+"""
+.. module:: writers
+"""
+
+from contextlib import contextmanager
+
+import fitsio
+import galsim
+import numpy as np
+
+from .util import make_dtype, adjust_value
+
+
+class FitsWriter:
+    """A writer object that writes to multiple FITS HDUs.
+
+    :param fits:        An already-open `fitsio.FITS` object.
+    :param base_name:   Base name to prepend to all object names, or `None`.
+    :param header:      Fields to be added to all FITS extension headers.
+    """
+
+    def __init__(self, fits, base_name, header):
+        self._fits = fits
+        self._base_name = base_name
+        self._header = header
+
+    @classmethod
+    @contextmanager
+    def open(cls, file_name):
+        """Return a context manager that opens the given FITS file and yields
+        a `FitsWriter` object.
+
+        :param file_name:   Name of the file to read.
+
+        :returns:  A context manager that yields a `FitsWriter` instance.
+
+        This also adds an empty primary HDU and sets up the returned writer
+        to add the Piff version to all HDUs (including the primary).
+        """
+        from . import __version__ as piff_version
+
+        header = {"piff_version": piff_version}
+        with fitsio.FITS(file_name, "rw", clobber=True) as f:
+            if len(f) == 1:
+                f.write(data=None, header=header)
+            yield cls(f, base_name=None, header=header.copy())
+
+    def write_struct(self, name, struct):
+        """Write a simple flat dictionary.
+
+        :param name:        Name used to save this struct.
+        :param struct:      A `dict` with `str` keys and `int`, `float`, `str`,
+                            `bool`, or `None` values.
+        """
+        cols = []
+        dtypes = []
+        for key, value in struct.items():
+            # Don't add values that are None to the table.
+            if value is None:
+                continue
+            dt = make_dtype(key, value)
+            value = adjust_value(value, dt)
+            cols.append([value])
+            dtypes.append(dt)
+        table = np.array(list(zip(*cols)), dtype=dtypes)
+        return self.write_table(name, table)
+
+    def write_table(self, name, array, metadata=None):
+        """Write a table via a numpy array.
+
+        :param name:      Name used to save this table.
+        :param metadata:  A `dict` of simple metadata to save with the table.
+                          Keys must be `str` (case may not be preserved) and
+                          values must be `int`, `float`, `str`, or `bool`.
+        """
+        if metadata:
+            header = self._header.copy()
+            header.update(metadata)
+        else:
+            header = self._header
+        self._fits.write_table(array, extname=self.get_full_name(name), header=header)
+
+    def write_wcs_map(self, name, wcs_map, pointing):
+        """Write a regular a map of WCS objects and an optoinal pointing coord.
+
+        :param name:      Name used to save this WCS map.
+        :param wcs_map:   A `dict` mapping `int` chipnum to `galsim.BaseWCS`.
+        :param pointing:  A `galsim.CelestialCoord`, or `None`.
+        :param logger:    A logger object for logging debug info.
+        """
+        import base64
+        import pickle
+        # Start with the chipnums
+        chipnums = list(wcs_map.keys())
+        cols = [chipnums]
+        dtypes = [("chipnums", int)]
+
+        # GalSim WCS objects can be serialized via pickle
+        wcs_str = [base64.b64encode(pickle.dumps(w)) for w in wcs_map.values()]
+        max_len = np.max([len(s) for s in wcs_str])
+        # Some GalSim WCS serializations are rather long.  In particular, the Pixmappy one
+        # is longer than the maximum length allowed for a column in a fits table (28799).
+        # So split it into chunks of size 2**14 (mildly less than this maximum).
+        chunk_size = 2**14
+        nchunks = max_len // chunk_size + 1
+        cols.append([nchunks] * len(chipnums))
+        dtypes.append(("nchunks", int))
+
+        # Update to size of chunk we actually need.
+        chunk_size = (max_len + nchunks - 1) // nchunks
+
+        chunks = [
+            [s[i : i + chunk_size] for i in range(0, max_len, chunk_size)]
+            for s in wcs_str
+        ]
+        cols.extend(zip(*chunks))
+        dtypes.extend(("wcs_str_%04d" % i, bytes, chunk_size) for i in range(nchunks))
+
+        if pointing is not None:
+            # Currently, there is only one pointing for all the chips, but write it out
+            # for each row anyway.
+            dtypes.extend((("ra", float), ("dec", float)))
+            ra = [pointing.ra / galsim.hours] * len(chipnums)
+            dec = [pointing.dec / galsim.degrees] * len(chipnums)
+            cols.extend((ra, dec))
+
+        data = np.array(list(zip(*cols)), dtype=dtypes)
+        self.write_table(name, data)
+
+    @contextmanager
+    def nested(self, name):
+        """Return a context manager that yields a new `FitsWriter` that nests
+        all names it is given within this one.
+
+        :param name:     Base name for all objects written with the returned object.
+
+        :returns:     A context manager that yields a nested writer object.
+        """
+        yield FitsWriter(
+            self._fits, base_name=self.get_full_name(name), header=self._header
+        )
+
+    def get_full_name(self, name):
+        """Return the full name of a data structure saved with the given name,
+        combining it with any base names added if this `FitsWriter` was created
+        by the `nested` method.
+
+        The FITS implementation concatenates with ``_`` as a delimiter, and
+        uses the full name as the EXTNAME header value.
+        """
+        return name if self._base_name is None else f"{self._base_name}_{name}"


Test patch:
diff --git a/tests/test_gp_interp.py b/tests/test_gp_interp.py
index 16f740e6..8b35ad9a 100644
--- a/tests/test_gp_interp.py
+++ b/tests/test_gp_interp.py
@@ -224,10 +224,10 @@ def check_gp(stars_training, stars_validation, kernel, optimizer,
 
     # Check I/O.
     file_name = os.path.join('output', 'test_gp.fits')
-    with fitsio.FITS(file_name,'rw',clobber=True) as fout:
-        interp.write(fout, extname='gp')
-    with fitsio.FITS(file_name,'r') as fin:
-        interp2 = piff.Interp.read(fin, extname='gp')
+    with piff.writers.FitsWriter.open(file_name) as w:
+        interp.write(w, 'gp')
+    with piff.readers.FitsReader.open(file_name) as r:
+        interp2 = piff.Interp.read(r, 'gp')
 
     stars_test = interp2.interpolateList(stars_validation)
     y_test = np.array([star.fit.params for star in stars_test])
diff --git a/tests/test_gsobject_model.py b/tests/test_gsobject_model.py
index 1c05d4b4..6fe50cce 100644
--- a/tests/test_gsobject_model.py
+++ b/tests/test_gsobject_model.py
@@ -163,10 +163,10 @@ def test_simple():
 
         # Also need to test ability to serialize
         outfile = os.path.join('output', 'gsobject_test.fits')
-        with fitsio.FITS(outfile, 'rw', clobber=True) as f:
-            model.write(f, 'psf_model')
-        with fitsio.FITS(outfile, 'r') as f:
-            roundtrip_model = piff.GSObjectModel.read(f, 'psf_model')
+        with piff.writers.FitsWriter.open(outfile) as w:
+            model.write(w, 'psf_model')
+        with piff.readers.FitsReader.open(outfile) as r:
+            roundtrip_model = piff.GSObjectModel.read(r, 'psf_model')
         assert model.__dict__ == roundtrip_model.__dict__
 
         # Check the deprecated name in config
@@ -903,10 +903,10 @@ def test_direct():
 
         # Also need to test ability to serialize
         outfile = os.path.join('output', 'gsobject_direct_test.fits')
-        with fitsio.FITS(outfile, 'rw', clobber=True) as f:
-            model.write(f, 'psf_model')
-        with fitsio.FITS(outfile, 'r') as f:
-            roundtrip_model = piff.GSObjectModel.read(f, 'psf_model')
+        with piff.writers.FitsWriter.open(outfile) as w:
+            model.write(w, 'psf_model')
+        with piff.readers.FitsReader.open(outfile) as r:
+            roundtrip_model = piff.GSObjectModel.read(r, 'psf_model')
         assert model.__dict__ == roundtrip_model.__dict__
 
     # repeat with fastfit=False
@@ -958,10 +958,10 @@ def test_direct():
 
         # Also need to test ability to serialize
         outfile = os.path.join('output', 'gsobject_direct_test.fits')
-        with fitsio.FITS(outfile, 'rw', clobber=True) as f:
-            model.write(f, 'psf_model')
-        with fitsio.FITS(outfile, 'r') as f:
-            roundtrip_model = piff.GSObjectModel.read(f, 'psf_model')
+        with piff.writers.FitsWriter.open(outfile) as w:
+            model.write(w, 'psf_model')
+        with piff.readers.FitsReader.open(outfile) as r:
+            roundtrip_model = piff.GSObjectModel.read(r, 'psf_model')
         assert model.__dict__ == roundtrip_model.__dict__
 
 @timer
diff --git a/tests/test_knn_interp.py b/tests/test_knn_interp.py
index ec30faf5..3605e160 100644
--- a/tests/test_knn_interp.py
+++ b/tests/test_knn_interp.py
@@ -165,9 +165,10 @@ def test_disk():
     knn.initialize(star_list)
     knn.solve(star_list)
     knn_file = os.path.join('output','knn_interp.fits')
-    with fitsio.FITS(knn_file,'rw',clobber=True) as f:
-        knn.write(f, 'knn')
-        knn2 = piff.KNNInterp.read(f, 'knn')
+    with piff.writers.FitsWriter.open(knn_file) as w:
+        knn.write(w, 'knn')
+    with piff.readers.FitsReader.open(knn_file) as r:
+        knn2 = piff.KNNInterp.read(r, 'knn')
     np.testing.assert_array_equal(knn.locations, knn2.locations)
     np.testing.assert_array_equal(knn.targets, knn2.targets)
     np.testing.assert_array_equal(knn.kwargs['keys'], knn2.kwargs['keys'])
@@ -246,9 +247,10 @@ def test_decam_disk():
     knn.misalign_wavefront(misalignment)
 
     knn_file = os.path.join('output','decam_wavefront.fits')
-    with fitsio.FITS(knn_file,'rw',clobber=True) as f:
-        knn.write(f, 'decam_wavefront')
-        knn2 = piff.des.DECamWavefront.read(f, 'decam_wavefront')
+    with piff.writers.FitsWriter.open(knn_file) as w:
+        knn.write(w, 'decam_wavefront')
+    with piff.readers.FitsReader.open(knn_file) as r:
+        knn2 = piff.des.DECamWavefront.read(r, 'decam_wavefront')
     np.testing.assert_array_equal(knn.locations, knn2.locations)
     np.testing.assert_array_equal(knn.targets, knn2.targets)
     np.testing.assert_array_equal(knn.keys, knn2.keys)
diff --git a/tests/test_optics.py b/tests/test_optics.py
index 87841909..263cdcdd 100644
--- a/tests/test_optics.py
+++ b/tests/test_optics.py
@@ -366,9 +366,10 @@ def test_disk():
     # save and load
     model = piff.Optical(lam=700.0, template='des')
     model_file = os.path.join('output','optics.fits')
-    with fitsio.FITS(model_file, 'rw', clobber=True) as f:
-        model.write(f, 'optics')
-        model2 = piff.Optical.read(f, 'optics')
+    with piff.writers.FitsWriter.open(model_file) as w:
+        model.write(w, 'optics')
+    with piff.readers.FitsReader.open(model_file) as r:
+        model2 = piff.Optical.read(r, 'optics')
 
     for key in model.kwargs:
         assert key in model2.kwargs, 'key %r missing from model2 kwargs'%key
diff --git a/tests/test_outliers.py b/tests/test_outliers.py
index 2366e7d7..11ffb4ae 100644
--- a/tests/test_outliers.py
+++ b/tests/test_outliers.py
@@ -159,11 +159,11 @@ def test_base():
     # Mock this by pretending that MADOutliers is the only subclass of Outliers.
     filename = os.path.join('input','D00240560_r_c01_r2362p01_piff.fits')
     with mock.patch('piff.Outliers.valid_outliers_types', {'MAD': piff.outliers.MADOutliers}):
-        with fitsio.FITS(filename,'r') as f:
-            np.testing.assert_raises(ValueError, piff.Outliers.read, f, extname='psf_outliers')
+        with piff.readers.FitsReader.open(filename) as r:
+            np.testing.assert_raises(ValueError, piff.Outliers.read, r, 'psf_outliers')
 
-    with fitsio.FITS(filename,'r') as f:
-        out = piff.Outliers.read(f, extname='psf_outliers')
+    with piff.readers.FitsReader.open(filename) as r:
+        out = piff.Outliers.read(r, 'psf_outliers')
         print(out)
 
     # Check that registering new types works correctly
diff --git a/tests/test_pixel.py b/tests/test_pixel.py
index 0c7f3f51..87842b54 100644
--- a/tests/test_pixel.py
+++ b/tests/test_pixel.py
@@ -359,9 +359,9 @@ def test_basis_interp():
     np.testing.assert_raises(RuntimeError, basis.solve, [star])
     np.testing.assert_raises(RuntimeError, basis.interpolate, star)
     file_name = os.path.join('output','test_basis_interp.fits')
-    with fitsio.FITS(file_name,'rw',clobber=True) as fout:
+    with piff.writers.FitsWriter.open(file_name) as w:
         with np.testing.assert_raises(RuntimeError):
-            basis.write(fout, extname='basis')
+            basis.write(w, 'basis')
 
     # Other options for order
     basis = piff.BasisPolynomial(order=2)
diff --git a/tests/test_poly.py b/tests/test_poly.py
index 7dfcec17..9925b3f4 100644
--- a/tests/test_poly.py
+++ b/tests/test_poly.py
@@ -339,10 +339,10 @@ def quadratic_func(pos):
     extname = "interp"
     dirname = 'output'
     filename=os.path.join(dirname, fname)
-    with fitsio.FITS(filename,'rw',clobber=True) as f:
-        interp.write(f, extname=extname)
-    with fitsio.FITS(filename, "r") as f2:
-        interp2 = piff.Polynomial.read(f2, extname=extname)
+    with piff.writers.FitsWriter.open(filename) as w:
+        interp.write(w, extname)
+    with piff.readers.FitsReader.open(filename) as r:
+        interp2 = piff.Polynomial.read(r, extname)
 
     # The type and other parameters should now have been overwritten and updated
     assert interp2.poly_type == interp.poly_type
@@ -396,9 +396,9 @@ def test_poly_raise():
 
     # Cannot write before running fit.
     filename = 'output/test_invalid.fits'
-    with fitsio.FITS(filename,'rw',clobber=True) as f:
+    with piff.writers.FitsWriter.open(filename) as w:
         with np.testing.assert_raises(RuntimeError):
-            interp.write(f, extname='junk')
+            interp.write(w, 'junk')
 
 
 @timer
diff --git a/tests/test_simple.py b/tests/test_simple.py
index 1ef3c4ff..171939f6 100644
--- a/tests/test_simple.py
+++ b/tests/test_simple.py
@@ -665,12 +665,12 @@ def test_model():
     # Mock this by pretending that Gaussian is the only subclass of Model.
     filename = os.path.join('input','D00240560_r_c01_r2362p01_piff.fits')
     with mock.patch('piff.Model.valid_model_types', {'Gaussian': piff.Gaussian}):
-        with fitsio.FITS(filename,'r') as f:
-            np.testing.assert_raises(ValueError, piff.Model.read, f, extname='psf_model')
+        with piff.readers.FitsReader.open(filename) as r:
+            np.testing.assert_raises(ValueError, piff.Model.read, r, 'psf_model')
 
     # But normally this file can be read...
-    with fitsio.FITS(filename,'r') as f:
-        model = piff.Model.read(f, extname='psf_model')
+    with piff.readers.FitsReader.open(filename) as r:
+        model = piff.Model.read(r, 'psf_model')
         print(model)
 
     # We don't have any abstract model base classes (as we do for interp for instance),
@@ -713,21 +713,21 @@ def test_interp():
     np.testing.assert_raises(NotImplementedError, interp.interpolateList, [None])
 
     filename1 = os.path.join('output','test_interp.fits')
-    with fitsio.FITS(filename1,'rw',clobber=True) as f:
-        np.testing.assert_raises(NotImplementedError, interp._finish_write, f, extname='interp')
+    with piff.writers.FitsWriter.open(filename1) as w:
+        np.testing.assert_raises(NotImplementedError, interp._finish_write, w.nested('interp'))
     filename2 = os.path.join('input','D00240560_r_c01_r2362p01_piff.fits')
-    with fitsio.FITS(filename2,'r') as f:
-        np.testing.assert_raises(NotImplementedError, interp._finish_read, f, extname='interp')
+    with piff.readers.FitsReader.open(filename2) as r:
+        np.testing.assert_raises(NotImplementedError, interp._finish_read, r.nested('interp'))
 
     # Invalid to read a type that isn't a piff.Interp type.
     # Mock this by pretending that Mean is the only subclass of Interp.
     with mock.patch('piff.Interp.valid_interp_types', {'Mean': piff.Mean}):
-        with fitsio.FITS(filename2,'r') as f:
-            np.testing.assert_raises(ValueError, piff.Interp.read, f, extname='psf_interp')
+        with piff.readers.FitsReader.open(filename2) as r:
+            np.testing.assert_raises(ValueError, piff.Interp.read, r, 'psf_interp')
 
     # But normally this file can be read...
-    with fitsio.FITS(filename2,'r') as f:
-        interp = piff.Interp.read(f, extname='psf_interp')
+    with piff.readers.FitsReader.open(filename2) as r:
+        interp = piff.Interp.read(r, 'psf_interp')
         print(interp)
 
     # BasisInterp is an abstract base class, so it shouldn't be in the list of valid types
diff --git a/tests/test_star.py b/tests/test_star.py
index 20c382dd..7b22c31e 100644
--- a/tests/test_star.py
+++ b/tests/test_star.py
@@ -478,12 +478,12 @@ def test_io():
 
         file_name = os.path.join('output','star_io.fits')
         print('Writing stars to ',file_name)
-        with fitsio.FITS(file_name,'rw',clobber=True) as fout:
-            piff.Star.write(stars, fout, extname='stars')
+        with piff.writers.FitsWriter.open(file_name) as w:
+            piff.Star.write(stars, w, 'stars')
 
         print('Reading from ',file_name)
-        with fitsio.FITS(file_name,'r') as fin:
-            stars2 = piff.Star.read(fin, extname='stars')
+        with piff.readers.FitsReader.open(file_name) as r:
+            stars2 = piff.Star.read(r, 'stars')
 
         for s1, s2 in zip(stars,stars2):
             assert s1.data['x'] == s2.data['x']
@@ -560,12 +560,12 @@ def test_multifit_io():
 
     file_name = os.path.join('output','star_multifit_io.fits')
     print('Writing stars to ',file_name)
-    with fitsio.FITS(file_name,'rw',clobber=True) as fout:
-        piff.Star.write(stars, fout, extname='stars')
+    with piff.writers.FitsWriter.open(file_name) as w:
+        piff.Star.write(stars, w, 'stars')
 
     print('Reading from ',file_name)
-    with fitsio.FITS(file_name,'r') as fin:
-        stars2 = piff.Star.read(fin, extname='stars')
+    with piff.readers.FitsReader.open(file_name) as r:
+        stars2 = piff.Star.read(r, 'stars')
 
     for s1, s2 in zip(stars,stars2):
         assert s1.data['x'] == s2.data['x']