Wavelength recalibration with the sky lines (Mendel OH bands)¶
Introduction¶
Data can be recalibrated in wavelength easily when one OH sky lines are visible in most parts of the field of view. Cubes observed in the red part of the spectrum (e.g. SN3 filter) are especially interesting in this respect.
The general idea is to extract integrated spectra at different positions in the field of view and measure the velocity of the sky lines (which should be 0 if the calibration was perfect) (see Martin et al. 2017a).
Then the correction map can be infered at each pixel by fitting a model.
The calibration model is based on a simple modeling of the interferometer.
First step: checking the calibration¶
[1]:
# import base class for the manipulation of a SITELLE spectral cube: HDFCube
from orcs.process import SpectralCube
import pylab as pl
[2]:
# load spectral cube
cube = SpectralCube('/home/thomas/M31_SN3.merged.cm1.1.0.hdf5')
master.fa09a|INFO| Cube is level 3
master.fa09a|INFO| shape: (2048, 2064, 840)
master.fa09a|INFO| wavenumber calibration: True
master.fa09a|INFO| flux calibration: True
master.fa09a|INFO| wcs calibration: True
[3]:
# extract and plot a spectrum of a large integrated region. Sky lines should be visible.
spectrum = cube.get_spectrum(50, 50, 20)
pl.figure(figsize=(10,6))
spectrum.plot()
pl.xlim(14500, 15500)
[3]:
(14500, 15500)
fitting the sky spectrum¶
we can also fit the integrated spectrum because there si a large number of sky lines and because we are not interested in a perfect fit (only the velocity is required), we can use a simple ‘sinc’ model with a fixed fwhm. All the lines are set to share the same velocity parameter and we set an inital guess of the velocity around 80 km/s. This general bias of 80 km/s is known (Martin et al. 2017b) and comes from the error made on the real wavelength of the calibration laser (which is falsely considered to be at 453.5 nm).
[4]:
sky_lines_cm1 = cube.get_sky_lines()
fit_res = spectrum.fit(sky_lines_cm1,
fmodel='sinc',
pos_def='1',
fwhm_def='fixed',
nofilter=False,
pos_cov=80)
print('Velocity: ', fit_res['velocity_gvar'][0])
pl.figure(figsize=(10,6))
spectrum.plot()
fit_res.get_spectrum().plot(c='red')
pl.xlim(14500, 15500)
Velocity: 77.5(1.5)
[4]:
(14500, 15500)
Second step: Mapping the sky velocity¶
This process, like other parallel processes may work better if run from a python script with the command. The following script can be used:
from orcs.process import SpectralCube
if __name__ == "__main__": # this is very important when running parallel processes
cube = SpectralCube('/home/thomas/M31_SN3.merged.cm1.1.0.hdf5')
cube.map_sky_velocity(80, div_nb=10, exclude_reg_file_path='m31_exclude.reg', no_fit=False, plot=False) # the mean velocity bias is set around 80 km/s
Warning: this process can take a long time because a spectrum is extracted and fitted for each point of a 40x40 grid by default. For a basic example it is recommended to limit the number of grid division to 10x10.
[7]:
# once the fit is done, you can see the results with the following
cube = SpectralCube('/home/thomas/M31_SN3.merged.cm1.1.0.hdf5') # it is always better to reload the data cube right before running a parallel process
cube.map_sky_velocity(80, div_nb=10, exclude_reg_file_path='m31_exclude.reg', no_fit=False) # the mean velocity bias is set around 80 km/s
master.fa09a|INFO| Cube is level 3
master.fa09a|INFO| shape: (2048, 2064, 840)
master.fa09a|INFO| wavenumber calibration: True
master.fa09a|INFO| flux calibration: True
master.fa09a|INFO| wcs calibration: True
master.fa09a|WARNING| fitting process already done ! set no_fit=True if you do not want to redo it
master.fa09a|INFO| X range: 0 2048, Y range: 0 2064
master.fa09a|INFO| Radius: 30
master.fa09a|INFO| excluding region from file m31_exclude.reg
master.fa09a|INFO| 58 regions to fit
master.fa09a|INFO| 88 sky lines to fit
loading region: Shape : circle ( Number(1861.81818),Number(1876.36364),Number(30.00000) )
master.fa09a|INFO| reloading cube
master.fa09a|INFO| Cube is level 3
master.fa09a|INFO| shape: (2048, 2064, 840)
master.fa09a|INFO| wavenumber calibration: True
master.fa09a|INFO| flux calibration: True
master.fa09a|INFO| wcs calibration: True
master.fa09a|INFO| Init of the parallel processing server with 32 threads
[==========] [100%] [completed in 4m33s]
master.fa09a|INFO| Closing parallel processing server
master.fa09a|INFO| reloading cube
master.fa09a|INFO| Cube is level 3
master.fa09a|INFO| shape: (2048, 2064, 840)
master.fa09a|INFO| wavenumber calibration: True
master.fa09a|INFO| flux calibration: True
master.fa09a|INFO| wcs calibration: True
master.fa09a|INFO| parallel processing closed
master.fa09a|INFO| Velocity of the first line (km/s): 77.3(2.0)
master.fa09a|INFO| Velocity of the first line (km/s): 75.4(2.0)
master.fa09a|INFO| Velocity of the first line (km/s): 71.3(1.9)
master.fa09a|INFO| Velocity of the first line (km/s): 70.7(2.3)
master.fa09a|INFO| Velocity of the first line (km/s): 71.7(2.8)
master.fa09a|INFO| Velocity of the first line (km/s): 66.7(2.5)
master.fa09a|INFO| Velocity of the first line (km/s): 64.1(2.3)
master.fa09a|INFO| Velocity of the first line (km/s): 58.1(2.3)
master.fa09a|INFO| Velocity of the first line (km/s): 62.6(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 61.0(2.5)
master.fa09a|INFO| Velocity of the first line (km/s): 74.0(1.8)
master.fa09a|INFO| Velocity of the first line (km/s): 75.5(1.9)
master.fa09a|INFO| Velocity of the first line (km/s): 77.2(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 72.5(3.1)
master.fa09a|INFO| Velocity of the first line (km/s): 69.7(3.3)
master.fa09a|INFO| Velocity of the first line (km/s): 55.3(2.8)
master.fa09a|INFO| Velocity of the first line (km/s): 65.0(2.9)
master.fa09a|INFO| Velocity of the first line (km/s): 58.9(3.0)
master.fa09a|INFO| Velocity of the first line (km/s): 62.7(2.4)
master.fa09a|INFO| Velocity of the first line (km/s): 60.3(2.3)
master.fa09a|INFO| Velocity of the first line (km/s): 76.4(2.0)
master.fa09a|INFO| Velocity of the first line (km/s): 77.9(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 73.4(2.6)
master.fa09a|INFO| Velocity of the first line (km/s): 48.3(2.6)
master.fa09a|INFO| Velocity of the first line (km/s): 61.7(2.8)
master.fa09a|INFO| Velocity of the first line (km/s): 78.6(1.9)
master.fa09a|INFO| Velocity of the first line (km/s): 75.4(2.5)
master.fa09a|INFO| Velocity of the first line (km/s): 104.3(3.2)
master.fa09a|INFO| Velocity of the first line (km/s): 65.0(2.7)
master.fa09a|INFO| Velocity of the first line (km/s): 75.6(2.4)
master.fa09a|INFO| Velocity of the first line (km/s): 78.8(2.9)
master.fa09a|INFO| Velocity of the first line (km/s): 65.2(2.4)
master.fa09a|INFO| Velocity of the first line (km/s): 80.8(2.6)
master.fa09a|INFO| Velocity of the first line (km/s): 74.4(3.0)
master.fa09a|INFO| Velocity of the first line (km/s): 69.8(2.3)
master.fa09a|INFO| Velocity of the first line (km/s): 78.5(2.4)
master.fa09a|INFO| Velocity of the first line (km/s): 65.4(2.5)
master.fa09a|INFO| Velocity of the first line (km/s): 65.9(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 74.6(2.5)
master.fa09a|INFO| Velocity of the first line (km/s): 58.8(2.7)
master.fa09a|INFO| Velocity of the first line (km/s): 67.0(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 70.6(2.0)
master.fa09a|INFO| Velocity of the first line (km/s): 80.9(2.7)
master.fa09a|INFO| Velocity of the first line (km/s): 80.1(3.5)
master.fa09a|INFO| Velocity of the first line (km/s): 70.8(3.1)
master.fa09a|INFO| Velocity of the first line (km/s): 78.1(2.6)
master.fa09a|INFO| Velocity of the first line (km/s): 70.2(1.8)
master.fa09a|INFO| Velocity of the first line (km/s): 74.5(1.4)
master.fa09a|INFO| Velocity of the first line (km/s): 79.4(2.9)
master.fa09a|INFO| Velocity of the first line (km/s): 79.3(2.7)
master.fa09a|INFO| Velocity of the first line (km/s): 77.5(3.3)
master.fa09a|INFO| Velocity of the first line (km/s): 77.7(3.1)
master.fa09a|INFO| Velocity of the first line (km/s): 76.1(3.5)
master.fa09a|INFO| Velocity of the first line (km/s): 74.7(2.9)
master.fa09a|INFO| Velocity of the first line (km/s): 75.2(2.2)
master.fa09a|INFO| Velocity of the first line (km/s): 75.7(1.7)
master.fa09a|INFO| Velocity of the first line (km/s): 71.6(1.5)
master.fa09a|INFO| Velocity of the first line (km/s): 73.8(1.3)
master.fa09a|INFO| x: [ 186.18181818 186.18181818 186.18181818 186.18181818 186.18181818
186.18181818 186.18181818 186.18181818 186.18181818 186.18181818
372.36363636 372.36363636 372.36363636 372.36363636 372.36363636
372.36363636 372.36363636 372.36363636 372.36363636 372.36363636
558.54545455 558.54545455 558.54545455 558.54545455 558.54545455
744.72727273 744.72727273 744.72727273 744.72727273 930.90909091
930.90909091 930.90909091 1117.09090909 1117.09090909 1117.09090909
1303.27272727 1303.27272727 1303.27272727 1489.45454545 1489.45454545
1489.45454545 1489.45454545 1675.63636364 1675.63636364 1675.63636364
1675.63636364 1675.63636364 1675.63636364 1861.81818182 1861.81818182
1861.81818182 1861.81818182 1861.81818182 1861.81818182 1861.81818182
1861.81818182 1861.81818182 1861.81818182]
master.fa09a|INFO| y: [ 187.63636364 375.27272727 562.90909091 750.54545455 938.18181818
1125.81818182 1313.45454545 1501.09090909 1688.72727273 1876.36363636
187.63636364 375.27272727 562.90909091 750.54545455 938.18181818
1125.81818182 1313.45454545 1501.09090909 1688.72727273 1876.36363636
187.63636364 375.27272727 562.90909091 1688.72727273 1876.36363636
187.63636364 375.27272727 1688.72727273 1876.36363636 187.63636364
375.27272727 1876.36363636 187.63636364 1688.72727273 1876.36363636
187.63636364 1688.72727273 1876.36363636 187.63636364 1501.09090909
1688.72727273 1876.36363636 187.63636364 375.27272727 1313.45454545
1501.09090909 1688.72727273 1876.36363636 187.63636364 375.27272727
562.90909091 750.54545455 938.18181818 1125.81818182 1313.45454545
1501.09090909 1688.72727273 1876.36363636]
master.fa09a|INFO| v: [ 77.30690122 75.42723292 71.29603104 70.74703947 71.71676527
66.74086248 64.11377386 58.053346 62.58297622 61.04474574
74.01054084 75.52157459 77.22636344 72.5418529 69.69014005
55.32754889 64.99891349 58.9375885 62.70848958 60.32799584
76.36498222 77.90017749 73.41305969 48.25166759 61.66747683
78.61262775 75.37893981 104.25308835 65.04902046 75.61380729
78.84037233 65.16314474 80.77022818 74.38134783 69.81025738
78.51941433 65.42526138 65.85118873 74.60799346 58.79111662
66.95923754 70.63770242 80.9106735 80.12792716 70.77220517
78.1282894 70.221003 74.52473849 79.4390245 79.2954783
77.52699254 77.7411457 76.13299069 74.74027681 75.16744211
75.65534916 71.59209525 73.79256937]
master.fa09a|WARNING| /home/thomas/Astro/Python/ORB/Orcs/orcs/utils.py:212: RuntimeWarning: invalid value encountered in greater
vel[vel > np.nanpercentile(vel, 95)] = np.nan
master.fa09a|INFO| First laser wavelentgh calibration estimation: 543.3684614763755 nm
master.fa09a|INFO| Angle at the center of the frame: 15.494273258242744
master.fa09a|WARNING| /home/thomas/Astro/Python/ORB/Orb/orb/utils/image.py:1191: RuntimeWarning: invalid value encountered in less
calib_laser_map[np.nonzero(calib_laser_map < value_min)] = np.nan
master.fa09a|INFO| > Binning calibration map
master.fa09a|WARNING| /home/thomas/Astro/Python/ORB/Orb/orb/utils/image.py:1074: RuntimeWarning: Mean of empty slice
return np.squeeze(np.nanmean(np.nanmean(im_view, axis=3), axis=1))
master.fa09a|INFO| > Calibration laser map fit
master.fa09a|INFO| > First fit on the central portion of the calibration laser map (30.0% of the total size)
master.fa09a|INFO| > Calibration laser map fit parameters:
distance to mirror: 23.681739781375743 cm
X angle from the optical axis to the center: 0.0 degrees (Fixed)
Y angle from the optical axis to the center: 15.501669386747308 degrees
Tip-tilt angle of the detector along X: 0.0 degrees (Fixed)
Tip-tilt angle of the detector along Y: 0.0 degrees (Fixed)
Rotation angle of the detector: 0.0 degrees (Fixed)
Calibration laser wavelength: 543.3684614763755 nm (Fixed)
Error on fit: mean -1.2212972352744125e-06, std 0.049408481925003994 (in nm)
Error on fit: mean -0.0006742923017409128, std 27.27899322170293 (in km/s)
master.fa09a|INFO| > Second fit on the central portion of the calibration laser map (30.0% of the total size)
master.fa09a|INFO| > Calibration laser map fit parameters:
distance to mirror: 23.66483634860634 cm
X angle from the optical axis to the center: -0.44808541110452904 degrees
Y angle from the optical axis to the center: 15.495141935947101 degrees
Tip-tilt angle of the detector along X: 0.41142114590561524 degrees
Tip-tilt angle of the detector along Y: -1.41769820240921 degrees
Rotation angle of the detector: 0.0 degrees (Fixed)
Calibration laser wavelength: 543.3684614763755 nm (Fixed)
Error on fit: mean -3.798378046820534e-07, std 0.0007481455641519504 (in nm)
Error on fit: mean -0.0002097128366541575, std 0.41305980224865047 (in km/s)
master.fa09a|INFO| > Third fit on a larger portion of the map (50.0% of the total size)
master.fa09a|INFO| > Calibration laser map fit parameters:
distance to mirror: 23.702051086393453 cm
X angle from the optical axis to the center: -0.44853898242361845 degrees
Y angle from the optical axis to the center: 15.495504360876263 degrees
Tip-tilt angle of the detector along X: 0.25839602728725264 degrees
Tip-tilt angle of the detector along Y: -0.8882041482610137 degrees
Rotation angle of the detector: 0.0 degrees (Fixed)
Calibration laser wavelength: 543.3684614763755 nm (Fixed)
Error on fit: mean -2.4518726781407214e-06, std 0.003041168391745898 (in nm)
Error on fit: mean -0.0013537072089970703, std 1.6790641750624251 (in km/s)
master.fa09a|INFO| > Zernike polynomials fit of the residual wavefront
master.fa09a|WARNING| /home/thomas/Astro/Python/ORB/Orb/orb/ext/zern.py:300: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
wf_zern_vec = np.linalg.lstsq((zern_basismat[:, grid_vec] * weight).T, wf_w.ravel())[0]
master.fa09a|INFO| Standard deviation of the residual: 1.8507972727041484e-05
master.fa09a|INFO| Median relative error (err/val)): 0.05%
master.fa09a|INFO| > final error (std on 50.0% of the total size): 1.763e-05 nm, 9.733e-03 km/s
master.fa09a|INFO| > New calibration laser map fit parameters:
distance to mirror: 23.67734145891187 cm
X angle from the optical axis to the center: 0.1594504146076437 degrees
Y angle from the optical axis to the center: 15.503103983197793 degrees
Tip-tilt angle of the detector along X: 4.033689007793616 degrees
Tip-tilt angle of the detector along Y: -0.3016934072811556 degrees
Rotation angle of the detector: 2.2124798963045356 degrees
Calibration laser wavelength: 543.3680110441915 nm
master.fa09a|INFO| fit residual std (in km/s):
master.fa09a|INFO| median error on the data (in km/s)
master.fa09a|INFO| Data written as ./M31_SN3/M31_SN3.calibration_laser_map.fits in 0.33 s
master.fa09a|INFO| Data written as ./M31_SN3/M31_SN3.wavefront_map.fits in 0.23 s
master.fa09a|INFO| Data written as ./M31_SN3/M31_SN3.skymap.fits in 0.26 s
Ouputs¶
The most important output is the skymap which gives the modelized velocity of the sky in the field of view and permit to correct the measured velocities in the cube.This map is named *skymap.fits
. You can add it to the obtained velocity maps to correct the velocity.
In this particular case this is M31_SN3/M31_SN3.skymap.fits
The histogram above shows the difference between the measured velocity and the fitted velocity at each valid data point.
You can see that there are some data points at the center of the field which were excluded from the fit. This comes from the fact that M31 is very bright at the center and the sky lines are lost in the noise produced by the high continuum background. Modelling the interferometer makes possible to estimate the velocity of the lines at the center of the field where no velocity measurement could be done. And the result is stricking, it works particularly well even if the case of M31 (Martin et al. 2017a)
[14]:
import orb.utils.io
skymap = orb.utils.io.read_fits('./M31_SN3/M31_SN3.skymap.fits')
pl.imshow(skymap.T, origin='bottom')
pl.grid()
cb = pl.colorbar()
cb.set_label('velocity error (km/s)')
[ ]: