Source code for besl.bgps_peak_flux_extract

"""
======================
BGPS Peak Flux Extract
======================

Extract peak values from BGPS maps

"""

import numpy as _np
import pyfits as _pyfits
import pywcs as _pywcs

[docs]class Dirs(object): """ Object to hold directories for interactive path editing """ def __init__(self): self.root_dir = '/mnt/eld_data/' self.bgps1_img_dir = self.root_dir + 'BGPS/Images/v1.0.2/all/' self.bgps2_img_dir = self.root_dir + 'BGPS/Images/v2.0.0/' self.bgps1_img_filen = 'v1.0.2_{}_13pca_map50.fits' self.bgps1_med_filen = 'v1.0.2_{}_13pca_medmap50.fits' self.bgps1_rms_filen = 'v1.0.2_{}_13pca_noisemap50.fits' self.bgps2_img_filen = 'v2.0_ds2_{}_13pca_map20.fits' self.bgps2_med_filen = 'v2.0_ds2_{}_13pca_medmap20.fits' self.bgps2_rms_filen = 'v2.0_ds2_{}_13pca_noisemap20.fits'
d = Dirs()
[docs]def extract_peak_bgps_props(out_filen='bgps_pk_extract'): """ Extract peak flux and noise values from the BGPS maps in Jy/beam. Citation: Ginsburg et al. (2013). Parameters ---------- out_filen : string output catalog file in CSV format Returns ------- molcat_pk : pandas.DataFrame Output catalog in a pandas DataFrame object """ # TODO add extraction from v1 maps from besl.coord import eq2gal from besl.catalog import read_bgps, read_bgps_bounds, read_molcat, \ select_bgps_field # read in catalogs bgps_bounds = read_bgps_bounds() molcat = read_molcat() # add new columns new_cols = ['map_pk', 'medmap_pk', 'map_rms'] for col in new_cols: molcat[col] = _np.nan # match clumps for clump in molcat.cnum.values: # index parameters cindex = _np.argwhere(molcat.cnum == clump)[0][0] ra, dec = molcat.ix[cindex]['hht_ra'], molcat.ix[cindex]['hht_dec'] glon, glat = eq2gal(ra, dec) # bgps v2 field = select_bgps_field(ra, dec, coord_type='eq') bgps2_img = _pyfits.open(d.bgps2_img_dir + d.bgps2_img_filen.format(field)) bgps2_med = _pyfits.open(d.bgps2_img_dir + d.bgps2_med_filen.format(field)) bgps2_rms = _pyfits.open(d.bgps2_img_dir + d.bgps2_rms_filen.format(field)) bgps_img_wcs = _pywcs.WCS(bgps2_img[0].header) pix_x, pix_y = _np.floor(bgps_img_wcs.wcs_sky2pix(glon, glat, 0)) pix_x, pix_y = _np.floor(pix_x)[0], _np.floor(pix_y)[0] molcat['map_pk'].ix[cindex] = bgps2_img[0].data[pix_y, pix_x] molcat['medmap_pk'].ix[cindex] = bgps2_med[0].data[pix_y, pix_x] molcat['map_rms'].ix[cindex] = bgps2_rms[0].data[pix_y, pix_x] molcat_pk = molcat[['cnum', 'map_pk', 'medmap_pk', 'map_rms']] molcat_pk.to_csv(out_filen + '.csv', index=False, na_rep='-999', float_format='%.6f') print '-- Catalog written to {}.csv'.format(out_filen) return molcat_pk