Source code for grizli.catalog

"""
Catalog table tools
"""
import os
import inspect

from collections import OrderedDict
import glob
import traceback

import numpy as np

import astropy.io.fits as pyfits
import astropy.units as u
import astropy.coordinates as coord
from astropy.table import Table

from . import utils

__all__ = ["table_to_radec", 
           "table_to_regions", 
           "randomize_segmentation_labels", 
           "get_ukidss_catalog",
           "get_sdss_catalog",
           "get_twomass_catalog",
           "get_irsa_catalog",
           "get_gaia_radec_at_time",
           "get_gaia_DR2_vizier_columns",
           "get_gaia_DR2_vizier",
           "gaia_dr2_conesearch_query",
           "get_gaia_DR2_catalog",
           "gen_tap_box_query",
           "query_tap_catalog",
           "get_hubble_source_catalog",
           "get_nsc_catalog",
           "get_desdr1_catalog",
           "get_skymapper_catalog",
           "get_vexas_catalog",
           "get_vhs_catalog",
           "get_panstarrs_catalog",
           "get_radec_catalog"]

[docs]def table_to_radec(table, output='coords.radec'): """Make a "radec" ascii file with ra, dec columns from a table object """ if 'X_WORLD' in table.colnames: rc, dc = 'X_WORLD', 'Y_WORLD' elif 'x_world' in table.colnames: rc, dc = 'x_world', 'y_world' elif 'ra' in table.colnames: rc, dc = 'ra', 'dec' else: raise ValueError("Couldn't identify sky coordinates (x_world, ra)") table[rc].format = '.8f' table[dc].format = '.8f' table[rc, dc].write(output, format='ascii.commented_header', overwrite=True)
[docs]def table_to_regions(table, output='ds9.reg', comment=None, header='global color=green width=1', size=0.5, use_ellipse=False, use_world=True, verbose=True, unit_offset_colums=['x'], scale_major=1, extra=None): """Make a DS9 region file from a table object """ fp = open(output, 'w') fp.write(header+'\n') if ('X_WORLD' in table.colnames) & (use_world): xc, yc = 'X_WORLD', 'Y_WORLD' is_world = True elif ('ra' in table.colnames) & (use_world): xc, yc = 'ra', 'dec' is_world = True elif ('x_world' in table.colnames) & (use_world): xc, yc = 'x_world', 'y_world' is_world = True elif 'x_image' in table.colnames: xc, yc = 'x_image','y_image' is_world = False elif 'x' in table.colnames: xc, yc = 'x','y' is_world = False else: raise ValueError("Couldn't identify ra/dec or x/y columns") xdata = table[xc]*1 ydata = table[yc]*1 if xc in unit_offset_colums: if verbose: print(f'Unit offset for columns {xc}, {yc}') xdata += 1 ydata += 1 if is_world: fp.write('icrs\n') sec='"' else: fp.write('image\n') sec='' # GAIA if 'solution_id' in table.colnames: e = np.sqrt(table['ra_error']**2+table['dec_error']**2)/1000. e = np.maximum(e, 0.1) else: e = np.ones(len(table))*size if use_ellipse: if is_world: if 'a_world' in table.colnames: amaj = table['a_world']*3600*scale_major amin = table['b_world']*3600*scale_major etheta = table['theta_world']/np.pi*180#+90 else: use_ellipse = False else: if 'a_image' in table.colnames: amaj = table['a_image']*scale_major amin = table['b_image']*scale_major etheta = table['theta_image']/np.pi*180#+90 else: use_ellipse = False if verbose: print(f'{output}: x = {xc}, y={yc}, ellipse={use_ellipse}') if use_ellipse: regstr = 'ellipse({x:.7f}, {y:.7f}, {a:.3f}{sec}, {b:.3f}{sec}, {theta:.1f})' lines = [regstr.format(x=xdata[i], y=ydata[i], a=amaj[i], b=amin[i], sec=sec, theta=etheta[i]) for i in range(len(table))] else: regstr = 'circle({0:.7f}, {1:.7f}, {2:.3f}{3})' lines = [regstr.format(xdata[i], ydata[i], e[i], sec) for i in range(len(table))] if comment is not None: for i in range(len(table)): lines[i] += ' # text={{{0}}}'.format(comment[i]) if extra is not None: for i, li in enumerate(lines): if '#' not in li: lines[i] += ' #' lines[i] += ' '+extra[i] # newline lines = [l+'\n' for l in lines] fp.writelines(lines) fp.close()
[docs]def randomize_segmentation_labels(seg, random_seed=1, fill_value=np.nan): """ Randomize labels on a segmentation image for easier visualization Parameters ---------- seg : str or array Segmentation filename or integer array random_seed : int Random number seed for `numpy.random.seed` fill_value : scalar Value to insert where ``seg == 0`` Returns ------- rand_ids : array 1D array of the randomized labels with size `max(seg)+1` rand_seg : array 2D rray with randomized labels """ np.random.seed(random_seed) if isinstance(seg, str): seg_im = pyfits.open(seg) seg_data = seg_im[0].data else: seg_data = seg imax = seg_data.max() rand_ids = np.append(fill_value, np.argsort(np.random.rand(imax))+1) return rand_ids, rand_ids[seg_data]
[docs]def get_ukidss_catalog(ra=165., dec=34.8, radius=3, database='UKIDSSDR9PLUS', programme_id='LAS'): """Query for objects in the UKIDSS catalogs Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin Returns ------- table : `~astropy.table.Table` Result of the query """ from astroquery.ukidss import Ukidss coo = coord.SkyCoord(ra*u.deg, dec*u.deg) table = Ukidss.query_region(coo, radius=radius*u.arcmin, database=database, programme_id=programme_id) return table
[docs]def get_sdss_catalog(ra=165.86, dec=34.829694, radius=3): """Query for objects in the SDSS photometric catalog Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin Returns ------- table : `~astropy.table.Table` Result of the query """ from astroquery.sdss import SDSS coo = coord.SkyCoord(ra*u.deg, dec*u.deg) fields = ['ra', 'dec', 'raErr', 'decErr', 'petroMag_r', 'petroMagErr_r'] # print fields fields = None table = SDSS.query_region(coo, radius=radius*u.arcmin, spectro=False, photoobj_fields=fields) return table
[docs]def get_twomass_catalog(ra=165.86, dec=34.829694, radius=3, catalog='allwise_p3as_psd'): return get_irsa_catalog(ra=ra, dec=dec, radius=radius, catalog='fp_psc', wise=False, twomass=True)
[docs]def get_irsa_catalog(ra=165.86, dec=34.829694, tab=None, radius=3, catalog='allwise_p3as_psd', wise=False, twomass=False, ROW_LIMIT=500000, TIMEOUT=3600): """Query for objects in the `AllWISE <http://wise2.ipac.caltech.edu/docs/release/allwise/>`_ source catalog Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin Returns ------- table : `~astropy.table.Table` Result of the query """ from astroquery.irsa import Irsa Irsa.ROW_LIMIT = ROW_LIMIT Irsa.TIMEOUT = TIMEOUT #all_wise = 'wise_allwise_p3as_psd' #all_wise = 'allwise_p3as_psd' if wise: catalog = 'allwise_p3as_psd' elif twomass: catalog = 'fp_psc' coo = coord.SkyCoord(ra*u.deg, dec*u.deg) table = Irsa.query_region(coo, catalog=catalog, spatial="Cone", radius=radius*u.arcmin, get_query_payload=False) return table
#
[docs]def get_gaia_radec_at_time(gaia_tbl, date=2015.5, format='decimalyear'): """ Use `~astropy.coordinates.SkyCoord.apply_space_motion` to compute GAIA positions at a specific observation date Parameters ---------- gaia_tbl : `~astropy.table.Table` GAIA table query, e.g., provided by `get_gaia_DR2_catalog`. date : e.g., float Observation date that can be parsed with `~astropy.time.Time` format : str Date format, see `~astropy.time.Time.FORMATS`. Returns ------- coords : `~astropy.coordinates.SkyCoord` Projected sky coordinates. """ from astropy.time import Time from astropy.coordinates import SkyCoord # Distance and radial_velocity are dummy numbers needed # to get the space motion correct try: # Try to use pyia import pyia g = pyia.GaiaData(gaia_tbl) coord = g.get_skycoord(distance=1*u.kpc, frame='icrs', radial_velocity=0.*u.km/u.second) except: # From table itself if 'ref_epoch' in gaia_tbl.colnames: ref_epoch = Time(gaia_tbl['ref_epoch'].data, format='decimalyear') else: ref_epoch = Time(2015.5, format='decimalyear') coord = SkyCoord(ra=gaia_tbl['ra'], dec=gaia_tbl['dec'], pm_ra_cosdec=gaia_tbl['pmra'], pm_dec=gaia_tbl['pmdec'], obstime=ref_epoch, frame='icrs', distance=1*u.kpc, radial_velocity=0.*u.km/u.second) new_obstime = Time(date, format=format) coord_at_time = coord.apply_space_motion(new_obstime=new_obstime) return(coord_at_time)
#GAIA_DR2_COLUMNS = None GAIA_DR2_COLUMNS = ['RA_ICRS','DE_ICRS','Epoch','e_RA_ICRS','e_DE_ICRS', 'pmRA','e_pmRA','pmDE','e_pmDE', 'NAL','NAC','Solved','APF','WAL']
[docs]def get_gaia_DR2_vizier_columns(): """ Get translation of Vizier GAIA DR2 columns. The subset of desired output columns can be specified in the GAIA_DR2_COLUMNS global variable, which, if `None`, will be all available columns. """ file = os.path.join(os.path.dirname(__file__), 'data/gaia_dr2_vizier_columns.txt') lines = open(file).readlines()[1:] gdict = OrderedDict() for line in lines: viz_col = line.split()[0] gaia_col = line.split()[1][1:-1] if GAIA_DR2_COLUMNS is not None: if viz_col in GAIA_DR2_COLUMNS: gdict[viz_col] = gaia_col else: gdict[viz_col] = gaia_col return gdict
[docs]def get_gaia_DR2_vizier(ra=165.86, dec=34.829694, radius=3., max=100000, catalog="I/355/gaiadr3", server='vizier.u-strasbg.fr', use_mirror=False, keys=None, mjd=None, clean_mjd=True): """ Query GAIA catalog from Vizier DR2: ``catalog="I/345/gaia2"`` EDR3: ``catalog="I/350/gaiaedr3"`` DR3: ``catalog="I/355/gaiadr3"`` """ from astroquery.vizier import Vizier import astropy.units as u import astropy.coordinates as coord coo = coord.SkyCoord(ra=ra, dec=dec, unit=(u.deg, u.deg), frame='icrs') gdict = get_gaia_DR2_vizier_columns() if keys is None: keys = list(gdict.keys()) try: # Hack, Vizier object doesn't seem to allow getting all keys # simultaneously (astroquery v0.3.7) N = 9 for i in range(len(keys)//N+1): v = Vizier(catalog=catalog, columns=['+_r']+keys[i*N:(i+1)*N]) v.VIZIER_SERVER = server v.ROW_LIMIT = max tab = v.query_region(coo, radius="{0}m".format(radius), catalog=catalog)[0] if i == 0: result = tab else: for k in tab.colnames: #print(i, k) result[k] = tab[k] for k in gdict: if k in result.colnames: result.rename_column(k, gdict[k]) except: utils.log_exception(utils.LOGFILE, traceback) return False if mjd is not None: rd = get_gaia_radec_at_time(result, date=mjd, format='mjd') result['ra_time'] = rd.ra.deg result['dec_time'] = rd.dec.deg result.meta['cootime'] = mjd, 'Specified MJD for ra/dec_time' if clean_mjd: ok = np.isfinite(rd.ra.deg + rd.dec.deg) result = result[ok] return result
def gaia_catalog_for_assoc(assoc_name='j191132p1652_hen-2-427-f335m_00007'): """ Get GAIA catalog """ from grizli.aws import db from grizli import utils res = db.SQL(f"""select t_min, t_max, filter, footprint from assoc_table where assoc_name = '{assoc_name}' """) sr = utils.SRegion(res['footprint'][0]) for fp in res['footprint'][1:]: sri = utils.SRegion(fp) sr.xy.extend(sri.xy) un = sr.union() radius = np.sqrt(un.sky_area()[0]).value gaia = get_gaia_DR2_vizier(*un.centroid[0], radius=radius) rd = get_gaia_radec_at_time(gaia, np.mean(res['t_min']), format='mjd') gaia['ra_orig'] = gaia['ra'] gaia['ra_dec'] = gaia['dec'] gaia['ra'] = rd.ra.deg gaia['dec'] = rd.dec.deg coo = np.array([gaia['ra'], gaia['dec']]) in_footprint = un.path[0].contains_points(coo.T) gaia['in_assoc'] = in_footprint gaia.meta['assocnam'] = assoc_name gaia['valid'] = np.isfinite(gaia['ra'] + gaia['dec']) return gaia
[docs]def gaia_dr2_conesearch_query(ra=165.86, dec=34.829694, radius=3., max=100000): """ Generate a query string for the TAP servers TBD Parameters ---------- ra, dec : float RA, Dec in decimal degrees radius : float Search radius, in arc-minutes. Returns ------- query : str Query string """ query = "SELECT TOP {3} * FROM gaiadr2.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',{0},{1},{2:.2f}))=1".format(ra, dec, radius/60., max) return query
[docs]def get_gaia_DR2_catalog(ra=165.86, dec=34.829694, radius=3., use_mirror=True, max_wait=20, max=100000, output_file='gaia.fits'): """Query GAIA DR2 astrometric catalog Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin use_mirror : bool If True, use the mirror at `gaia.ari.uni-heidelberg.de`. Otherwise use `gea.esac.esa.int`. Returns ------- table : `~astropy.table.Table` Result of the query """ try: import httplib from urllib import urlencode except: # python 3 import http.client as httplib from urllib.parse import urlencode # import http.client in Python 3 # import urllib.parse in Python 3 import time from xml.dom.minidom import parseString host = "gea.esac.esa.int" port = 80 pathinfo = "/tap-server/tap/async" if use_mirror: host = "gaia.ari.uni-heidelberg.de" pathinfo = "/tap/async" # ------------------------------------- # Create job query = gaia_dr2_conesearch_query(ra=ra, dec=dec, radius=radius, max=max) # "SELECT TOP 100000 * FROM gaiadr2.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',{0},{1},{2:.2f}))=1".format(ra, dec, radius/60.) print(query) params = urlencode({ "REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "fits", "PHASE": "RUN", "QUERY": query }) headers = { "Content-type": "application/x-www-form-urlencoded", "Accept": "text/plain" } connection = httplib.HTTPConnection(host, port) connection.request("POST", pathinfo, params, headers) # Status response = connection.getresponse() print("Status: " + str(response.status), "Reason: " + str(response.reason)) # Server job location (URL) location = response.getheader("location") print("Location: " + location) # Jobid jobid = location[location.rfind('/')+1:] print("Job id: " + jobid) connection.close() # ------------------------------------- # Check job status, wait until finished tcount = 0 while True: connection = httplib.HTTPConnection(host, port) connection.request("GET", pathinfo+"/"+jobid) response = connection.getresponse() data = response.read() # XML response: parse it to obtain the current status dom = parseString(data) if use_mirror: phaseElement = dom.getElementsByTagName('phase')[0] else: phaseElement = dom.getElementsByTagName('uws:phase')[0] phaseValueElement = phaseElement.firstChild phase = phaseValueElement.toxml() print("Status: " + phase) # Check finished if phase == 'COMPLETED': break #wait and repeat time.sleep(0.2) tcount += 0.2 if (phase == 'ERROR') | (tcount > max_wait): return False # print "Data:" # print data connection.close() # ------------------------------------- # Get results connection = httplib.HTTPConnection(host, port) connection.request("GET", pathinfo+"/"+jobid+"/results/result") response = connection.getresponse() data = response.read() outputFileName = output_file + (not use_mirror)*".gz" try: outputFile = open(outputFileName, "w") outputFile.write(data) except: # Python 3 outputFile = open(outputFileName, "wb") outputFile.write(data) outputFile.close() connection.close() print("Data saved in: " + outputFileName) if not use_mirror: # ESA archive returns gzipped try: os.remove(output_file) except: pass os.system('gunzip {output_file}.gz'.format(output_file=output_file)) table = Table.read(output_file, format='fits') return table
[docs]def gen_tap_box_query(ra=165.86, dec=34.829694, radius=3., corners=None, max=100000, db='ls_dr7.tractor_primary', columns=['*'], rd_colnames=['ra', 'dec'], wcs_pad=0.5): """ Generate a query string for the NOAO Legacy Survey TAP server Parameters ---------- ra, dec : float RA, Dec in decimal degrees radius : float Search radius, in arc-minutes. corners : 4-tuple, `~astropy.wcs.WCS` or None ra_min, ra_max, dec_min, dec_max of a query box to use instead of `radius`. Or if a `~astropy.wcs.WCS` object, get limits from the `~astropy.wcs.WCS.calc_footprint` method Returns ------- query : str Query string """ rmi = radius/60/2 cosd = np.cos(dec/180*np.pi) if max is not None: maxsel = 'TOP {0}'.format(max) else: maxsel = '' if corners is not None: if hasattr(corners, 'calc_footprint'): foot = corners.calc_footprint() left = foot[:,0].min() right = foot[:,0].max() bottom = foot[:,1].min() top = foot[:,1].max() dx = (right-left) dy = (top-bottom) left -= wcs_pad*dx right += wcs_pad*dx bottom -= wcs_pad*dy top += wcs_pad*dy elif len(corners) != 4: msg = 'corners needs 4 values (ra_min, ra_max, dec_min, dec_max)' raise ValueError(msg) else: left, right, bottom, top = corners else: left = ra - rmi / cosd right = ra + rmi / cosd bottom = dec - rmi top = dec + rmi fmt = dict(rc=rd_colnames[0], dc=rd_colnames[1], left=left, right=right, top=top, bottom=bottom, maxsel=maxsel, db=db, output_columns=', '.join(columns)) if not np.isfinite(ra+dec): query = "SELECT {maxsel} {output_columns} FROM {db} " else: query = ("SELECT {maxsel} {output_columns} FROM {db} WHERE " + "{rc} > {left} AND {rc} < {right} AND " + "{dc} > {bottom} AND {dc} < {top} ") return query.format(**fmt)
[docs]def query_tap_catalog(ra=165.86, dec=34.829694, radius=3., corners=None, max_wait=20, db='ls_dr9.tractor', columns=['*'], extra='', rd_colnames=['ra', 'dec'], tap_url='https://datalab.noirlab.edu/tap', max=1000000, clean_xml=True, verbose=True, des=False, gaia=False, nsc=False, vizier=False, skymapper=False, hubble_source_catalog=False, tap_kwargs={}, **kwargs): """Query NOAO Catalog holdings Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin corners : 4-tuple, `~astropy.wcs.WCS` or None ra_min, ra_max, dec_min, dec_max of a query box to use instead of `radius`. Or if a `WCS` object, get limits from the `~astropy.wcs.WCS.calc_footprint` method db : str Parent database (https://datalab.noirlab.edu/query.php). columns : list of str List of columns to output. Default ['*'] returns all columns. extra : str String to add to the end of the positional box query, e.g., 'AND mag_auto_i > 16 AND mag_auto_i < 16.5'. rd_colnames : str, str Column names in `db` corresponding to ra/dec (degrees). tap_url : str TAP hostname des : bool Query `des_dr1.main` from NOAO. gaia : bool Query `gaiadr2.gaia_source` from http://gea.esac.esa.int. nsc : bool Query the NOAO Source Catalog (Nidever et al. 2018), `nsc_dr1.object`. vizier : bool Use the VizieR TAP server at http://tapvizier.u-strasbg.fr/TAPVizieR/tap, see http://tapvizier.u-strasbg.fr/adql/about.html. hubble_source_catalog : bool Query the Hubble Source Catalog (v3). If no 'NumImages' criteria is found in `extra`, then add an additional requirement: >>> extra += 'AND NumImages > 1' Returns ------- table : `~astropy.table.Table` Result of the query """ from astroquery.utils.tap.core import TapPlus # DES DR1 if des: if verbose: print('Query DES DR1 from NOAO') db = 'des_dr1.main' tap_url = 'https://datalab.noirlab.edu/tap' # NOAO source catalog, seems to have some junk if nsc: if verbose: print('Query NOAO source catalog') db = 'nsc_dr1.object' tap_url = 'https://datalab.noirlab.edu/tap' extra += ' AND nsc_dr1.object.flags = 0' # GAIA DR2 if gaia: if verbose: print('Query GAIA DR2 from ESA') db = 'gaiadr2.gaia_source' tap_url = 'http://gea.esac.esa.int/tap-server/tap' # VizieR TAP server if vizier: if verbose: print('Query {0} from VizieR TAP server'.format(db)) tap_url = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap' rd_colnames = ['RAJ2000', 'DEJ2000'] if skymapper: if verbose: print('Query {0} from VizieR TAP server'.format(db)) tap_url = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap' rd_colnames = ['RAICRS', 'DEICRS'] if hubble_source_catalog: if db is None: db = 'dbo.SumPropMagAutoCat' elif 'dbo' not in db: db = 'dbo.SumPropMagAutoCat' tap_url = 'http://vao.stsci.edu/HSCTAP/tapservice.aspx' rd_colnames = ['MatchRA', 'MatchDec'] if 'NumImages' not in extra: extra += 'AND NumImages > 1' tap = TapPlus(url=tap_url, **tap_kwargs) query = gen_tap_box_query(ra=ra, dec=dec, radius=radius, max=max, db=db, columns=columns, rd_colnames=rd_colnames, corners=corners) job = tap.launch_job(query+extra, dump_to_file=True, verbose=verbose) try: table = job.get_results() if clean_xml: if hasattr(job, 'outputFile'): jobFile = job.outputFile else: jobFile = job.get_output_file() os.remove(jobFile) # Provide ra/dec columns for c, cc in zip(rd_colnames, ['ra', 'dec']): if (c in table.colnames) & (cc not in table.colnames): table[cc] = table[c] table.meta['TAPURL'] = tap_url, 'TAP URL' table.meta['TAPDB'] = db, 'TAP database name' table.meta['TAPQUERY'] = query+extra, 'TAP query' table.meta['RAQUERY'] = ra, 'Query central RA' table.meta['DECQUERY'] = dec, 'Query central Dec' table.meta['RQUERY'] = radius, 'Query radius, arcmin' if hubble_source_catalog: for col in table.colnames: if table[col].dtype == 'object': print('Reformat column: {0}'.format(col)) strcol = list(table[col]) table[col] = strcol except: if hasattr(job, 'outputFile'): jobFile = job.outputFile else: jobFile = job.get_output_file() print('Query failed, check {0} for error messages'.format(jobFile)) table = None return table
# Limit Hubble Source Catalog query to brighter sources in limited bands HSCv3_FILTER_LIMITS = {'W3_F160W': 23.5, 'W3_F140W': 23.5, 'W3_F125W': 23.5, 'W3_F110W': 23.5, 'W3_F098M': 23.5, 'W3_F105W': 23.5, 'A_F814W': 23.5, 'W3_F814W': 23.5, 'A_F606W': 23.5, 'W3_F606W': 23.5, 'A_F850LP': 23.5, 'W3_F850LP': 23.5, 'A_F775W': 23.5, 'W3_F775W': 23.5} HSCv3_COLUMNS = ['MatchRA', 'MatchDec', 'CI', 'CI_Sigma', 'KronRadius', 'KronRadius_Sigma', 'Extinction', 'TargetName', 'NumImages', 'NumFilters', 'NumVisits', 'DSigma']
[docs]def get_hubble_source_catalog(ra=0., dec=0., radius=3, corners=None, max=int(1e7), extra=' AND NumImages > 0', kron_max=0.45, dsigma_max=100, clip_singles=10*u.arcsec, verbose=True, columns=HSCv3_COLUMNS, filter_limits=HSCv3_FILTER_LIMITS): """ Query NOAO Source Catalog, which is aligned to GAIA DR1. The default `extra` query returns well-detected sources in red bands. filter_limits : query on individual HSC filter magnitudes """ import astropy.table msg = 'Query NOAO Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')' print(msg.format(ra=ra, dec=dec, radius=radius)) if kron_max is not None: extra += ' AND KronRadius < {0}'.format(kron_max) if dsigma_max is not None: extra += ' AND DSigma < {0}'.format(dsigma_max) if filter_limits is not None: limit_list = ['{0} < {1}'.format(f, filter_limits[f]) for f in filter_limits] filter_selection = ' AND ({0})'.format(' OR '.join(limit_list)) extra += filter_selection columns += [f for f in filter_limits] db = 'dbo.SumPropMagAutoCat p join dbo.SumMagAutoCat m on p.MatchID = m.MatchID' else: db = 'dbo.SumPropMagAutoCat' tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners, max=max, extra=extra, hubble_source_catalog=True, verbose=verbose, db=db, columns=columns) if clip_singles not in [None, False]: rr = tab['NumImages'] > 1 if (rr.sum() > 0) & ((~rr).sum() > 0): r0, r1 = tab[rr], tab[~rr] idx, dr = utils.GTable(r0).match_to_catalog_sky(r1) new = dr > clip_singles xtab = astropy.table.vstack([r0, r1[new]]) if verbose: msg = ('HSCv3: Remove {0} NumImages == 1 sources ' + 'with tolerance {1}') print(msg.format((~new).sum(), clip_singles)) return tab
[docs]def get_nsc_catalog(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (rerr < 0.08 OR ierr < 0.08 OR zerr < 0.08) AND raerr < 0.2 AND decerr < 0.2', verbose=True): """ Query NOAO Source Catalog, which is aligned to GAIA DR1. The default `extra` query returns well-detected sources in red bands. """ msg = 'Query NOAO Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')' print(msg.format(ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners, extra=extra, nsc=True, verbose=verbose, max=max) return tab
[docs]def get_desdr1_catalog(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (e_rmag < 0.15 OR e_imag < 0.15)', verbose=True): """ Query DES DR1 Catalog from Vizier The default `extra` query returns well-detected sources in one or more red bands. """ msg = 'Query DES Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')' print(msg.format(ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners, extra=extra, db='"II/357/des_dr1"', vizier=True, verbose=verbose, max=max) return tab
def get_desdr1_catalog_old(ra=0., dec=0., radius=3, corners=None, max=100000, extra=' AND (magerr_auto_r < 0.15 OR magerr_auto_i < 0.15)', verbose=True): """ Query DES DR1 Catalog. The default `extra` query returns well-detected sources in one or more red bands. """ msg = 'Query DES Source Catalog ({ra:.5f},{dec:.5f},{radius:.1f}\')' print(msg.format(ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, corners=corners, extra=extra, des=True, verbose=verbose, max=max) return tab
[docs]def get_vhs_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='', table='vhs_dr4'): """ VEXAS DR2 from vizier https://vizier.cds.unistra.fr/viz-bin/VizieR?-source=II/369 table: 'vexasds', 'vexasps', 'vexassm' """ msg = 'Query VHS DR4 {table} catalog ({ra},{dec},{radius})' print(msg.format(table=table, ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2, corners=corners, extra=extra, vizier=True, db=f'"II/359/{table}"', verbose=verbose, max=max_records) return tab
[docs]def get_vexas_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='', table='vexasdes'): """ VEXAS DR2 from vizier https://vizier.cds.unistra.fr/viz-bin/VizieR?-source=II/369 table: 'vexasds', 'vexasps', 'vexassm' """ msg = 'Query VEXAS DR2 {table} catalog ({ra},{dec},{radius})' print(msg.format(table=table, ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2, corners=corners, extra=extra, vizier=True, db=f'"II/369/{table}"', verbose=verbose, max=max_records) return tab
[docs]def get_skymapper_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra=''): """ Get Skymapper DR1 from Vizier """ msg = 'Query Skymapper DR1 catalog ({ra},{dec},{radius})' print(msg.format(ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2, corners=corners, extra=extra, skymapper=True, db='"II/358/smss"', verbose=verbose, max=max_records) return tab
def get_legacysurveys_catalog(ra=0., dec=0., radius=3., verbose=True, db='ls_dr9.tractor', sn_lim=('r',10), **kwargs): """ Query LegacySurveys TAP catalog """ if verbose: msg = 'Query LegacySurveys ({db}) catalog ({ra},{dec},{radius:.2f})' print(msg.format(ra=ra, dec=dec, radius=radius, db=db)) if (sn_lim is not None): if len(sn_lim) >= 2: b = sn_lim[0] extra = f' AND flux_{b}*SQRT(flux_ivar_{b}) > {sn_lim[1]}' if 'extra' in kwargs: kwargs['extra'] += extra else: kwargs['extra'] = extra tab = query_tap_catalog(ra=ra, dec=dec, radius=radius, db=db, **kwargs) # if (sn_lim is not None) & (len(tab) > 1): # if len(sn_lim) >= 2: # sn = tab[f'flux_{sn_lim[0]}'] # sn *= np.sqrt(tab[f'flux_ivar_{sn_lim[0]}']) # # valid = sn > sn_lim[1] # return tab[valid] return tab
[docs]def get_panstarrs_catalog(ra=0., dec=0., radius=3., corners=None, max_records=500000, verbose=True, extra='AND "II/349/ps1".e_imag < 0.2 AND "II/349/ps1".e_RAJ2000 < 0.15 AND "II/349/ps1".e_DEJ2000 < 0.15'): """ Get PS1 from Vizier """ msg = 'Query PanSTARRS catalog ({ra},{dec},{radius})' print(msg.format(ra=ra, dec=dec, radius=radius)) tab = query_tap_catalog(ra=ra, dec=dec, radius=radius*2, corners=corners, extra=extra, vizier=True, db='"II/349/ps1"', verbose=verbose, max=max_records) return tab
[docs]def get_radec_catalog(ra=0., dec=0., radius=3., product='cat', verbose=True, reference_catalogs=['GAIA', 'LS_DR9', 'PS1', 'Hubble', 'NSC', 'SDSS', 'WISE', 'DES'], use_self_catalog=False, **kwargs): """Decide what reference astrometric catalog to use First search SDSS, then WISE looking for nearby matches. Parameters ---------- ra, dec : float Center of the query region, decimal degrees radius : float Radius of the query, in arcmin product : str Basename of the drizzled product. If a locally-created catalog with filename that startswith `product` is found, use that one instead of the external (low precision) catalogs so that you're matching HST-to-HST astrometry. reference_catalogs : list Order in which to query reference catalogs. Options are 'GAIA', 'PS1' (STScI PanSTARRS), 'SDSS', 'WISE', 'NSC' (NOAO Source Catalog), 'DES' (Dark Energy Survey DR1), 'Hubble' (Hubble Source Catalog v3), 'LS_DR9' (LegacySurveys DR9). Returns ------- radec : str Filename of the RA/Dec list derived from the parent catalog ref_catalog : str, {'SDSS', 'WISE', 'VISIT'} Provenance of the `radec` list. """ query_functions = {'SDSS': get_sdss_catalog, 'GAIA_TAP': get_gaia_DR2_catalog, 'PS1': get_panstarrs_catalog, 'WISE': get_irsa_catalog, '2MASS': get_twomass_catalog, 'GAIA_Vizier': get_gaia_DR2_vizier, 'GAIA': get_gaia_DR2_vizier, 'NSC': get_nsc_catalog, 'DES': get_desdr1_catalog, 'Hubble': get_hubble_source_catalog, 'Skymapper': get_skymapper_catalog, 'VEXAS': get_vexas_catalog, 'LS_DR9': get_legacysurveys_catalog} # Try queries has_catalog = False ref_catalog = 'None' ref_cat = [] for ref_src in reference_catalogs: try: if ref_src == 'GAIA': try: ref_cat = query_functions[ref_src](ra=ra, dec=dec, radius=radius, use_mirror=False) except: try: ref_cat = query_functions[ref_src](ra=ra, dec=dec, radius=radius) except: ref_cat = False # Try GAIA mirror at Heidelberg if ref_cat is False: ref_cat = query_functions[ref_src](ra=ra, dec=dec, radius=radius, use_mirror=True) else: ref_cat = query_functions[ref_src](ra=ra, dec=dec, radius=radius) # # # ref_cat = query_functions[ref_src](ra=ra, dec=dec, # radius=radius) valid = np.isfinite(ref_cat['ra']+ref_cat['dec']) ref_cat = ref_cat[valid] if len(ref_cat) < 2: raise ValueError table_to_regions(ref_cat, output='{0}_{1}.reg'.format(product, ref_src.lower())) ref_cat['ra', 'dec'].write('{0}_{1}.radec'.format(product, ref_src.lower()), format='ascii.commented_header', overwrite=True) radec = '{0}_{1}.radec'.format(product, ref_src.lower()) ref_catalog = ref_src has_catalog = True if len(ref_cat) > 0: break except: print('{0} query failed'.format(ref_src)) has_catalog = False if (ref_src.startswith('GAIA')) & ('date' in kwargs) & has_catalog: if kwargs['date'] is not None: if 'date_format' in kwargs: date_format = kwargs['date_format'] else: date_format = 'mjd' gaia_tbl = ref_cat # utils.GTable.gread('gaia.fits') coo = get_gaia_radec_at_time(gaia_tbl, date=kwargs['date'], format=date_format) coo_tbl = utils.GTable() coo_tbl['ra'] = coo.ra coo_tbl['dec'] = coo.dec ok = np.isfinite(coo_tbl['ra']) & np.isfinite(coo_tbl['dec']) coo_tbl.meta['date'] = kwargs['date'] coo_tbl.meta['datefmt'] = date_format msg = 'Apply observation ({0},{1}) to GAIA catalog' print(msg.format(kwargs['date'], date_format)) table_to_regions(coo_tbl[ok], output='{0}_{1}.reg'.format(product, ref_src.lower())) coo_tbl['ra', 'dec'][ok].write('{0}_{1}.radec'.format(product, ref_src.lower()), format='ascii.commented_header', overwrite=True) if not has_catalog: return False # WISP, check if a catalog already exists for a given rootname and use # that if so. cat_files = glob.glob('-f1'.join(product.split('-f1')[:-1]) + '-f*.cat*') if (len(cat_files) > 0) & (use_self_catalog): ref_cat = utils.GTable.gread(cat_files[0]) root = cat_files[0].split('.cat')[0] ref_cat['X_WORLD', 'Y_WORLD'].write('{0}.radec'.format(root), format='ascii.commented_header', overwrite=True) radec = '{0}.radec'.format(root) ref_catalog = 'VISIT' if verbose: msg = '{0} - Reference RADEC: {1} [{2}] N={3}' print(msg.format(product, radec, ref_catalog, len(ref_cat))) return radec, ref_catalog