Source Detection and Fit

A lot of efforts have been put to create robust and fast methods for finding and fitting sources in the datacubes.

Starting from a 2D map where sources should shine over the background (namely a detection frame), we show here how to use the methods to find the positions and fit both spectrally and spatially the source candidates.

Details on the methods implementation can be found in sitelle.fit and sitelle.source.

Note

A lot of these methods are parallelizable in order to be used in scripts. See here how to do that.

In [1]:
%load_ext autoreload
%autoreload 2

%matplotlib notebook
import numpy as np
import pandas as pd

from sitelle.constants import FITS_DIR
from sitelle.process import SpectralCubePatch as SpectralCube

from orb.utils import io

from sitelle.fit import fit_SN2
from sitelle.plot import plot_map, plot_scatter, plot_spectra
from sitelle.source import extract_point_source, get_sources

SN2_ORCS = SpectralCube(FITS_DIR/'orig/M31_SN2.merged.cm1.1.0.hdf5')
SN2_ORCS.set_wcs(FITS_DIR/'M31_SN2.1.0.ORCS/M31_SN2.1.0.wcs.deep_frame.fits')
SN2_ORCS.set_dxdymaps(FITS_DIR/'M31_SN2.1.0.ORCS/M31_SN2.1.0.wcs.dxmap.fits',
                      FITS_DIR/'M31_SN2.1.0.ORCS/M31_SN2.1.0.wcs.dymap.fits')
SN2_ORCS.correct_wavelength(FITS_DIR/'M31_SN2.1.0.ORCS/M31_SN2.1.0.skymap.fits')
/usr/local/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
INFO| Data shape : (2048, 2064, 556)
WARNING| /Users/Barth/Documents/M31/orcs/orcs/core.py:1596: UserWarning: Malformed spectral cube. The number of steps in the header (900) does not correspond to the real size of the data cube (556)
  warnings.warn('Malformed spectral cube. The number of steps in the header ({}) does not correspond to the real size of the data cube ({})'.format(step_nb, self.dimz))

INFO| Cube is in WAVENUMBER (cm-1)
INFO| Cube is CALIBRATED in wavenumber
WARNING| /Users/Barth/Documents/M31/orb/orb/core.py:416: UserWarning: Parameter already defined
  warnings.warn('Parameter already defined')

INFO| Cube is in WAVENUMBER (cm-1)
INFO| Cube is CALIBRATED in wavenumber
WARNING| /Users/Barth/Documents/M31/orcs/orcs/core.py:1450: UserWarning: dxmap reshaped from (200, 200) to (2048, 2064)
  format(dxmap.shape, self.dimx, self.dimy))

WARNING| /Users/Barth/Documents/M31/orcs/orcs/core.py:1457: UserWarning: dymap reshaped from (200, 200) to (2048, 2064)
  format(dymap.shape, self.dimx, self.dimy))

(see here for details on the notebook header)

Detection

We assume that a detection frame has already been built. We go over it to detect positions of the source candidates, using sitelle.source.get_sources. See sitelle.source for more details.

In [2]:
# We load and look at the detection frame
SN2_detection_frame = io.read_fits(FITS_DIR/'SN2/detection_frame.fits')
In [3]:
snr_threshold = 8 # Fixed signal-to-noise ratio threshold
SN2_sources = get_sources(SN2_detection_frame,
                          threshold=snr_threshold,
                          mode='SEGM')

SN2_sources
INFO| Detecting
WARNING| /usr/local/lib/python2.7/site-packages/photutils/segmentation/detect.py:123: RuntimeWarning: invalid value encountered in greater
  check_normalization=True) > threshold)

