search

esheldon / fitsio

A python package for FITS input/output wrapping cfitsio
GNU General Public License v2.0
133 stars 57 forks source link

Slicing without scaling #262

Open at88mph opened 5 years ago

at88mph commented 5 years ago

I have a use case where I'm trying to perform some slicing:

fits = fitsio.FITS('myfile.fits', 'r')
image_hdu = fits[1]
cutout = image_hdu[10:60]

The myfile.fits contains BZERO/BSCALE header cards. I would like to slice on the pixels without having them scaled. Can that be done?

I think cfitsio supports it but I'm not sure if I can go that low level through fitsio. Is there a solution with numpy's memmap?

Thank you!

at88mph commented 5 years ago

I should clarify. I know CFITSIO has the ability to pull the underlying pixels unmutated, but I don't think it's exposed via fitsio since as soon as I perform any kind of read operation the pixels are scaled.

Can someone confirm this? Thank you.

beckermr commented 5 years ago

I changed the issue title since it confused me. :)

FWIW, as far as I know, this package aims to be compliant with the FITS standard and uses cfitsio internally, but it does not necessarily aim to implement all of the features in cfitsio. I know there is some subtlety in that statement, but I think it is a reasonable point of view.

Do you have a suggested API for this?

at88mph commented 5 years ago

Thank you for the reply.

Intuitively, I wanted to run it like: fitsio.read('myfile.fits', 'r')[:,30:65]

But that read does what the current one does already, so another argument would probably be needed to make use of the fits_set_bscale and fits_set_bscale routines as per the documentation at https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node26.html.

Something like:

fits = fitsio.FITS('myfile.fits', 'r', ignore_scaling=True)
hdu = fits[0]
header = hdu.read_header()  # BZERO and BSCALE are preserved but ignored.
cutout = hdu[:,30:65]

According to https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node76.html#ffpscl, it those cfitsio routines do not need to get set once and then set back again; it only affects the data read in.

esheldon commented 5 years ago

That sounds fine to me. Would you be willing to submit a pull request?

Note that one may not want this option for all hdus in a file, so it may be we should have this as an option in the read method as well, to over-ride whatever is sent to the constructor.

beckermr commented 5 years ago

According to https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node76.html#ffpscl, it those cfitsio routines do not need to get set once and then set back again; it only affects the data read in.

I read the docs there but did not see any info about how this effects subsequent reads of other extensions. Can you elaborate a bit?

beckermr commented 5 years ago

My read is that we would either have to force the whole file to be read this way, or we need to set and then unset these parameters when reading specific extensions.

beckermr commented 5 years ago

We might be better off parsing the header values ourselves and adjusting the output by hand.

at88mph commented 5 years ago

What if we just exposed the fits_set_bscale and fits_set_tscale routines and left it to the caller to set them as they see fit?

Index #1 at the https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node76.html#ffpscl documentation says that it affects reading/writing from/to the FITS file, which leads me to believe that it's an all-in operation. I'm not sure how cfitsio holds state while it's working with a FITS file such that one could set it on and off arbitrarily while working with different extensions. Would there be a use case for not scaling one extension while scaling another from within a file?

I'm happy to provide a PR once a design is agreed upon.

at88mph commented 5 years ago

I've begun a Work In Progress PR. The API I went with is something like:

fits = fitsio.FITS('myMEFfile.fits', 'r')
hdu1 = fits[1]
hdu1.ignore_scaling = True
header = hdu1.read_header()  # BZERO and BSCALE are preserved but ignored.
cutout = hdu1[:,30:65]
...
hdu3 = fits[3]
print(hdu3.ignore_scaling). # Prints False
...

The default setting is False, but with the option to ignore it per HDU. Am I going the correct direction?

I have NOT added it to the write functions, yet.

esheldon commented 3 years ago

Sorry we let this one slip for so long

I think we'll take what you have with the per-hdu settings. We can decide later if we want to also have a constructor keyword, if it ever comes up. Please do make a test or a few for this as well.

Please do add