Skip to content

Commit

Permalink
Fix problems in the construction of the merged WCS in cube
Browse files Browse the repository at this point in the history
  • Loading branch information
@sergiopasra
sergiopasra committed Aug 8, 2020
1 parent ed43e06 commit 4903c77
Show file tree
Hide file tree
Showing 3 changed files with 146 additions and 64 deletions.
96 changes: 62 additions & 34 deletions megaradrp/processing/cube.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,9 +336,9 @@ def create_cube_from_rss(rss, p=1, target_scale_arcsec=1.0, conserve_flux=True):
cube[0].data = np.moveaxis(cube_data, 2, 0)
cube[0].data.astype('float32')

# Merge headers
merge_wcs(rss['FIBERS'].header, rss[0].header, out=cube[0].header)
# Update values of WCS
sky_header = rss['FIBERS'].header.copy()
spec_header = rss[0].header
# Update values of sky WCS
# CRPIX1, CRPIX2
# CDELT1, CDELT2
# minx, miny
Expand All @@ -347,14 +347,16 @@ def create_cube_from_rss(rss, p=1, target_scale_arcsec=1.0, conserve_flux=True):
crpix_x = -refx / target_scale - j1min
crpix_y = -refy / target_scale - i1min
# Map the center of original field
#
#
cube[0].header['CRPIX1'] = crpix_x
cube[0].header['CRPIX2'] = crpix_y
cube[0].header['CDELT1'] = -target_scale_arcsec / (3600.0)
cube[0].header['CDELT2'] = target_scale_arcsec / (3600.0)
sky_header['CRPIX1'] = crpix_x
sky_header['CRPIX2'] = crpix_y
sky_header['CDELT1'] = -target_scale_arcsec / (3600.0)
sky_header['CDELT2'] = target_scale_arcsec / (3600.0)

# Merge headers
# 2D from FIBERS
# WL from PRIMARY
merge_wcs(sky_header, spec_header, out=cube[0].header)

# done
return cube

Expand All @@ -371,54 +373,77 @@ def recompute_wcs(hdr):


def merge_wcs(hdr_sky, hdr_spec, out=None):
"""Merge sky WCS with spectral WCS"""
"""Merge sky WCS with spectral WCS
Works only with main WCS and B WCS
"""
if out is None:
hdr = hdr_spec.copy()
else:
hdr = out

allw = astropy.wcs.find_all_wcs(hdr_spec)
for w in allw:
ss = w.wcs.alt
wcsnames = [w.wcs.alt for w in allw]
for ss in wcsnames:
merge_wcs_alt(hdr_sky, hdr_spec, hdr, spec_suffix=ss)

return hdr


def merge_wcs_alt(hdr_sky, hdr_spec, out, spec_suffix=''):
def merge_wcs_alt(hdr_sky, hdr_spec, out, spec_suffix=' '):
"""Merge sky WCS with spectral WCS"""

hdr = out
s = spec_suffix
sf = s
if spec_suffix == ' ':
sf = ''
else:
sf = spec_suffix
# Extend header for third axis
c_crpix = 'Pixel coordinate of reference point'
c_cunit = 'Units of coordinate increment and value'
hdr.set('CUNIT1{}'.format(sf), comment=c_cunit, after='CDELT1{}'.format(sf))
hdr.set('CUNIT2{}'.format(sf), comment=c_cunit, after='CUNIT1{}'.format(sf))
hdr.set('CUNIT3{}'.format(sf), value='', comment=c_cunit, after='CUNIT2{}'.format(sf))
hdr.set('CRPIX2{}'.format(sf), value=1, comment=c_crpix, after='CRPIX1{}'.format(sf))
hdr.set('CRPIX3{}'.format(sf), value=1, comment=c_crpix, after='CRPIX2{}'.format(sf))
hdr.set('CDELT3{}'.format(sf), after='CDELT2{}'.format(sf))
hdr.set('CTYPE3{}'.format(sf), after='CTYPE2{}'.format(sf))
hdr.set('CRVAL3{}'.format(sf), after='CRVAL2{}'.format(sf))
wcsname_s = 'WCSNAME{}'.format(sf)
if wcsname_s in hdr:
prev = wcsname_s
else:
prev = 'CTYPE1{}'.format(sf)

hdr.set('WCSAXES{}'.format(sf), value=3, before=prev)
if sf != '':
hdr.set('WCSNAME{}'.format(sf), value='', after='PC3_3')
hdr.set('CTYPE1{}'.format(sf), value='', after='WCSNAME{}'.format(sf))
hdr.set('CRPIX1{}'.format(sf), value=1.0, after='CTYPE1{}'.format(sf))
hdr.set('CRVAL1{}'.format(sf), value=0.0, after='CRPIX1{}'.format(sf))
hdr.set('CDELT1{}'.format(sf), value=1.0, after='CRVAL1{}'.format(sf))
hdr.set('CUNIT1{}'.format(sf), value='deg', comment=c_cunit, after='CDELT1{}'.format(sf))
hdr.set('CTYPE2{}'.format(sf), after='CUNIT1{}'.format(sf))
if sf != '':
hdr.set('CRPIX2{}'.format(sf), value=1.0, after='CTYPE2{}'.format(sf))
hdr.set('CRVAL2{}'.format(sf), value=0.0, after='CRPIX2{}'.format(sf))
hdr.set('CDELT2{}'.format(sf), value=1.0, after='CRVAL2{}'.format(sf))
hdr.set('CUNIT2{}'.format(sf), value='deg', comment=c_cunit, after='CDELT2{}'.format(sf))
hdr.set('CRPIX2{}'.format(sf), value=1, comment=c_crpix, after='CTYPE2{}'.format(sf))
hdr.set('CTYPE3{}'.format(sf), after='CUNIT2{}'.format(sf))
hdr.set('CRPIX3{}'.format(sf), value=1, comment=c_crpix, after='CTYPE3{}'.format(sf))
hdr.set('CRVAL3{}'.format(sf), after='CRPIX3{}'.format(sf))
hdr.set('CDELT3{}'.format(sf), after='CRVAL3{}'.format(sf))
hdr.set('CUNIT3{}'.format(sf), comment=c_cunit, after='CDELT3{}'.format(sf))
c_pc = 'Coordinate transformation matrix element'
hdr.set('PC1_1{}'.format(sf), value=1.0, comment=c_pc, after='CRVAL3{}'.format(sf))
hdr.set('PC1_1{}'.format(sf), value=1.0, comment=c_pc, after='CUNIT3{}'.format(sf))
hdr.set('PC1_2{}'.format(sf), value=0.0, comment=c_pc, after='PC1_1{}'.format(sf))
hdr.set('PC2_1{}'.format(sf), value=0.0, comment=c_pc, after='PC1_2{}'.format(sf))
hdr.set('PC2_2{}'.format(sf), value=1.0, comment=c_pc, after='PC2_1{}'.format(sf))
hdr.set('PC3_3{}'.format(sf), value=1.0, comment=c_pc, after='PC2_2{}'.format(sf))