INFO| Deblending
INFO| Retieving properties
INFO| Filtering Quantity columns
Out[3]:
id xcentroid ycentroid sky_centroid sky_centroid_icrs source_sum source_sum_err background_sum background_mean background_at_centroid ... eccentricity orientation ellipticity elongation covar_sigx2 covar_sigxy covar_sigy2 cxx cxy cyy
0 1 1799.015929 55.714739 None None 71.180574 None None None None ... 0.792046 -1.155019 0.389539 1.638106 0.475824 -0.232195 0.899208 2.404618 1.241849 1.272426
1 2 1993.039044 63.526273 None None 1125.244804 None None None None ... 0.779919 -0.816430 0.374119 1.597748 3.966933 -1.778681 4.187997 0.311380 0.264492 0.294944
2 3 1523.544932 64.680744 None None 207.789919 None None None None ... 0.784809 -1.186372 0.380262 1.613586 0.891429 -0.405525 1.729777 1.255715 0.588773 0.647124
3 4 1128.344588 74.832441 None None 182.670688 None None None None ... 0.663103 1.474573 0.251472 1.335956 0.840420 0.062619 1.483131 1.193637 -0.100792 0.676377
4 5 380.656979 96.441849 None None 1505.593478 None None None None ... 0.810000 1.006158 0.413570 1.705233 3.270872 1.824305 4.995344 0.383931 -0.280424 0.251392
5 6 1725.319049 137.258052 None None 533.696530 None None None None ... 0.713386 -1.023478 0.299229 1.426999 1.845152 -0.663517 2.529558 0.598406 0.313930 0.436499
6 7 1728.159251 183.295429 None None 68.944996 None None None None ... 0.719897 -0.519687 0.305919 1.440754 0.651145 -0.166774 0.455075 1.694840 1.242238 2.425064
7 8 973.679613 199.060425 None None 1606.124847 None None None None ... 0.712291 1.404945 0.298116 1.424736 2.245869 0.366329 4.373040 0.451430 -0.075633 0.231842
8 9 881.669442 203.291408 None None 198.724241 None None None None ... 0.753854 1.247243 0.342958 1.521973 0.810944 0.284004 1.562628 1.316956 -0.478707 0.683450
9 10 1505.959018 205.014543 None None 468.669756 None None None None ... 0.744913 -1.205079 0.332839 1.498888 1.308970 -0.470044 2.356416 0.822904 0.328296 0.457117
10 11 1746.763690 208.003640 None None 773.474223 None None None None ... 0.691606 -1.008548 0.277725 1.384514 2.320459 -0.761242 3.048965 0.469396 0.234390 0.357240
11 12 597.085607 214.813913 None None 265.430982 None None None None ... 0.559345 0.865776 0.171065 1.206367 1.204609 0.227251 1.278309 0.858952 -0.305400 0.809430
12 13 463.268744 218.559349 None None 77.706591 None None None None ... 0.700721 1.229102 0.286565 1.401668 0.450618 0.123828 0.754757 2.323947 -0.762546 1.387482
13 14 926.868065 254.273642 None None 148.410118 None None None None ... 0.665099 -1.457984 0.253245 1.339127 0.569972 -0.050072 1.006264 1.762175 0.175371 0.998138
14 15 956.568499 255.670040 None None 74.877226 None None None None ... 0.867288 -1.456329 0.502194 2.008813 0.245308 -0.081271 0.942850 4.196345 0.723425 1.091793
15 16 346.295450 259.514618 None None 146.303972 None None None None ... 0.648809 0.649715 0.239049 1.314145 0.907008 0.217412 0.786027 1.180816 -0.653219 1.362560
16 17 1624.377497 264.910053 None None 1088.191316 None None None None ... 0.733938 -1.068608 0.320783 1.472284 2.380641 -0.923070 3.554657 0.467085 0.242584 0.312818
17 18 1485.591063 292.288122 None None 863.582306 None None None None ... 0.672237 -1.230141 0.259664 1.350738 1.826939 -0.434366 2.898346 0.567588 0.170125 0.357772
18 19 858.487701 292.747787 None None 329.848964 None None None None ... 0.697446 1.360308 0.283362 1.395405 0.958410 0.178115 1.754021 1.063464 -0.215982 0.581085
19 20 1805.641201 300.530231 None None 109.346757 None None None None ... 0.855054 -0.728189 0.481461 1.928496 1.038436 -0.557735 0.910246 1.435349 1.758962 1.637488
20 21 1754.809356 302.778597 None None 535.540122 None None None None ... 0.581543 -0.920014 0.186484 1.229233 1.751503 -0.363270 1.951977 0.593860 0.221039 0.532869
21 22 1457.107482 314.458123 None None 357.555907 None None None None ... 0.669662 -1.097012 0.257334 1.346500 1.163532 -0.328487 1.635760 0.911106 0.365930 0.648079
22 23 1605.143219 314.640817 None None 461.419309 None None None None ... 0.661483 -1.039160 0.250039 1.333404 1.399873 -0.396597 1.840977 0.760783 0.327787 0.578497
23 24 43.725146 316.604942 None None 1656.579232 None None None None ... 0.749669 0.700779 0.338187 1.511001 4.249578 1.536000 3.724659 0.276537 -0.228081 0.315510
24 25 1411.918205 326.043996 None None 181.876904 None None None None ... 0.759465 -1.554233 0.349452 1.537166 0.565909 -0.012768 1.336463 1.767450 0.033770 0.748405
25 26 1644.461975 332.235644 None None 190.637693 None None None None ... 0.611188 -1.027621 0.208514 1.263447 0.973379 -0.221529 1.206550 1.072150 0.393706 0.864953
26 27 1470.304241 337.345915 None None 193.180665 None None None None ... 0.560127 -1.362723 0.171593 1.207137 0.840690 -0.076197 1.185506 1.196469 0.153804 0.848465
27 28 1301.394119 345.526065 None None 434.407504 None None None None ... 0.722568 -1.269971 0.308700 1.446551 1.076584 -0.303727 1.961359 0.971298 0.300821 0.533143
28 29 1770.531013 365.323643 None None 104.171702 None None None None ... 0.568166 -1.133353 0.177086 1.215194 0.597070 -0.100612 0.765157 1.712798 0.450439 1.336537
29 30 759.853129 375.135418 None None 154.626239 None None None None ... 0.258135 0.640725 0.033891 1.035080 0.763071 0.024961 0.748210 1.311925 -0.087533 1.337984
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
262 263 1240.798481 1771.740614 None None 110.350568 None None None None ... 0.464460 0.913862 0.114406 1.129185 0.801533 0.096697 0.852344 1.264922 -0.287006 1.189515
263 264 764.911556 1772.462578 None None 219.895532 None None None None ... 0.551196 -1.376755 0.165624 1.198501 0.705766 -0.057345 0.986310 1.423625 0.165543 1.018692
264 265 1794.985256 1773.832564 None None 816.964128 None None None None ... 0.669823 0.750642 0.257479 1.346763 2.354356 0.665883 2.261632 0.463327 -0.272831 0.482323
265 266 304.018033 1784.224528 None None 708.398924 None None None None ... 0.656510 -0.686795 0.245683 1.325702 2.001495 -0.511620 1.797048 0.538841 0.306816 0.600143
266 267 39.544627 1784.183048 None None 94.053607 None None None None ... 0.842093 -0.441039 0.460668 1.854147 1.057070 -0.332306 0.509973 1.189719 1.550476 2.466046
267 268 1075.005076 1790.272941 None None 40.952773 None None None None ... 0.735696 -0.005918 0.322688 1.476424 0.432548 -0.001385 0.198444 2.311935 0.032281 5.039316
268 269 1252.449152 1792.283980 None None 60.673446 None None None None ... 0.743792 1.268415 0.331589 1.496086 0.247415 0.078477 0.474506 4.265567 -1.410940 2.224132
269 270 802.334658 1796.206982 None None 68.155110 None None None None ... 0.427861 -0.415040 0.096155 1.106385 0.489190 -0.034058 0.426904 2.055615 0.327993 2.355533
270 271 1489.923974 1820.615449 None None 474.607814 None None None None ... 0.434906 1.129648 0.099524 1.110524 1.346863 0.116355 1.538340 0.747349 -0.113054 0.654327
271 272 573.296475 1821.237799 None None 395.789768 None None None None ... 0.590216 -0.869429 0.192755 1.238781 1.206272 -0.260008 1.294499 0.866516 0.348090 0.807458
272 273 383.131914 1826.290120 None None 107.846309 None None None None ... 0.726492 -1.243639 0.312825 1.455234 0.451794 -0.137766 0.811011 2.334312 0.793057 1.300388
273 274 929.104841 1827.577246 None None 643.263681 None None None None ... 0.445952 -1.347671 0.104943 1.117247 1.303512 -0.068989 1.591905 0.768922 0.066646 0.629622
274 275 1089.532051 1830.495321 None None 40.864115 None None None None ... 0.196264 0.836791 0.019449 1.019835 0.248973 0.004873 0.249978 4.018037 -0.156664 4.001877
275 276 450.812324 1840.295728 None None 109.663206 None None None None ... 0.767995 -0.507498 0.359544 1.561387 0.821824 -0.239205 0.524677 1.402981 1.279262 2.197550
276 277 1029.409056 1842.112202 None None 1210.795727 None None None None ... 0.218697 -0.561640 0.024207 1.024808 2.079917 -0.045458 2.036278 0.481023 0.021477 0.491332
277 278 1575.544536 1855.967439 None None 60.739555 None None None None ... 0.770802 -1.562028 0.362925 1.569673 0.248017 -0.003183 0.610985 4.032259 0.042010 1.636811
278 279 463.372137 1862.210632 None None 53.053791 None None None None ... 0.773794 1.245381 0.366563 1.578689 0.233651 0.091642 0.474334 4.630795 -1.789356 2.281070
279 280 678.868929 1867.289480 None None 1230.612734 None None None None ... 0.538046 -0.866491 0.157085 1.186359 2.188163 -0.375743 2.311123 0.470129 0.152868 0.445117
280 281 158.754594 1872.333563 None None 999.239968 None None None None ... 0.737969 -0.759790 0.325165 1.481844 2.664631 -0.977056 2.564461 0.436229 0.332405 0.453268
281 282 1009.443103 1888.507987 None None 42.009249 None None None None ... 0.287295 -0.859887 0.042158 1.044013 0.246763 -0.010572 0.249936 4.059833 0.343451 4.008285
282 283 1182.760405 1896.044937 None None 1130.890674 None None None None ... 0.357309 1.170091 0.066014 1.070679 1.925024 0.098987 2.116757 0.520726 -0.048702 0.473560
283 284 1336.260532 1895.544638 None None 182.206861 None None None None ... 0.294425 1.375911 0.044326 1.046381 0.850207 0.015277 0.924586 1.176534 -0.038880 1.081887
284 285 635.360438 1896.123060 None None 104.264509 None None None None ... 0.626063 -0.434015 0.220228 1.282426 0.732641 -0.117720 0.533216 1.415126 0.624844 1.944386
285 286 398.101122 1924.588694 None None 164.168910 None None None None ... 0.536544 -1.318896 0.156127 1.185013 0.790535 -0.075250 1.063551 1.273543 0.180215 0.946622
286 287 936.891443 1925.494030 None None 270.677859 None None None None ... 0.337355 -0.642285 0.058623 1.062273 0.968963 -0.055147 0.936503 1.035501 0.121954 1.071393
287 288 655.032965 1968.495371 None None 407.097320 None None None None ... 0.699799 -0.997779 0.285661 1.399895 1.268982 -0.432737 1.660433 0.864900 0.450815 0.660997
288 289 1020.374868 1974.576975 None None 186.977612 None None None None ... 0.482543 -0.953825 0.124128 1.141719 0.840719 -0.109312 0.917281 1.208179 0.287955 1.107337
289 290 820.784621 1978.316132 None None 1234.485713 None None None None ... 0.565465 -1.220576 0.175228 1.212456 1.834287 -0.263300 2.458929 0.553681 0.118575 0.413030
290 291 327.985887 1978.494752 None None 42.813448 None None None None ... 0.916559 -0.564968 0.600100 2.500625 0.473593 -0.236970 0.249972 4.016890 7.615899 7.610318
291 292 1335.509554 1984.502901 None None 40.828840 None None None None ... 0.198388 -0.789527 0.019876 1.020279 0.249909 -0.005017 0.249992 4.003074 0.160682 4.001747

