search

astrostat / pylira

A Python package for Bayesian low-counts image reconstruction and analysis
9 stars 7 forks source link

Which data structure with WCS support to use? #22

Open adonath opened 2 years ago

adonath commented 2 years ago

Typically the astronomical data to be used with pylira is given as FITS files with associated WCS information. Currently pylira does not support handling any WCS information for I/O or plotting. There are already a few existing data structures one could use:

vkashyap commented 2 years ago

Why not steal the header from the input image and paste it on to the outputs? They have the same sizes and orientation.

On Mon, Nov 15, 2021 at 11:51 AM Axel Donath @.***> wrote:

Typically the astronomical data to be used with pylira is given as FITS files with associated WCS information. Currently pylira does not support handling any WCS information for I/O or plotting. There are already a few existing data structures one could use:

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/astrostat/pylira/issues/22, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAOWYQM6MZJR7Z35NSUG4A3UME27PANCNFSM5ICE7ZZA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

adonath commented 2 years ago

Yes, I agree this works as well (basically option 3). So when working with pylira, users would write:

from astropy.io import fits
from pylira import LIRADeconvolver

hdulist = fits.open("counts.fits")
counts_hdu = hdulist[0]

# read psf, exposure and baseline here...

deconvolve = LIRADeconvolver()

result = deconvolve.run({"counts": counts_hdu.data, ...})

hdu_out = fits.PrimaryHDU(data=result, header=counts_hdu.header)

hdulist = fits.HDUList([hdu_out])
hdulist.writeto("result.fits")

Which is acceptable, but maybe a bit inconvenient.

But I agree LIRA itself does not need to know anything about the WCS it is only needed for I/O and visualization and as long as we provide an interface based on numpy arrays users can do whatever they want. E.g. for Gammapy I'd use:

from gammapy.maps import Map
from pylira import LIRADeconvolver

counts = Map.read("counts.fits")

# read psf, exposure and baseline here...

deconvolve = LIRADeconvolver()

data = deconvolve.run({"counts": counts.data, ...})

result = counts.copy(data=data)
result.write("result.fits")
adonath commented 2 years ago

For now we can do the following: handle all the data passed from and to LIRA as a dict of numpy arrays and allow for an optional astropy.wcs.WCS object, which can then optionally be used for plotting and serialization.