Repository navigation
wcs _TAB option: example please #7883
Description
Activity
Hello. Do you not have a FITS file already? If not, what format is your data in? What instrument did you use to get the data? Usually, an instrument comes with its own calibration software.
Also, unless you are stuck with Python 2, we recommend using the latest release of 3.x series.
- I’m in the process of commissioning a new mode on the Robert Stobie Spectrograph of the Southern African Large Telescope (I did the spectrograph optical design). Right now I am writing the calibration/ analysis software for this mode. I need to define a fits file data product that deals with the non-linear prism dispersion in this mode, and it seems that the wcs -TAB option is most appropriate. So I am trying to put data that has already been processed from raw detector fits files into reduced fits files that incorporate wavelength calibration of the spectra. All my attempts to do this result in garbage. I thought I could find dumb mistakes or version problems by comparing with a –TAB that does work, hence the title. Failing this, I show a failed example here. Below is the info() and header for a test fits file I made. The current SALT pipeline uses python 2.7.3. I would definitely like to avoid updating the whole pipeline!>> hdul.info()Filename: testfits.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 167 () 1 SCI 1 ImageHDU 24 (1601, 514, 2) float32 2 VAR 1 ImageHDU 9 (1601, 514, 2) float32 3 BPM 1 ImageHDU 9 (1601, 514, 2) uint8 4 TMAP 1 ImageHDU 9 (1601, 514, 2) uint8 5 TGT 1 BinTableHDU 43 108R x 16C [19A, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D] 6 WCS-table 1 BinTableHDU 13 12R x 2C [K, D]>> print repr(hdul['SCI'].header)XTENSION= 'IMAGE ' / Image extension BITPIX = -32 / array data type NAXIS = 3 / number of array dimensions NAXIS1 = 1601 NAXIS2 = 514 NAXIS3 = 2 PCOUNT = 0 / number of parameters GCOUNT = 1 / number of groups EXTNAME = 'SCI ' / Extension name EXTVER = 1 / Extension number VAREXT = 2 / Variance frame extension BPMEXT = 3 / BPM Extension CRPIX1 = 0 / WCS: X reference pixel CRPIX2 = 0 / WCS: Y reference pixel CRVAL1 = 0 / WCS: X reference coordinate value CRVAL2 = 4.0 / WCS: Y reference coordinate value CDELT1 = 1 / WCS: X pixel size CDELT2 = 4.0 / WCS: Y pixel size CTYPE1 = 'WAVE-TAB' / X type CTYPE2 = 'pixel ' / Y type CUNIT1 = 'Angstrom' PS1_0 = 'WCS-table' PS1_1 = 'WaveCoord' PS1_2 = 'WaveIndex'>> hdul['WCS-table'].dataFITS_rec([( 0, 3200.), ( 1, 3225.), ( 2, 3250.), ( 3, 3275.), ( 4, 3300.), ( 5, 3325.), ( 6, 3350.), ( 7, 3375.), ( 8, 3400.), ( 9, 3435.), (10, 3470.), (11, 3505.)], dtype=(numpy.record, [('WaveIndex', '>i8'), ('WaveCoord', '>f8')])) If I do a specWcs = wcs.WCS(hdul['SCI'],hdul).sub(1) print specWcs gives: WCS Keywords Number of WCS axes: 1 CTYPE : 'WAVE-TAB' CRVAL : 0.0 CRPIX : 0.0 PC1_1 : 1.0 CDELT : 1.0 NAXIS : 1601 514 2 print specWcs.all_pix2world(np.arange(12),0) gives the bogus answer: [array([ 6.90074734e-310, 1.59840355e-316, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000])] print specWcs.wcs_pix2world(np.arange(12),0) similarly gives: [array([ 6.90074734e-310, 1.65180137e-316, 1.84706718e-316, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000])] print specWcs.p4_pix2foc(np.arange(12),0) gives the proper answer: [array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])] Clearly the second step, the WAVE-TAB step, is getting bogus information somewhere. Thanks! From: P. L. Lim [mailto:[email protected]] Sent: Wednesday, October 17, 2018 8:18 AM To: astropy/astropy Cc: KENNETH H NORDSIECK; Author Subject: Re: [astropy/astropy] wcs _TAB option: example please (#7883) Hello. Do you not have a FITS file already? If not, what format is your data in? What instrument did you use to get the data? Usually, an instrument comes with its own calibration software. Also, unless you are stuck with Python 2, we recommend using the latest release of 3.x series. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#7883 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ADRf7VWScPUk6pmZlqyfZ8D3xtYOyjazks5uly4HgaJpZM4XYSx1>.
@knordsieck There's no support for the -TAB convention in astropy yet. There's some initial scoping of work being done to see what needs to be implemented.
The lookup table distortion is a different structure - it's the one described in the distortion paper. It's essentially an image extension with an optional interpolation applied.
- Aha! I was assuming the Greisen (2006) definition of the –TAB option was just an update of the Calabretta (2004) ( “distortion paper”) table lookup option described in section 3.4. Apparently not. Perhaps a little documentation update would clarify this .. So astropy does support the Calabretta version? Do you have a working example of a fits file using it? From: Nadia Dencheva [mailto:[email protected]] Sent: Wednesday, October 17, 2018 11:55 AM To: astropy/astropy Cc: KENNETH H NORDSIECK; Mention Subject: Re: [astropy/astropy] wcs _TAB option: example please (#7883) @knordsieck<https://github.com/knordsieck> There's no support for the -TAB convention in astropy yet. There's some initial scoping of work being done to see what needs to be implemented. The lookup table distortion is a different structure - it's the on described in the distortion paper<https://fits.gsfc.nasa.gov/wcs/dcs_20040422.pdf>. It's essentially an image extension with an optional interpolation applied. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#7883 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ADRf7X-ibxY1aSMb5hZODodSIYT_fbgLks5ul2DJgaJpZM4XYSx1>.
Since you mentioned SALT and spectroscopy, cc @crawfordsm
astropy implements the
Lookuptable option from the distortion paper. An example file is in the astropy source directory, astropy/wcs/tests/data/dist_lookup.fits.gzIt is used in HST imaging modes to describe part of the distortion in ACS and WFC3. Another source of files like that is the STScI archive.
- Hmm. This will take some study. The example you give is pretty complex. Is there a further writeup on it, and is there a couple of astropy commands I can try on it to verify that my python and astropy versions are adequate? From: Nadia Dencheva [mailto:[email protected]] Sent: Wednesday, October 17, 2018 1:58 PM To: astropy/astropy Cc: KENNETH H NORDSIECK; Mention Subject: Re: [astropy/astropy] wcs _TAB option: example please (#7883) astropy implements the Lookup table option from the distortion paper. An example file is in the astropy source directory, astropy/wcs/tests/data/dist_lookup.fits.gz<https://github.com/astropy/astropy/blob/master/astropy/wcs/tests/data/dist_lookup.fits.gz> It is used in HST imaging modes to describe part of the distortion in ACS and WFC3. Another source of files like that is the STScI archive. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#7883 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ADRf7X2z8IcSGOgYHZ8_pp8BHuWdXe20ks5ul321gaJpZM4XYSx1>.
Yeah, it's not straightforward. In astropy, using the file from the tests directory (I used astropy master but this should work with the release as well):
f = fits.open('dist_lookup.fits.gz') w = wcs.WCS(f[1].header, f) w.p4_pix2foc(1,1,1) Out[5]: [array(0.96995641), array(1.01104154)]p4_pix2focevaluates only the lookup table distortion.f.info() Filename: dist_lookup.fits.gz No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 4 () 1 SCI 1 ImageHDU 170 (100, 100) float32 2 D2IMARR 1 ImageHDU 15 (4096, 1) float32 3 WCSDVARR 1 ImageHDU 15 (65, 33) float32 4 WCSDVARR 2 ImageHDU 15 (65, 33) float32Lookup table extensions have
EXTNAME=WCSDVARR. This file has two Lookup table extensions - one for each axis. Their size is smaller than the image inSCI(65, 33). A correction value is computed by bilinear interpolation. The best explanation for the format is the distortion paper. The paper explains what the keywords in the SCI header which point to the WCSDVARR extensions are.
There are more details to all this...Again, in HST these are used for correcting residual distortions in images.
This is different from the -TAB convention which is stored in a binary table, uses different keywords to connect the SCI extension with the lookup table extension and is approved convention. Hope some of this is helpful.- OK, I get the same answer with my current versions, so that’s ok. If the distortion lookup table is done in the first, p4 step, I will have to deal with the detector binning outside of the wcs. Not ideal, but it might work. I’ll leave this open until I have got farther, to keep you in the string in case of more questions. From: Nadia Dencheva [mailto:[email protected]] Sent: Wednesday, October 17, 2018 4:43 PM To: astropy/astropy Cc: KENNETH H NORDSIECK; Mention Subject: Re: [astropy/astropy] wcs _TAB option: example please (#7883) Yeah, it's not straightforward. In astropy, using the file from the tests directory (I used astropy master but this should work with the release as well): f = fits.open('dist_lookup.fits.gz') w = wcs.WCS(f[1].header, f) w.p4_pix2foc(1,1,1) Out[5]: [array(0.96995641), array(1.01104154)] p4_pix2foc evaluates only the lookup table distortion. f.info() Filename: dist_lookup.fits.gz No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 4 () 1 SCI 1 ImageHDU 170 (100, 100) float32 2 D2IMARR 1 ImageHDU 15 (4096, 1) float32 3 WCSDVARR 1 ImageHDU 15 (65, 33) float32 4 WCSDVARR 2 ImageHDU 15 (65, 33) float32 Lookup table extensions have EXTNAME=WCSDVARR. This file has two Lookup table extensions - one for each axis. Their size is smaller than the image in SCI (65, 33). A correction value is computed by bilinear interpolation. The best explanation for the format is the distortion paper. The paper explains what the keywords in the SCI header which point to the WCSDVARR extensions are. There are more details to all this... Again, in HST these are used for correcting residual distortions in images. This is different from the -TAB convention which is stored in a binary table, uses different keywords to connect the SCI extension with the lookup table extension and is approved convention. Hope some of this is helpful. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#7883 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ADRf7dxLYXyQwdrVWwvCO5PbK3Y4Ipt-ks5ul6RSgaJpZM4XYSx1>.
Joining onto this after trying
specutilsfor this (see: astropy/specutils#528 ): I too would like to be able to write a-TABspectral WCS, per C&G 2006.@nden says above that there is no support for the
-TABconvention, but later says "astropy implements the Lookup table option from the distortion paper." I thought that's what the spectral-TABformat was?This is my specutils example for how I was going to write a per-pixel table spectral WCS (which doesn't work correctly in specutils currently). If there is a way to do this in astropy, I would be interested to see it.
import astropy.units as u import specutils n_pixels=1000 wavelength = np.linspace(300, 1100, n_pixels)*u.nm spectrum = np.random.random(n_pixels) * u.count spec=specutils.Spectrum1D(flux=spectrum, spectral_axis=wavelength) spec.write('outputfile.fits')
@parejkoj The -TAB convention is described in WCS Paper III (the paper for spectral coordinates). The LookupTable extension is a simpler convention described in the unpublished Distortion WCS paper. The latter is implemented in astropy.wcs and wcslib, the former is being worked on now.
What you are attempting to do can be done using gwcs:
import astropy.units as u import gwcs from astropy.modeling.tabular import Tabular1D n_pixels=1000 wavelength = np.linspace(300, 1100, n_pixels)*u.nm forward = Tabular1D(np.arange(n_pixels), lookup_table=wavelength, name='dispersion') pixel_frame = gwcs.CoordinateFrame(naxes=1, axes_type=('detector',), axes_order=(0,)) spectral_frame = gwcs.SpectralFrame(axes_order=(0,), unit=wavelength.unit) pipeline = [(pixel_frame, forward), (spectral_frame, None)] wcs = gwcs.wcs.WCS(pipeline)Of course, this only creates the WCS object. And it won't save it to fits. Are you sure specutils can save a gwcs object to fits file?
No, if I pass the gwcs object to specutils, it writes it out in the same way: the examples given on the above specutils issue behave the same whether you pass in a gwcs object, or the table as the spectral axis directly.
Is there an Issue or PR for the
-TABthat you say "is being worked on now"? Is there a way I could contribute to that?In the jwst pipeline we save this as pure asdf or as AsdfInFits (the asdf part is saved as a FITS extension). I don't know what specutils does but I doubt it can save to -TAB, that's just not implemented in astropy.
The -TAB in astropy.wcs involves wrapping parts of wcslib. We are only now starting to look at what needs to be done and how feasible it is.
I have been unable to get the _TAB transformation class ("table lookup distortions") to work. I need it because of quite nonlinear spectral distortions I need to manage in objective prism images. Could someone point me to a successful use of this class? An actual fits file with the correct header keywords and table extension would be very valuable.
I am currently using astropy 2.0.2. I s this a problem?
Thanks