292 rows × 35 columns

Details about all the outputs can be found in sitelle.source.get_sources method.

We can plot these detections on the detection frame for verification.

Warning

This uses astropy convention, i.e. ycentroid is actually abscisse and xcentroid columns… To inverse this:

>>> SN2_sources.rename(columns={'ycentroid':'xcentroid', 'xcentroid':'ycentroid'}, inplace=True)
In [4]:
SN2_sources.rename(columns={'ycentroid':'xcentroid', 'xcentroid':'ycentroid'}, inplace=True)

f,ax = plot_map(SN2_detection_frame, cmap='jet', colorbar=True, vmax=10.)
plot_scatter(SN2_sources['xcentroid'], SN2_sources['ycentroid'], ax=ax, c='k', marker='o', facecolors='none')

ax.set_xlim(400,1000)
ax.set_ylim(500,1100)
Out[4]:
(500, 1100)

Cosmetics

We need to transform the data a bit to be able to fit spectra.

In [5]:
### Filtering out of bounds detections
from sitelle.calibration import filter_star_list
source_list = filter_star_list(SN2_sources.as_matrix(['xcentroid', 'ycentroid']))
SN2_sources['xcentroid'] = source_list[:,0]
SN2_sources['ycentroid'] = source_list[:,1]

### Dropping null rows, i.e. bad detections
for colname in SN2_sources.columns:
    if pd.isnull(SN2_sources[colname]).all():
        SN2_sources.drop(colname, axis=1, inplace=True)
