Last active
January 23, 2020 15:38
-
-
Save tomr-stargazer/7a4efed0f8a1f5d410066ca87515f692 to your computer and use it in GitHub Desktop.
Output of attempting to load CLASS-generated FITS spectrum into specutils
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ipython --pylab | |
| Python 3.6.0 |Anaconda custom (x86_64)| (default, Dec 23 2016, 13:19:00) | |
| Type "copyright", "credits" or "license" for more information. | |
| IPython 5.1.0 -- An enhanced Interactive Python. | |
| ? -> Introduction and overview of IPython's features. | |
| %quickref -> Quick reference. | |
| help -> Python's own help system. | |
| object? -> Details about 'object', use 'object??' for extra details. | |
| Using matplotlib backend: MacOSX | |
| In [1]: import specutils | |
| INFO:root:Successfully loaded reader "cubetest1". | |
| INFO:root:Successfully loaded reader "tabular-fits". | |
| INFO:root:Successfully loaded reader "wcs1d-fits". | |
| In [3]: s1d = specutils.wcs_fits_reader.wcs1d_fits('test_L1455IRS1_DCN2-1.fits') | |
| INFO:root:Spectrum file looks like wcs1d-fits | |
| WARNING: Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. [astropy.io.fits.header] | |
| WARNING:astropy:Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. | |
| WARNING: FITSFixedWarning: 'celfix' made the change 'Linear transformation matrix is singular'. [astropy.wcs.wcs] | |
| WARNING:astropy:FITSFixedWarning: 'celfix' made the change 'Linear transformation matrix is singular'. | |
| --------------------------------------------------------------------------- | |
| SingularMatrixError Traceback (most recent call last) | |
| <ipython-input-3-165e2b31c6c4> in <module>() | |
| ----> 1 s1d = specutils.wcs_fits_reader.wcs1d_fits('test_L1455IRS1_DCN2-1.fits') | |
| 2 | |
| /Users/tsrice/anaconda3/lib/python3.6/site-packages/specutils-0.1.dev727-py3.6.egg/specutils/io/registers.py in wrapper(*args, **kwargs) | |
| 32 @wraps(func) | |
| 33 def wrapper(*args, **kwargs): | |
| ---> 34 return func(*args, **kwargs) | |
| 35 return wrapper | |
| 36 return decorator | |
| /Users/tsrice/anaconda3/lib/python3.6/site-packages/specutils-0.1.dev727-py3.6.egg/specutils/io/default_loaders/wcs_fits_reader.py in wcs1d_fits(file_name, spectral_axis_unit, **kwargs) | |
| 37 with fits.open(file_name, **kwargs) as hdulist: | |
| 38 header = hdulist[0].header | |
| ---> 39 wcs = WCS(header) | |
| 40 | |
| 41 if 'BUNIT' in header: | |
| /Users/tsrice/anaconda3/lib/python3.6/site-packages/astropy/wcs/wcs.py in __init__(self, header, fobj, key, minerr, relax, naxis, keysel, colsel, fix, translate_units, _do_set) | |
| 501 | |
| 502 if _do_set: | |
| --> 503 self.wcs.set() | |
| 504 | |
| 505 for fd in close_fds: | |
| SingularMatrixError: ERROR 3 in wcsset() at line 2213 of file cextern/wcslib/C/wcs.c: | |
| Linear transformation matrix is singular. | |
| ERROR 3 in linset() at line 607 of file cextern/wcslib/C/lin.c: | |
| PCi_ja matrix is singular. | |
| In [5]: from astropy.io.fits import getdata | |
| In [6]: data, header = getdata("test_L1455IRS1_DCN2-1.fits", header=True) | |
| WARNING: Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. [astropy.io.fits.header] | |
| WARNING:astropy:Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. | |
| In [10]: data | |
| Out[10]: | |
| array([[[[-0.00654329, 0.01778164, -0.01153884, ..., -0.01729367, | |
| 0.01069437, -0.01851475]]]]) | |
| In [11]: data.shape | |
| Out[11]: (1, 1, 1, 37275) | |
| In [7]: header | |
| Out[7]: | |
| SIMPLE = T / | |
| BITPIX = -64 | |
| NAXIS = 4 / | |
| NAXIS1 = 37275 / | |
| NAXIS2 = 1 / | |
| NAXIS3 = 1 / | |
| NAXIS4 = 1 / | |
| BLOCKED = T / | |
| DATAMIN = -0.146224617958E+000 / | |
| DATAMAX = 0.140489089489E+001 / | |
| BUNIT = 'K ' / | |
| CTYPE1 = 'FREQ ' / | |
| CRVAL1 = 0.000000000000E+000 / Offset frequency | |
| CDELT1 = 0.488281250000E+005 / Frequency resolution | |
| CRPIX1 = 0.176140019531E+005 / | |
| CTYPE2 = 'RA---GLS' / | |
| EQUINOX = 0.200000000000E+004 / | |
| CRVAL2 = 0.519129166667E+002 / | |
| CDELT2 = 0.000000000000E+000 / | |
| CRPIX2 = 0.000000000000E+000 / | |
| CTYPE3 = 'DEC--GLS' / | |
| CRVAL3 = 0.302175000000E+002 / | |
| CDELT3 = 0.000000000000E+000 / | |
| CRPIX3 = 0.000000000000E+000 / | |
| CTYPE4 = 'STOKES ' / | |
| CRVAL4 = 1.0000000000000 / | |
| CDELT4 = 0.0000000000000 / | |
| CRPIX4 = 0.0000000000000 / | |
| TELESCOP= '30ME1-LI-F0-' / | |
| OBJECT = 'L1455IRS1 ' / | |
| LINE = 'DCN(2-1) ' / Line name | |
| RESTFREQ= 0.144828001000E+012 / Rest frequency | |
| VELO-LSR= 0.500000000000E+004 / Velocity of reference channel | |
| DELTAV = -0.101073712111E+003 / Velocity spacing of channels | |
| IMAGFREQ= 0.157326768103E+012 / Image frequency | |
| TSYS = 0.845444717407E+002 / System temperature | |
| OBSTIME = 0.165860595703E+004 / Integration time | |
| SCAN-NUM= 0.100000000000E+001 / Scan number | |
| TAU-ATM = 0.288365986198E-001 / Atmospheric opacity | |
| NPHASE = 1 / Number of frequency phases | |
| DELTAF1 = 0.000000000000E+000 / Frequency offset Phase 1 | |
| PTIME1 = 0.498978763819E+000 / Duration of Phase 1 | |
| WEIGHT1 = 0.000000000000E+000 / Weight of Phase 1 | |
| GAINIMAG= 0.501190014184E-001 / Image sideband gain ratio | |
| BEAMEFF = 0.930000007153E+000 / Beam efficiency | |
| FORWEFF = 0.930000007153E+000 / Image sideband gain ratio | |
| ORIGIN = 'CLASS-Grenoble ' / | |
| DATE = '2017-09-21T00:00:00.000' / Date written | |
| TIMESYS = 'UTC ' / | |
| DATE-OBS= '2014-12-26T18:48:39.588' / Date observed | |
| DATE-RED= '2017-09-21T00:00:00.000' / Date reduced | |
| ELEVATIO= 0.592332514636E+002 / Telescope elevation | |
| AZIMUTH = 0.000000000000E+000 / Telescope azimuth | |
| UT = ' 18:48:39.588' / Universal time at start | |
| LST = ' 00:55:48.714' / Sideral time at start |
@keflavich thanks for this - I'm able to load the file in pyspeckit (and get access to the spectral axis via e.g. spectrum.xarr) so this is definitely a direction I can move in. I'll follow examples in the documentation for pyspeckit.Spectrum, e.g. this chunk:
Examples
--------
>>> sp = pyspeckit.Spectrum(data=np.random.randn(100),
xarr=np.linspace(-50, 50, 100), error=np.ones(100)*0.1,
xarrkwargs={'unit':'km/s', 'refX':4.829, 'refX_unit':'GHz',
'xtype':'VLSR-RAD'}, header={})
>>> xarr = pyspeckit.units.SpectroscopicAxis(np.linspace(-50,50,100),
units='km/s', refX=6562.83, refX_unit='angstroms')
>>> data = np.random.randn(100)*5 + np.random.rand(100)*100
>>> err = np.sqrt(data/5.)*5. # Poisson noise
>>> sp = pyspeckit.Spectrum(data=data, error=err, xarr=xarr, header={})
>>> # if you already have a simple fits file
>>> sp = pyspeckit.Spectrum('test.fits')but I'll be following closely in case either spectral-cube or specutils develop ways to deal with this weird FITS file.
Yeah, this one is a real oddball. It would be helpful if you could post the file as an Issue on specutils, and perhaps even help us figure out how to correctly read it.
In pyspeckit, there is also a CLASS file reader:
http://pyspeckit.readthedocs.io/en/latest/classfiles.html
http://pyspeckit.readthedocs.io/en/latest/guide_class.html
@keflavich unfortunately, I keep hitting errors when attempting to load the file (in its ".30m" class format) into pyspeckit, in either "obsblock" or "spectra" mode.
In [9]: thing = class_to_obsblocks("FTS_2014_December.30m", telescope='*-F*', line='HCN-3MM')
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-9-849e89ad1e76> in <module>()
----> 1 thing = class_to_obsblocks("FTS_2014_December.30m", telescope='*-F*', line='HCN-3MM')
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in wrapper(*arg, **kwargs)
43 def wrapper(*arg,**kwargs):
44 t1 = time.time()
---> 45 res = func(*arg,**kwargs)
46 t2 = time.time()
47 log.info('%s took %0.5g s' % (func.__name__, (t2-t1)))
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in class_to_obsblocks(filename, telescope, line, datatuple, source, imagfreq, DEBUG, **kwargs)
1481 """
1482 if datatuple is None:
-> 1483 spectra,header,indexes = read_class(filename,DEBUG=DEBUG, **kwargs)
1484 else:
1485 spectra,header,indexes = datatuple
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in wrapper(*arg, **kwargs)
43 def wrapper(*arg,**kwargs):
44 t1 = time.time()
---> 45 res = func(*arg,**kwargs)
46 t2 = time.time()
47 log.info('%s took %0.5g s' % (func.__name__, (t2-t1)))
TypeError: read_class() got an unexpected keyword argument 'DEBUG'In [12]: spec = class_to_spectra("L1455IRS1_FTS_HCN-3MM.30m")
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
<ipython-input-12-341edc7f7896> in <module>()
----> 1 spec = class_to_spectra("L1455IRS1_FTS_HCN-3MM.30m")
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in wrapper(*arg, **kwargs)
43 def wrapper(*arg,**kwargs):
44 t1 = time.time()
---> 45 res = func(*arg,**kwargs)
46 t2 = time.time()
47 log.info('%s took %0.5g s' % (func.__name__, (t2-t1)))
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in class_to_spectra(filename, datatuple, **kwargs)
1593 """
1594 if datatuple is None:
-> 1595 spectra,header,indexes = read_class(filename, **kwargs)
1596 else:
1597 spectra,header,indexes = datatuple
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in wrapper(*arg, **kwargs)
43 def wrapper(*arg,**kwargs):
44 t1 = time.time()
---> 45 res = func(*arg,**kwargs)
46 t2 = time.time()
47 log.info('%s took %0.5g s' % (func.__name__, (t2-t1)))
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in read_class(filename, downsample_factor, sourcename, telescope, line, posang, verbose, flag_array)
1389 log.info("Reading...")
1390 selection = [ii
-> 1391 for source in sourcename
1392 for tel in telescope
1393 for li in line
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in <listcomp>(.0)
1395 telescope=tel,
1396 line=li,
-> 1397 posang=posang)]
1398
1399 sphdr = classobj.read_observations(selection)
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in select_spectra(self, all, line, linere, linereflags, number, scan, offset, source, sourcere, sourcereflags, range, quality, telescope, telescopere, telescopereflags, subscan, entry, posang, frequency, section, user, include_old_versions)
1314 else True) and
1315 (h['XVER'] > 0 if not include_old_versions else True)
-> 1316 for h in self.headers
1317 ]
1318
/Users/tsrice/anaconda3/lib/python3.6/site-packages/pyspeckit/spectrum/readers/read_class.py in <listcomp>(.0)
1314 else True) and
1315 (h['XVER'] > 0 if not include_old_versions else True)
-> 1316 for h in self.headers
1317 ]
1318
KeyError: 'XVER'I'll pop over to the Issues in specutils and see if it's useful to anyone.
If it's relevant:
In [13]: pyspeckit.__version__
Out[13]: '0.1.20'
Shouldn't CDELT equal to 0 result in a singular matrix?
Yeah, setting CDELT2, CDELT3, and CDELT4 to 1 allows the WCS object to be created.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I don't know exactly what's going on here, but it is really weird to have
CRVAL1 = 0!We might need to special-case CLASS FITS files. Have you tried loading this with
pyspeckit? I spent some time working with CLASS spectra there; you can even load CLASS 1D spectra directly into pyspeckit.