# Mapping, which keyword comes from each header
mappings = [('CRPIX3', 'CRPIX1', s, 0),
('CDELT3', 'CDELT1', s, 0),
('CRVAL3', 'CRVAL1', s, 0),
('CTYPE3', 'CTYPE1', s, 0),
mappings = [('CRPIX3', 'CRPIX1', sf, 0),
('CDELT3', 'CDELT1', sf, 0),
('CRVAL3', 'CRVAL1', sf, 0),
('CTYPE3', 'CTYPE1', sf, 0),
('CRPIX1', 'CRPIX1', '', 1),
('CDELT1', 'CDELT1', '', 1),
('CRVAL1', 'CRVAL1', '', 1),
('CTYPE1', 'CTYPE1', '', 1),
('CUNIT1', 'CUNIT1', '', 1),
('CUNIT3', 'CUNIT1', sf, 0),
('PC1_1', 'PC1_1', '', 1),
('PC1_2', 'PC1_2', '', 1),
('CRPIX2', 'CRPIX2', '', 1),
Expand All @@ -429,10 +454,13 @@ def merge_wcs_alt(hdr_sky, hdr_spec, out, spec_suffix=''):
('PC2_1', 'PC2_1', '', 1),
('PC2_2', 'PC2_2', '', 1),
('LONPOLE', 'LONPOLE', '', 1),
('RADESYS', 'READESYS', '', 1),
('specsys', 'SPECSYS', s, 0),
('ssysobs', 'SSYSOBS', s, 0),
('velosys', 'VELOSYS', s, 0)
('LATPOLE', 'LATPOLE', '', 1),
('RADESYS', 'RADESYS', '', 1),
('EQUINOX', 'EQUINOX', '', 1),
('WCSNAME', 'WCSNAME', sf, 0),
('specsys', 'SPECSYS', sf, 0),
('ssysobs', 'SSYSOBS', sf, 0),
('velosys', 'VELOSYS', sf, 0)
]

hdr_in = {}
Expand Down
59 changes: 29 additions & 30 deletions megaradrp/processing/tests/test_cube.py
View file
Original file line number Diff line number Diff line change
@@ -1,45 +1,22 @@

import pytest

import numpy as np
import astropy.wcs

from megaradrp.tests.simpleobj import create_spec_header, create_sky_header
from megaradrp.tests.simpleobj import create_spec_header2, create_sky_header2
from megaradrp.processing.wavecalibration import header_add_barycentric_correction

from ..cube import create_cube, merge_wcs, merge_wcs_alt
from ..cube import create_cube, merge_wcs


def test_create_cube_raise():
with pytest.raises(ValueError):
create_cube(None, None, 3)


def test_merge_wcs():
hdr1 = create_spec_header()
hdr1 = header_add_barycentric_correction(hdr1)
hdr2 = create_sky_header()
res = merge_wcs(hdr2, hdr1)
cunit3 = res['CUNIT3']
assert cunit3 == ''


def test_merge_wcs_2():
import astropy.wcs
hdr_sky = create_sky_header()
hdr_spec = create_spec_header()
hdr_spec = header_add_barycentric_correction(hdr_spec)
allw = astropy.wcs.find_all_wcs(hdr_spec)
out = hdr_spec.copy()
for w in allw:
ss = w.wcs.alt
merge_wcs_alt(hdr_sky, hdr_spec, out, spec_suffix=ss)

assert True


def test_merge2_wcs():
hdr_sky = create_sky_header()
hdr_spec = create_spec_header()
def test_sub_wcs():
hdr_sky = create_sky_header2()
hdr_spec = create_spec_header2()
hdr_spec = header_add_barycentric_correction(hdr_spec)
wcs_sky = astropy.wcs.WCS(header=hdr_sky)
wcs_spec = astropy.wcs.WCS(header=hdr_spec, key='B')
Expand All @@ -53,4 +30,26 @@ def test_merge2_wcs():
wcs3.wcs.ssysobs = wcs_spec.wcs.ssysobs
wcs3.wcs.velosys = wcs_spec.wcs.velosys
hdr3 = wcs3.to_header(key='B')
assert True
assert True


def test_merge_wcs():
hdr_spec = create_spec_header2()
hdr_spec = header_add_barycentric_correction(hdr_spec)
hdr_sky = create_sky_header2()
out = hdr_spec.copy()
merge_wcs(hdr_sky, hdr_spec, out=out)

wcs0 = astropy.wcs.WCS(header=out)
wcsB = astropy.wcs.WCS(header=out, key='B')

assert list(wcs0.wcs.ctype) == ['RA---TAN', 'DEC--TAN', 'AWAV']
assert np.allclose([out['CRVAL1'], out['CRVAL2'], out['CRVAL3']],
[hdr_sky['CRVAL1'], hdr_sky['CRVAL2'], hdr_spec['CRVAL1']])
assert np.allclose([out['CDELT1'], out['CDELT2'], out['CDELT3']],
[hdr_sky['CDELT1'], hdr_sky['CDELT2'], hdr_spec['CDELT1']])
assert list(wcsB.wcs.ctype) == ['RA---TAN', 'DEC--TAN', 'AWAV']
assert np.allclose([out['CRVAL1B'], out['CRVAL2B'], out['CRVAL3B']],
[hdr_sky['CRVAL1'], hdr_sky['CRVAL2'], hdr_spec['CRVAL1B']])
assert np.allclose([out['CDELT1B'], out['CDELT2B'], out['CDELT3B']],
[hdr_sky['CDELT1'], hdr_sky['CDELT2'], hdr_spec['CDELT1B']])
55 changes: 55 additions & 0 deletions megaradrp/tests/simpleobj.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,34 @@ def create_sky_header():
return wcsl.to_header()


def create_sky_header2():
hdr = fits.Header()
hdr['DATE-OBS'] = '2017-08-23T21:38:30.55'
# The values of CRPIX and CDELT
# are not in pixels
# CRPIX is in mm off the center of the focal plane
# and CDELT in deg / mm
hdr['WCSAXES'] = 2
hdr['CTYPE1'] = 'RA---TAN'
hdr['CRPIX1'] = 0.0
hdr['CRVAL1'] = 255.876802773371
hdr['CDELT1'] = -0.000336666666666667
hdr['CUNIT1'] = 'deg'
hdr['CRPIX2'] = 0
hdr['CRVAL2'] = 45.6801440902947
hdr['CDELT2'] = 0.000336666666666667
hdr['CTYPE2'] = 'DEC--TAN'
hdr['CUNIT2'] = 'deg'
hdr['PC1_1'] = 0.0104368794387827
hdr['PC1_2'] = 0.999945534290533
hdr['PC2_1'] = -0.999945534290533
hdr['PC2_2'] = 0.0104368794387827
hdr['LONPOLE'] = 180.0
hdr['LATPOLE'] = hdr['CRVAL2']
hdr['RADESYS'] = 'FK5'
hdr['EQUINOX'] = 2000.
return hdr

def create_spec_header():
hdr = fits.Header()
hdr['DATE-OBS'] = '2017-08-23T21:38:30.55'
Expand All @@ -53,6 +81,33 @@ def create_spec_header():
hdr['CDELT1'] = 1.86
hdr['CUNIT1'] = 'nm'

hdr['CRPIX2'] = 0
hdr['CRVAL2'] = 0
hdr['CDELT2'] = 1
hdr['CTYPE2'] = ''
return hdr


def create_spec_header2():
hdr = fits.Header()
hdr['DATE-OBS'] = '2020-06-24T02:53:22.03'
# GTC
hdr['OBSGEO-X'] = 5327285.0921
hdr['OBSGEO-Y'] = -1718777.1125
hdr['OBSGEO-Z'] = 3051786.7327
hdr['OBSGEO-L'] = -17.881600
hdr['OBSGEO-B'] = 28.760600
hdr['OBSGEO-H'] = 2322.994

hdr['RADEG'] = 255.876802773371
hdr['DECDEG'] = 45.6801440902947

hdr['CTYPE1'] = 'AWAV'
hdr['CRPIX1'] = 1
hdr['CRVAL1'] = 6030
hdr['CDELT1'] = 0.31
hdr['CUNIT1'] = 'Angstrom'

hdr['CRPIX2'] = 0
hdr['CRVAL2'] = 0
hdr['CDELT2'] = 1
Expand Down

0 comments on commit 4903c77

Please sign in to comment.