SN2_sources.dropna(inplace=True)

### Rounding positionss for spectra extraction
SN2_sources['xpos'] = SN2_sources['xcentroid'].apply(lambda x: int(round(x)))
SN2_sources['ypos'] = SN2_sources['ycentroid'].apply(lambda y: int(round(y)))

### Sorting detections on the brightness
SN2_sources.sort_values('source_sum', ascending=False, inplace=True)
SN2_sources.reset_index(drop=True, inplace=True)

## Computing physical positions
radec = SN2_ORCS.pix2world(SN2_sources.as_matrix(['xcentroid', 'ycentroid']))
SN2_sources['ra'] = radec[:,0]
SN2_sources['dec'] = radec[:,1]

Spectral Fitting

We are going to focus on one source of these detections

In [24]:
source = SN2_sources.iloc[1]
fit_res, fit_params = fit_SN2(source, SN2_ORCS, return_fit_params=True)
In [25]:
# We display the original spectra and the fit
x,y = source[['xpos', 'ypos']].astype(int)

axis, spectra = extract_point_source(x,y, SN2_ORCS)
fitted_spectra = fit_params['fitted_vector']
f,ax = plot_spectra(axis,spectra, label='Original')
plot_spectra(axis, fitted_spectra, ax=ax, label='Fit')
Out[25]:
(<Figure size 640x480 with 2 Axes>,
 <matplotlib.axes._subplots.AxesSubplot at 0x11f07b750>)
