On that last note, I would add at *some point* `astropy.io.fits` is going to start returning some kind of `NDData` object rather than just plain Numpy arrays.  It will probably be a subclass specially designed to meet the special needs of FITS.  But part of this interface might include ways to provide mappings from FITS headers to data models accessible through the `NDData` object's meta attribute.  I hesitate to bring this up because I can't make any promises right now as to when that would happen, but it's just some kind of big picture stuff to keep in mind.