Writing Quantity to FITS file

[astrofrog]

I often encounter the following issue:

In [1]: from astropy import units as u                                                                                                                 

In [2]: from astropy.io import fits                                                                                                                    

In [3]: fits.writeto('test.fits', [[0]] * u.Jy)                                                                                                        
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-3-51c0417c4f82> in <module>
----> 1 fits.writeto('test.fits', [[0]] * u.Jy)

~/python/dev/lib/python3.7/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    519                             kwargs[new_name[i]] = value
    520 
--> 521             return function(*args, **kwargs)
    522 
    523         return wrapper

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/convenience.py in writeto(filename, data, header, output_verify, overwrite, checksum)
    418         hdu = PrimaryHDU(data, header=header)
    419     hdu.writeto(filename, overwrite=overwrite, output_verify=output_verify,
--> 420                 checksum=checksum)
    421 
    422 

~/python/dev/lib/python3.7/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    519                             kwargs[new_name[i]] = value
    520 
--> 521             return function(*args, **kwargs)
    522 
    523         return wrapper

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py in writeto(self, name, output_verify, overwrite, checksum)
    371         hdulist = HDUList([self])
    372         hdulist.writeto(name, output_verify, overwrite=overwrite,
--> 373                         checksum=checksum)
    374 
    375     @classmethod

~/python/dev/lib/python3.7/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    519                             kwargs[new_name[i]] = value
    520 
--> 521             return function(*args, **kwargs)
    522 
    523         return wrapper

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py in writeto(self, fileobj, output_verify, overwrite, checksum)
    927             for hdu in self:
    928                 hdu._prewriteto(checksum=checksum)
--> 929                 hdu._writeto(hdulist._file)
    930                 hdu._postwriteto()
    931         hdulist.close(output_verify=output_verify, closed=closed)

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py in _writeto(self, fileobj, inplace, copy)
    673 
    674         with _free_space_check(self, dirname):
--> 675             self._writeto_internal(fileobj, inplace, copy)
    676 
    677     def _writeto_internal(self, fileobj, inplace, copy):

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py in _writeto_internal(self, fileobj, inplace, copy)
    679         if not inplace or self._new:
    680             header_offset, _ = self._writeheader(fileobj)
--> 681             data_offset, data_size = self._writedata(fileobj)
    682 
    683             # Set the various data location attributes on newly-written HDUs

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/base.py in _writedata(self, fileobj)
    609         if self._data_loaded or self._data_needs_rescale:
    610             if self.data is not None:
--> 611                 size += self._writedata_internal(fileobj)
    612             # pad the FITS data block
    613             if size > 0:

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/hdu/image.py in _writedata_internal(self, fileobj)
    626                         output.byteswap(True)
    627                         try:
--> 628                             fileobj.writearray(output)
    629                         finally:
    630                             output.byteswap(True)

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/file.py in writearray(self, array)
    363 
    364         if hasattr(self._file, 'write'):
--> 365             _array_to_file(array, self._file)
    366 
    367     def flush(self):

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/util.py in _array_to_file(arr, outfile)
    643     else:
    644         # Just pass the whole array to the write routine
--> 645         return write(arr, outfile)
    646 
    647     # Write one chunk at a time for systems whose fwrite chokes on large

~/python/dev/lib/python3.7/site-packages/astropy/io/fits/util.py in <lambda>(a, f)
    621 
    622     if isfile(outfile) and not isinstance(outfile, io.BufferedIOBase):
--> 623         write = lambda a, f: a.tofile(f)
    624     else:
    625         write = _array_to_file_like

~/python/dev/lib/python3.7/site-packages/astropy/units/quantity.py in tofile(self, fid, sep, format)
   1388 
   1389     def tofile(self, fid, sep="", format="%s"):
-> 1390         raise NotImplementedError("cannot write Quantities to file.  Write "
   1391                                   "array with q.value.tofile(...)")
   1392 

NotImplementedError: cannot write Quantities to file.  Write array with q.value.tofile(...)

and I was wondering whether it might make sense to actually allow this to work, and automatically set BUNIT in the FITS header based on the quantity's unit? (emitting a warning if header already defines a unit that gets overwritten).

[mhvk]

Trying to trace this a bit, I see that hdu.data becomes a Quantity, and that cannot be written directly.

I think it would not be too difficult to solve: just adjust data.setter in hdu.image._ImageBaseHDU to check for a unit on the input and adjust the header accordingly. One might do this very bluntly, by passing all possible input through u.Quantity and appropriately deal with dimensionless.

Alternatively, since hdu.data is already a possibly scaled version of the actual data stored in the fits file, one could also go even further, and generally attach the unit to it as well, ie., change it to a Quantity always. Though I think this might have problems with backwards compatibility...

[MSeifert04]

Is there any reason that you don't use NDData or one derivative that supports writing/reading FITS files with units?

I'm not exactly sure if we should support quantities (or units) in io.fits. I mostly understood io.fits as low-level wrapper for the fits format, not as the high-level "astropy interface", which would be astropy.nddata or astropy.table. These have the advantage that these classes can also be used for other file-formats, not just FITS.

[astrofrog]

@MSeifert04 - nothing beats fits.writeto when trying to write an array quickly to disk 😆 - I was just thinking that since fits.writeto is a convenience anyway, it could silently do the right thing. No strong opinion, but if we decide to not support it it might be good to have a clearer error message.

[pllim]

FITS is already complicated enough as is. If a higher level API (like Table or NDData) can do this, I would vote that we should not touch FITS.

[astrofrog]

Unfortunately we don't actually have FITS readers/writers for plain NDData objects (yet), but we really should add some!

Another thing we could do here that would work in general would be in fits.writeto to check for __array__ on objects that are passed in, since using that for a quantity would automatically drop the unit, and it would be good to support objects that implement __array__ in general?

[MSeifert04]

Unfortunately we don't actually have FITS readers/writers for plain NDData objects (yet), but we really should add some!

Ah yeah, they are only implemented on CCDData... right.

Another thing we could do here that would work in general would be in fits.writeto to check for `` on objects that are passed in, since using that for a quantity would automatically drop the unit, an__array__d it would be good to support objects that implement __array__ in general?

I think that we duck-type on the tofile method of the object passed in. Not sure what would regress if we changed that to the __array__ interface.

[MSeifert04]

Hm, I had a PR (#4799) trying to add some reader/writer for FITS.

However it probably should work if we associate the CCDData reader/writer with the NDIOMixin. 😄

[mhvk]

I agree with fits being for low-level stuff (so perhaps scrap my suggestion to change the data setter), but also that fits.writeto is a convenience method which might as well support this. It would seem to make sense to use NDData internally to make that work.

[saimn]

I think it would make sense to be able to write directly a quantity object with fits.writeto, and to read it maybe with a new keyword in fits.open. But if there is a way to achieve this with NDData or CCDData that probably better.