In [26]:
## We append the fit result to the original source data
source = source.append(fit_res)

## We define some threshold to say if a line is detected or not
from sitelle.constants import SN2_LINES
snr_lim = 5
for line_name in SN2_LINES:
    line_name = line_name.lower().replace('[', '').replace(']', '')
    source['%s_detected'%line_name] = source['snr_%s'%line_name] > snr_lim
    source['flux_%s'%line_name] = np.where(source['%s_detected'%line_name],
                                            source['flux_%s'%line_name],
                                            snr_lim*source['flux_%s_err'%line_name])

The sitelle.fit.fit_SN2 method gathers a lot of different steps and uses a lot of default values to facilitate the fit. It is largely tunable through the keywords, to modify the LSF to use for the fit, guess values etc…

In [9]:
## An example of tweaking the keywords
## We only fit the 0III5007 line with a gaussian
kwargs_spec = {'fmodel':'gaussian'}
fit_res, fit_params = fit_SN2(source, SN2_ORCS, return_fit_params=True,
                             lines=['[OIII]5007'], kwargs_spec=kwargs_spec)

fitted_spectra = fit_params['fitted_vector']
f,ax = plot_spectra(axis,spectra, label='Original')
plot_spectra(axis, fitted_spectra, ax=ax, label='Fit')
Out[9]:
(<Figure size 640x480 with 2 Axes>,
 <matplotlib.axes._subplots.AxesSubplot at 0x112b289d0>)

