#!/usr/bin/env python # -*- coding: utf-8 -*- # # query_lsd.py # # Copyright 2012 Greg # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, # MA 02110-1301, USA. # # import os, sys, argparse from os.path import abspath import matplotlib as mplib mplib.use('Agg') import matplotlib.pyplot as plt import healpy as hp import numpy as np np.seterr(all='ignore') try: import astropy.io.fits as pyfits except ImportError: import pyfits import h5py import lsd import iterators import hputils from astropy.table import Table def lb2pix(nside, l, b, nest=True): theta = np.pi/180. * (90. - b) phi = np.pi/180. * l return hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest) def pix2lb(nside, ipix, nest=True): theta, phi = hp.pixelfunc.pix2ang(nside, ipix, nest=True) l = 180./np.pi * phi b = 90. - 180./np.pi * theta return l, b def PS1_saturated(obj): ''' Return indices of saturated detections in bands. ''' sat = np.zeros((len(obj), 5), dtype=np.bool) PS1_sat_limit = [14.5, 14.5, 14.5, 14., 13.] for i,m in enumerate(PS1_sat_limit): idx = (obj['ps1mag'][:,i] < m) sat[idx, i] = 1 return sat def UKIDSS_saturated(ppErrBits): ''' Return indices of saturated detections in bands. ''' sat = np.zeros((len(ppErrBits), 3), dtype=np.bool) ukidss_sat_limit = [256,256,256] for i,m in enumerate(ukidss_sat_limit): idx = (ppErrBits[:,i] > m) sat[idx, i] = 1 return sat def tmass_ext(obj): ''' Returns indices of extended objects. ''' return (obj['ext_key'] > 0) def tmass_hq_phot(obj): ''' Return index of which detections in which bands are of high quality, as per the 2MASS recommendations: ''' # Photometric quality in each passband idx = (obj['ph_qual'] == '0') obj['ph_qual'][idx] = '000' ph_qual = np.array(map(list, obj['ph_qual'])) # Read flag in each passband idx = (obj['rd_flg'] == '0') obj['rd_flg'][idx] = '000' rd_flg = np.array(map(list, obj['rd_flg'])) #rd_flg = (rd_flg == '1') | (rd_flg == '3') # Contamination flag in each passband idx = (obj['cc_flg'] == '0') obj['cc_flg'][idx] = '000' cc_flg = np.array(map(list, obj['cc_flg'])) # Combine passband flags cond_1 = (ph_qual == 'A') | (rd_flg == '1') | (rd_flg == '3') cond_1 &= (cc_flg == '0') # Source quality flags cond_2 = (obj['use_src'] == 1) & (obj['gal_contam'] == 0)# & (obj['ext_key'] <= 0) # Combine all flags for each object hq = np.empty((len(obj), 3), dtype=np.bool) for i in range(3): hq[:,i] = cond_1[:,i] & cond_2 return hq def start_file(base_fname, index): fout = open('%s_%d.in' % (base_fname, index), 'wb') f.write(np.array([0], dtype=np.uint32).tostring()) return f def to_file(f, pix_index, nside, nest, EBV, data): close_file = False if type(f) == str: f = h5py.File(fname, 'a') close_file = True ds_name = '/photometry/pixel %d-%d' % (nside, pix_index) ds = f.create_dataset(ds_name, data=data, chunks=True, compression='gzip', compression_opts=3) ds[:] = data[:] N_stars = data.shape[0] t,p = hp.pixelfunc.pix2ang(nside, pix_index, nest=nest) t *= 180. / np.pi p *= 180. / np.pi gal_lb = np.array([p, 90. - t], dtype='f8') att_f8 = np.array([EBV], dtype='f8') att_u8 = np.array([pix_index], dtype='u8') att_u4 = np.array([nside, N_stars], dtype='u4') att_u1 = np.array([nest], dtype='u1') ds.attrs['healpix_index'] = att_u8[0] ds.attrs['nested'] = att_u1[0] ds.attrs['nside'] = att_u4[0] ds.attrs['l'] = gal_lb[0] ds.attrs['b'] = gal_lb[1] ds.attrs['EBV'] = att_f8[0] if close_file: f.close() return gal_lb def adaptive_subdivide(pix_idx, nside, obj, n_stars_max, n_stars_min=10, nside_max=2048): n_obj = len(obj) n_max = n_stars_max[nside] n_max = max(n_max, 12*n_stars_min) # print 'n_max: {} -> {} (EBV = {:.3f} mag)'.format(n_stars_max[nside], n_max, np.median(obj['EBV'])) # Subdivide pixel if (n_obj > n_max) and (nside < nside_max): sub_pix_idx = lb2pix(nside*2, obj['l'], obj['b'], nest=True) # Check that all pixels have more than minimum # of pixels ''' over_threshold = True for i in xrange(4 * pix_idx, 4 * pix_idx + 4): idx = (sub_pix_idx == i) if np.sum(idx) < n_stars_min: over_threshold = False break if not over_threshold: return [(nside, pix_idx, obj)] ''' # Return subdivided pixel ret = [] for i in xrange(4 * pix_idx, 4 * pix_idx + 4): idx = (sub_pix_idx == i) tmp = adaptive_subdivide(i, nside*2, obj[idx], n_stars_max, n_stars_min, nside_max, scale_by_opacity=scale_by_opacity) for pix in tmp: ret.append(pix) return ret else: return [(nside, pix_idx, obj)] def subdivider(keyvalue, nside, n_stars_max, n_stars_min, nside_max): pix_index, obj = keyvalue obj = lsd.colgroup.fromiter(obj, blocks=True) # Adaptively subdivide pixel ret = adaptive_subdivide(pix_index, nside, obj, n_stars_max, n_stars_min, nside_max) for subpixel in ret: sub_nside, sub_idx, sub_obj = subpixel yield ((sub_nside, sub_idx), sub_obj) def mapper(qresult, nside, nest, bounds): obj = lsd.colgroup.fromiter(qresult, blocks=True) if (obj is not None) and (len(obj) > 0): # Determine healpix index of each star theta = np.pi/180. * (90. - obj['b']) phi = np.pi/180. * obj['l'] pix_indices = hp.ang2pix(nside, theta, phi, nest=nest) # Group together stars having same index for pix_index, block_indices in iterators.index_by_key(pix_indices): yield (pix_index, obj[block_indices]) def reducer(keyvalue): pix_index, obj = keyvalue obj = lsd.colgroup.fromiter(obj, blocks=True) print("unwise secondary {}".format(np.sum(obj['unwise_obj_primary.primary']!=1))) #print('objlength',len(obj)) #set up data structure data = np.empty(len(obj), dtype=[('ra','f8'), ('dec','f8'),('EBV','f4'),('obj_id','u8'), ('l','f8'),('b','f8'), ('ps1mag','5f4'),('ps1magerr','5f4'),('ps1_nmag_ok','5i4'), ('tmass.designation','a16'),('tmassmag','3f4'),('tmassmagerr','3f4'), ('gaia.source_id','u8'), ('parallax','f8'), ('parallax_error','f8'), ('pmra','f8'),('pmra_error','f8'),('pmdec','f8'),('pmdec_error','f8'), ('visibility_periods_used','u4'),('astrometric_excess_noise','f4'), ('astrometric_sigma5d_max','f4'),('astrometric_excess_noise_sig','f4'), ('phot_bp_rp_excess_factor','f4'),('gaiamag','3f4'),('gaiamagerr','3f4'), ('unwise.objid','a16'),('unwisemag','2f4'),('unwisemagerr','2f4'),('unwise.fracflux','2f4'), ('unwise.flags_info','2i4'),('unwise.primary','i4'), ('ukidssmag','3f4'),('ukidssmagerr','3f4')]) tmassmag=np.vstack((obj['J'],obj['H'],obj['K'])).T tmassmagerr=np.vstack((obj['J_sig'],obj['H_sig'],obj['K_sig'])).T ukidssmag=np.vstack((obj['ukidss.j_1AperMag3'],obj['ukidss.hAperMag3'],obj['ukidss.kAperMag3'])).T ukidssmagerr=np.vstack((obj['ukidss.j_1AperMag3Err'],obj['ukidss.hAperMag3Err'],obj['ukidss.kAperMag3Err'])).T unwisemag=np.vstack((obj['W1'],obj['W2'])).T unwisemagerr=np.vstack((obj['W1_err'],obj['W2_err'])).T gaiamag=np.vstack((obj['gaia.phot_g_mean_mag'],obj['gaia.phot_bp_mean_mag'],obj['gaia.phot_rp_mean_mag'])).T gaiamagerr=np.vstack((obj['phot_g_mean_mag_err'],obj['phot_bp_mean_mag_err'],obj['phot_rp_mean_mag_err'])).T ## 2MASS Data ## good_det_tmass = np.empty((len(obj), 3), dtype=np.bool) hq_idx_tmass = (tmassmagerr < 0.2) & (tmassmag > 0) good_det_tmass[:,:] = hq_idx_tmass tmassmag[~good_det_tmass]=np.nan tmassmagerr[~good_det_tmass]=np.nan tmassmag[tmass_ext(obj)]=np.nan tmassmagerr[tmass_ext(obj)]=np.nan ## unWISE Data ## good_det_unwise = np.empty((len(obj), 2), dtype=np.bool) hq_index_unwise=(unwisemagerr < 0.2) & (obj['unwise_obj_primary.flags_unwise']==0) & (obj['unwise_obj_primary.fracflux'] > 0.5) & (unwisemag < 50.) & (unwisemag > 0) good_det_unwise[:,:] = hq_index_unwise unwisemag[~good_det_unwise]=np.nan unwisemagerr[~good_det_unwise]=np.nan #UKIDSS good detections# ukidssppErrBits=np.vstack((obj['ukidss.j_1ppErrBits'],obj['ukidss.hppErrBits'],obj['ukidss.kppErrBits'])).T good_det_ukidss = np.empty((len(obj), 3), dtype=np.bool) hq_idx_ukidss = (ukidssmagerr < 0.2) & (ukidssmag > 0) & ~UKIDSS_saturated(ukidssppErrBits) good_det_ukidss[:,:] = hq_idx_ukidss ukidssmag[~good_det_ukidss]=np.nan ukidssmagerr[~good_det_ukidss]=np.nan goodstar=(obj['ukidss.pStar'] > 0.9) ukidssmag[~goodstar]=np.nan ukidssmagerr[~goodstar]=np.nan data['ra'][:] = obj['ra'][:] data['dec'][:] = obj['dec'][:] data['l'][:]=obj['l'][:] data['b'][:]=obj['b'][:] data['EBV'][:] = obj['EBV'][:] data['obj_id'][:] = obj['obj_id'][:] data['ps1mag'][:,:] = obj['mean'][:,:] data['ps1magerr'][:,:] = obj['err'][:,:] data['ps1_nmag_ok'][:,:] = obj['nmag_ok'][:,:] data['tmass.designation'][:]=obj['tmass.designation'][:] data['tmassmag'][:,:] = tmassmag data['tmassmagerr'][:,:] = tmassmagerr data['gaia.source_id'][:] = obj['gaia.source_id'][:] data['parallax'][:] = obj['gaia.parallax'][:] data['parallax_error'][:] = obj['gaia.parallax_error'][:] data['pmra'][:] = obj['gaia.pmra'][:] data['pmra_error'][:] = obj['gaia.pmra_error'][:] data['pmdec'][:] = obj['gaia.pmdec'][:] data['pmdec_error'][:] = obj['gaia.pmdec_error'][:] data['visibility_periods_used'][:] = obj['gaia.visibility_periods_used'][:] data['astrometric_excess_noise'][:] = obj['gaia.astrometric_excess_noise'][:] data['astrometric_sigma5d_max'][:] = obj['gaia.astrometric_sigma5d_max'][:] data['astrometric_excess_noise_sig'][:] = obj['gaia.astrometric_excess_noise_sig'][:] data['phot_bp_rp_excess_factor'][:] = obj['gaia.phot_bp_rp_excess_factor'][:] data['gaiamag'][:,:] = gaiamag data['gaiamagerr'][:,:] = gaiamagerr data['unwise.objid'][:]=obj['unwise_obj_primary.unwise_objid'][:] data['unwisemag'][:] = unwisemag data['unwisemagerr'][:,:] = unwisemagerr data['unwise.fracflux'][:,:]=obj['unwise_obj_primary.fracflux'][:,:] data['unwise.flags_info'][:,:] = obj['unwise_obj_primary.flags_info'][:,:] data['ukidssmag'][:] = ukidssmag data['ukidssmagerr'][:,:] = ukidssmagerr data['unwise.primary'][:]=obj['unwise_obj_primary.primary'][:] ### PS1 Data ### #Scale PS1 errors err_scale = 1.3 err_floor = 0.02 data['ps1magerr'] = np.sqrt((err_scale * data['ps1magerr'])**2. + err_floor**2.) # Determine which bands have good detections good_det_ps1 = np.empty((len(data), 5), dtype=np.bool) good_det_ps1[:,:] = (data['ps1_nmag_ok'] > 0) & ~PS1_saturated(data) & (data['ps1magerr'] < 0.20) #print(data['ps1mag']) #print(~PS1_saturated(data)) #print(np.sum(data['ps1_nmag_ok'] > 0),np.sum(~PS1_saturated(data)),np.sum((data['ps1err'] < 0.20))) data['ps1mag'][~good_det_ps1]=np.nan data['ps1mag'][~good_det_ps1]=np.nan # Filter out stars that have detections in <4 bands (double check) ps1_mask = (np.sum(good_det_ps1, axis=1) >= 4) #print('ps1masksum',np.sum(ps1_mask)) data=data[ps1_mask] #Gaia Data### gaia_mask = ((np.isfinite(data['parallax'])) & (data['parallax_error'] != 0) & (data['gaiamag'][:,0] <= 21) & (data['visibility_periods_used'] >= 6) & ( data['astrometric_sigma5d_max'] <= 1.2 * np.max([np.ones((len(data))), 10**(0.2 * (data['gaiamag'][:,0] - 18))], axis=0))) data=data[gaia_mask] #print('gaiamasksum',np.sum(ps1_mask)) #print('datalength',len(data)) yield (pix_index, data) def main(): parser = argparse.ArgumentParser( prog='query_brutus_final.py', description='Generate input files for brutus runs', add_help=True) parser.add_argument('out', type=str, help='Output filename.') parser.add_argument('-nmin', '--nside-min', type=int, default=256, help='Lowest resolution in healpix nside parameter (default: 512).') parser.add_argument('-nmax', '--nside-max', type=int, default=256, help='Lowest resolution in healpix nside parameter (default: 512).') parser.add_argument('-rt', '--res-thresh', type=int, nargs='+', default=None, help='Maximum # of pixels for each healpix resolution (from lowest to highest).') parser.add_argument('-b', '--bounds', type=float, nargs=4, default=None, help='Restrict pixels to region enclosed by: l_min, l_max, b_min, b_max.') parser.add_argument('-min', '--min-stars', type=int, default=1, help='Minimum # of stars in pixel (default: 1).') parser.add_argument('-max', '--max-stars', type=int, default=50000, help='Maximum # of stars in file') parser.add_argument('-w', '--n-workers', type=int, default=5, help='# of workers for LSD to use.') values = parser.parse_args() #nPointlike = values.n_bands - 1 #if nPointlike == 0: # nPointlike = 1 # Handle adaptive pixelization parameters base_2_choices = [2**n for n in xrange(15)] if values.nside_min not in base_2_choices: raise ValueError('--nside-min is not a small power of two.') elif values.nside_max not in base_2_choices: raise ValueError('--nside-max is not a small power of two.') elif values.nside_max < values.nside_min: raise ValueError('--nside-max is less than --nside-min.') n_stars_max = None n_pixels_at_res = None nside_options = None if values.nside_min == values.nside_max: n_stars_max = {values.nside_max: 1} n_pixels_at_res = {values.nside_max: 0} nside_options = [values.nside_max] else: n_stars_max = {} n_pixels_at_res = {} nside_options = [] nside = values.nside_min k = 0 while nside < values.nside_max: n_stars_max[nside] = values.res_thresh[k] n_pixels_at_res[nside] = 0 nside_options.append(nside) nside *= 2 k += 1 n_stars_max[values.nside_max] = 1 n_pixels_at_res[values.nside_max] = 0 nside_options.append(values.nside_max) # Determine the query bounds query_bounds = None if values.bounds is not None: pix_scale = hp.pixelfunc.nside2resol(values.nside_min) * 180. / np.pi query_bounds = [] query_bounds.append(values.bounds[0] - 3.*pix_scale) #max([0., values.bounds[0] - 3.*pix_scale])) query_bounds.append(values.bounds[1] + 3.*pix_scale) #min([360., values.bounds[1] + 3.*pix_scale])) query_bounds.append(max([-90., values.bounds[2] - 3.*pix_scale])) query_bounds.append(min([90., values.bounds[3] + 3.*pix_scale])) query_bounds = lsd.bounds.rectangle(query_bounds[0], query_bounds[2], query_bounds[1], query_bounds[3], coordsys='gal') #query_bounds = (query_bounds, []) query_bounds = lsd.bounds.make_canonical(query_bounds) # Set up the query db = lsd.DB(os.environ['LSD_DB']) query = None # query = ("select ra, dec, equgal(ra, dec) as (l, b), SFD.EBV(l, b) as EBV, obj_id, 2.5/log(10.)*err/mean as err, -2.5*log10(mean) as mean, -2.5*log10(mean_ap) as mean_ap, nmag_ok, " # "tmass.designation, tmass.ph_qual as ph_qual, tmass.use_src as use_src, tmass.rd_flg as rd_flg, " # "tmass.ext_key as ext_key, tmass.gal_contam as gal_contam, tmass.cc_flg as cc_flg, " # "tmass.j_m as J, tmass.j_msigcom as J_sig, tmass.h_m as H, tmass.h_msigcom as H_sig, tmass.k_m as K, tmass.k_msigcom as K_sig, " # "gaia.source_id, gaia.parallax, gaia.parallax_error, gaia.pmra, gaia.pmra_error, gaia.pmdec, gaia.pmdec_error, gaia.visibility_periods_used, " # "gaia.astrometric_excess_noise, gaia.astrometric_sigma5d_max,gaia.astrometric_excess_noise_sig, gaia.phot_bp_rp_excess_factor, " # "gaia.phot_g_mean_mag, 1.0857*gaia.phot_g_mean_flux_error/gaia.phot_g_mean_flux as phot_g_mean_mag_err, gaia.phot_bp_mean_mag, " # "1.0857*gaia.phot_bp_mean_flux_error/gaia.phot_bp_mean_flux as phot_bp_mean_mag_err, gaia.phot_rp_mean_mag, 1.0857*gaia.phot_rp_mean_flux_error/gaia.phot_rp_mean_flux as phot_rp_mean_mag_err, " # "unwise_obj_primary.unwise_objid, 22.5-2.5*log10(clip(unwise_obj_primary.flux(0), 1e-30, inf)) as W1, " # "22.5-2.5*log10(clip(unwise_obj_primary.flux(1), 1e-30, inf)) as W2, " # "1.0857*(unwise_obj_primary.dflux(0)/clip(unwise_obj_primary.flux(0),1e-12,inf)) as W1_err, " # "1.0857*(unwise_obj_primary.dflux(1)/clip(unwise_obj_primary.flux(1),1e-12,inf)) as W2_err, " # "unwise_obj_primary.fracflux, unwise_obj_primary.flags_unwise, unwise_obj_primary.primary, unwise_obj_primary.primary12, unwise_obj_primary.flags_info, " # "ukidss.j_1AperMag3, ukidss.j_1AperMag3Err, ukidss.hAperMag3, ukidss.hAperMag3Err, ukidss.kAperMag3, ukidss.kAperMag3Err, " # "ukidss.pStar, ukidss.j_1ppErrBits, ukidss.hppErrBits, ukidss.kppErrBits " # "from ucal_fluxqz, tmass(outer, matchedto=ucal_fluxqz, dmax=1.0, nmax=1), " # "gaia_dr2_source(matchedto=ucal_fluxqz, dmax=1, nmax=1) as gaia, " # "unwise_obj_primary(outer, matchedto=ucal_fluxqz, dmax=1, nmax=1) as unwise_obj_primary, " # "ukidss_las_dr10(outer, matchedto=ucal_fluxqz, dmax=1, nmax=1) as ukidss " # "where (numpy.sum(nmag_ok, axis=1) >= 4) & (numpy.sum(mean - mean_ap < 0.1, axis=1) >= 4) " # "& (numpy.abs(b) > 10) & ((gaia.parallax - 2*gaia.parallax_error) < 0.5) & (mean[:,1] >= 19.) & (mean[:,1] < 20) ") query = db.query(query) # Initialize stats on pixels, # of stars, etc. l_min = np.inf l_max = -np.inf b_min = np.inf b_max = -np.inf N_stars = 0 N_pix = 0 N_min = np.inf N_max = -np.inf N_det_hist = np.zeros(8, dtype='i8') N_per_band = np.zeros(8, dtype='i8') N_gaia = 0 N_gaia_hq = 0 N_in_pixel = [] N_pix_too_sparse = 0 N_pix_out_of_bounds = 0 fnameBase = abspath(values.out) fnameSuffix = 'h5' if fnameBase.endswith('.h5'): fnameBase = fnameBase[:-3] elif fnameBase.endswith('.hdf5'): fnameBase = fnameBase[:-5] fnameSuffix = 'hdf5' f = None nFiles = 0 nInFile = 0 # Write each pixel to the same file nest = True star_count=0 for (pix_info, obj) in query.execute([(mapper, values.nside_min, nest, values.bounds), reducer, (subdivider, values.nside_min, n_stars_max, values.min_stars, values.nside_max)], bounds=query_bounds, nworkers=values.n_workers): nside, pix_index = pix_info star_count=star_count+len(obj) # Filter out pixels that are outside of bounds l_center, b_center = pix2lb(nside, pix_index, nest=nest) if values.bounds is not None: if not hputils.lb_in_bounds(l_center, b_center, values.bounds): N_pix_out_of_bounds += 1 continue # Filter out pixels that have too few stars if len(obj) < values.min_stars: N_pix_too_sparse += 1 print '({:.3f}, {:.3f}) : {:d} stars'.format(l_center, b_center, len(obj)) continue # Prepare output for pixel EBV = np.percentile(obj['EBV'][:], 95.) # Open output file if f == None: fname = '%s.%.5d.%s' % (fnameBase, nFiles, fnameSuffix) f = h5py.File(fname, 'w') nInFile = 0 nFiles += 1 # Write to file gal_lb = to_file(f, pix_index, nside, nest, EBV, obj) # Update stats N_pix += 1 stars_in_pix = len(obj) N_stars += stars_in_pix nInFile += stars_in_pix N_in_pixel.append(stars_in_pix) n_pixels_at_res[nside] += 1 #N_per_band += np.sum(obj['err'] < 1.e9, axis=0) #N_gaia += int(np.sum(obj['pi_err'] < 1.e9)) #N_gaia_hq += int(np.sum(obj['pi'] / obj['pi_err'] > 5.)) #N_det_tmp = np.sum(obj['err'] < 1.e9, axis=1) #for k in xrange(8): # N_det_hist[k] += np.sum(N_det_tmp == k+1) if gal_lb[0] < l_min: l_min = gal_lb[0] if gal_lb[0] > l_max: l_max = gal_lb[0] if gal_lb[1] < b_min: b_min = gal_lb[1] if gal_lb[1] > b_max: b_max = gal_lb[1] if stars_in_pix < N_min: N_min = stars_in_pix if stars_in_pix > N_max: N_max = stars_in_pix # Close file if size exceeds max_stars if nInFile >= values.max_stars: f.close() f = None if f is not None: f.close() if N_pix != 0: N_in_pixel = np.array(N_in_pixel) print '# of stars in footprint: %d.' % N_stars print '# of pixels in footprint: %d.' % N_pix print 'Stars per pixel:' print ' min: %d' % N_min print ' 5%%: %d' % (np.percentile(N_in_pixel, 5.)) print ' 50%%: %d' % (np.percentile(N_in_pixel, 50.)) print ' mean: %d' % (float(N_stars) / float(N_pix)) print ' 95%%: %d' % (np.percentile(N_in_pixel, 95.)) print ' max: %d' % N_max #print '# of stars detected in each band:' #for i,b in enumerate('grizyJHK'): # print ' %s: %d (%.3f %%)' % (b, N_per_band[i], float(N_per_band[i])/float(N_stars)*100.) #print '# of stars with n-band detections:' #for i in xrange(8): # print ' %d: %d (%.3f %%)' % (i+1, N_det_hist[i], float(N_det_hist[i])/float(N_stars)*100.) #print '# of stars with Gaia parallaxes: {:d}'.format(N_gaia) #print '# of stars with Gaia parallax/err > 5: {:d}'.format(N_gaia_hq) print '# of pixels at each nside resolution:' for nside in nside_options: area_per_pix = hp.pixelfunc.nside2pixarea(nside, degrees=True) #area_per_pix = 4.*np.pi * (180./np.pi)**2. / (12. * nside**2.) area = n_pixels_at_res[nside] * area_per_pix print ' %d: %d (%.2f deg^2)' % (nside, n_pixels_at_res[nside], area) pct_sparse = 100. * float(N_pix_too_sparse) / float(N_pix_too_sparse + N_pix) print '# of pixels too sparse: %d (%.3f %%)' % (N_pix_too_sparse, pct_sparse) pct_out_of_bounds = 100. * float(N_pix_out_of_bounds) / float(N_pix_out_of_bounds + N_pix) print '# of pixels out of bounds: %d (%.3f %%)' % (N_pix_out_of_bounds, pct_out_of_bounds) print '# of files: %d.' % nFiles else: print 'No pixels in specified bounds with sufficient # of stars.' print '# of pixels out of bounds: %d' % (N_pix_out_of_bounds) if (values.bounds is not None) and (np.sum(N_pix) != 0): print '' print 'Bounds of included pixel centers:' print '\t(l_min, l_max) = (%.3f, %.3f)' % (l_min, l_max) print '\t(b_min, b_max) = (%.3f, %.3f)' % (b_min, b_max) print('star count',star_count) return 0 if __name__ == '__main__': main()