Enables I/O of covariance/correlation matrices for the tabular fits format

[kbwestfall]

Addresses #1149 .

This requires the new Covariance subclass of astropy.nddata.NDUncertainty; see https://github.com/astropy/astropy/pull/16690. When on the covariance branch in astropy, this is functional code, but I have punted on some functionality, as you'll see with a number of NotImplementedError exceptions.

I'm opening this as a draft while https://github.com/astropy/astropy/pull/16690 is under review.