Spatial Fitting

We can also look spatially at this source, to evaluate its PSF etc.

In [33]:
analysis = analyse_source(source, SN2_ORCS, plot=True)
print analysis
erf_amplitude                                                   1.71839e-14
erf_amplitude_err                                               7.86383e-17
erf_fwhm                                                            4.44762
erf_ks_pvalue                                                      0.717872
erf_xfwhm                                                           2.98996
erf_yfwhm                                                           1.99384
flux_err_r                [7.247704606718017e-18, 1.5121543475326628e-17...
flux_fraction_3                                                    0.273046
flux_map_ks_pvalue                                              8.70415e-14
flux_r                    [0.0, 6.071588579643632e-16, 2.370023932326401...
model_flux_fraction_15                                             0.996758
modeled_flux_r            [-9.719107571595448e-19, 6.153336475257311e-16...
psf_amplitude                                                     1.605e-14
psf_ks_pvalue                                                   2.77239e-09
psf_snr                                                             82.0693
psf_xfwhm                                                           5.41384
psf_yfwhm                                                           3.58781
dtype: object
In [34]:
source = source.append(analysis)

Parallel Computations

After identification of the positions (i.e. one we have the SN2_sources DataFrame), one can perform the analysis steps on all the sources at once using fast parallelization.

Here goes:

Spectral Fit

In [ ]:
from sitelle.parallel import parallel_apply_over_df

v_min = -1200.
v_max = 300.
SN2_ORCS.inputparams = {}
# kwargs_spec = {'fmodel':'sincgauss', 'sigma_def':['1'], 'sigma_cov':20}
kwargs_spec = {}
fit_dataframe = parallel_apply_over_df(SN2_sources, fit_SN2, axis=1, args=(SN2_ORCS, ), v_min = v_min, v_max = v_max, kwargs_spec = kwargs_spec,modules = 'from sitelle.region import centered_square_region','from sitelle.fit import guess_source_velocity, refine_velocity_guess', 'from sitelle.utils import stats_without_lines','from orcs.utils import fit_lines_in_spectrum')

The modules argument is quite important : because the different threads in the parallelization do not share the common environnement we have built here by importing all the modules we needed, we have to explicitely say what needs to be imported so that the sitelle.fit.fit_SN2 method works correctly.

fit_dataframe is now a DataFrame of the same shape as SN2_sources, where each row corresponds to the fit parameters of each source

In [ ]:
## We merge it with the original parameters, and evaluate detections
SN2_fit_res = pd.merge(SN2_sources, fit_dataframe, left_index=True, right_index=True)
snr_lim = 3.

for line_name in SN2_LINES:
    line_name = line_name.lower().replace('[', '').replace(']', '')
    SN2_fit_res['%s_detected'%line_name] = SN2_fit_res['snr_%s'%line_name] > snr_lim
    SN2_fit_res['flux_%s'%line_name] = np.where(SN2_fit_res['%s_detected'%line_name],
                                            SN2_fit_res['flux_%s'%line_name],
                                            snr_lim*SN2_fit_res['flux_%s_err'%line_name])

Spatial Fit

In [ ]:
analysis = parallel_apply_over_df(SN2_sources, analyse_source, axis=1, args=(SN2_ORCS,))