Python MCUtils.distance方法代码示例

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import distance [as 别名]
def makemap(band,skypos,trange,skyrange,response=False,verbose=0,detsize=1.1):
	imsz = gxt.deg2pix(skypos,skyrange)
	photons = np.array(gQuery.getArray(gQuery.skyrect(band,
		skypos[0],skypos[1],trange[0],trange[1],skyrange[0],skyrange[1]),
		verbose=verbose),dtype='float64')
	try:
		events = {'t':photons[:,0 ]/tscale,'ra':photons[:,1],'dec':photons[:,2],
			'xi':photons[:,3],'eta':photons[:,4],
			'x':photons[:,5], 'y':photons[:,6]}
	except IndexError:
		if verbose>2:
			print 'No events found at {s} +/- {r} in {t}.'.format(
				s=skypos,r=skyrange,t=trange)
		return np.zeros(imsz)
	# Trim the data on detsize
	col, row = ct.xieta2colrow(events['xi'],events['eta'],band)
	ix = np.where((1.25/800.)*mc.distance(col,row,400,400)<=detsize)
	n = len(ix[0])
	m = len(col)
	#print 'With detsize {d} using {n} of {m} data.'.format(d=detsize,n=n,m=m)
	if n == 0:
		return np.zeros(imsz)
	for k in events.keys():
		events[k] = events[k][ix]
	events = ct.hashresponse(band,events)
	wcs = define_wcs(skypos,skyrange,width=False,height=False)
	coo = zip(events['ra'],events['dec'])
	foc = wcs.sip_pix2foc(wcs.wcs_world2pix(coo,1),1)
	weights = 1./events['response'] if response else None
	H,xedges,yedges=np.histogram2d(foc[:,1]-0.5,foc[:,0]-0.5,bins=imsz,
		range=([ [0,imsz[0]],[0,imsz[1]] ]),weights=weights)
	return H

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import distance [as 别名]
def construct_row(i, band, objid, mcat, data):
    # Note: mcat['skybg'] is in counts per second per square arcseconds
    #       where as gPhoton is reporting cps over the aperture area.
    return (
        objid,
        data["t0"][0],
        data["t1"][0],
        mcat[band]["expt"][i],
        data["exptime"][0],
        mcat["ra"][i],
        mcat["dec"][i],
        data["racent"][0],
        data["deccent"][0],
        mcat[band][4]["mag"][i],
        mcat[band][4]["err"][i],
        data["mag_bgsub_cheese"][0],
        data["mag_bgsub"][0],
        data["mag"][0],
        mc.distance(data["detxs"], data["detys"], 400, 400)[0],
        data["responses"][0],
        mcat[band]["skybg"][i],
        data["bg"]["simple"][0],
        data["bg"]["cheese"][0],
        data["bg"]["eff_area"],
    )

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import distance [as 别名]
def sigmaclip_bg(data,radius,annulus,skypos,maxiter=10,sigmaclip=3.,
				 gausslim=50.,verbose=0,pixsz=0.000416666666666667):
	"""Produce an estimate of background counts within an aperture (radius)
	using a sigma clipping method for extracting the background from an
	annulus.
	This attempts to reproduce the calcuations of the backcalc() function in
	mapaps/poissonbg.c of the mission pipeline. (Probably written by Ted Wyder.)
	"""
	# FIXME: Does not apply response!

	# This cut is now handled by the ugly loop below, which barely dodges a
	# conceptula issue about fractional pixels...
	#ix = np.where((d>annulus[0]) & (d<annulus[1]))

	imsz=gxt.deg2pix(skypos,[annulus[1]*2,annulus[1]*2])
	wcs=define_wcs(skypos,[annulus[1]*2,annulus[1]*2])
	foc_ra,foc_dec=wcs.sip_pix2foc(wcs.wcs_world2pix(data['ra'],data['dec'],1),1)
	H,xedges,yedges=np.histogram2d(foc_ra-0.5,foc_dec-0.5,bins=imsz,
				       range=([ [0,imsz[0]],[0,imsz[1]]]))

	# Convert Gaussian sigma to a probability
	problim = 0.5*scipy.special.erfc(sigmaclip/np.sqrt(2.0))

	# Mask out non-annulus regions... there's probalby a more pythonic way
	bgimg=np.copy(H)
	for i in range(H.shape[0]):
		for j in range(H.shape[1]):
			# Add a little buffer to account for pixel widths?
			# FIXME? including everything within the annulus...
#			if (mc.distance(H.shape[0]/2.,H.shape[1]/2.,i,j)<annulus[0]/pixsz or
			if	mc.distance(H.shape[0]/2.,H.shape[1]/2.,i,j)>annulus[1]/pixsz:#):

				bgimg[i,j]=-1

	ix=np.where(bgimg>=0)
	m,s=bgimg[ix].mean(),bgimg[ix].std()
	d = 1.
	for i in range(maxiter):
		if d<=10e-5 or m<2:
			continue
		if m>=gausslim:
			# Mask anything outside of 3 sigma from the mean (of unmasked data)
			klim=m+sigmaclip*np.sqrt(m)#s
			klo=m-sigmaclip*np.sqrt(m)#s
			if verbose:
				print 'Gaussian cut: {klo} to {klim}'.format(klo=klo,klim=klim)
		else:
			klim = scipy.special.gammainccinv(m,problim)
			klo = -1 # None
			if verbose:
				print 'Poisson cut: {klo} to {klim}'.format(klo=klo,klim=klim)
		ix = np.where((bgimg>=klim) | (bgimg<=klo))
		bgimg[ix]=-1
		ix=np.where(bgimg>=0)
		d = np.abs((bgimg[ix].mean()-m)/m)# - 1)
		m,s=bgimg[ix].mean(),bgimg[ix].std()
	ix = np.where(bgimg>=0)
	return mc.area(radius)*bgimg[ix].mean()/mc.area(pixsz)

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import distance [as 别名]
def getcurve(band, ra0, dec0, radius, annulus=None, stepsz=None, lcurve={},
             trange=None, tranges=None, verbose=0, coadd=False, minexp=1.,
             maxgap=1., maskdepth=20, maskradius=1.5,
             photonfile=None, detsize=1.1):
    skyrange = [np.array(annulus).max().tolist() if annulus else radius,
                np.array(annulus).max().tolist() if annulus else radius,]
    if verbose:
        mc.print_inline("Getting exposure ranges.")
    if tranges is None:
        tranges = dbt.fGetTimeRanges(band, [ra0, dec0], trange=trange,
                maxgap=maxgap, minexp=minexp, verbose=verbose, detsize=detsize)
    elif not np.array(tranges).shape:
        print "No exposure time at this location: [{ra},{dec}]".format(
                                                            ra=ra0,dec=dec0)
    # FIXME: Everything goes to hell if no exposure time is available...
    # TODO: Add an ability to specify or exclude specific time ranges
    if verbose:
        mc.print_inline("Moving to photon level operations.")
    # FIXME: This error handling is hideous.
    try:
        lcurve = quickmag(band, ra0, dec0, tranges, radius, annulus=annulus,
                          stepsz=stepsz, verbose=verbose, coadd=coadd,
                          maskdepth=maskdepth,
                          maskradius=maskradius,photonfile=photonfile)
        lcurve['cps'] = lcurve['sources']/lcurve['exptime']
        lcurve['cps_bgsub'] = (lcurve['sources']-
                               lcurve['bg']['simple'])/lcurve['exptime']
        lcurve['cps_bgsub_cheese'] = (lcurve['sources']-
                               lcurve['bg']['cheese'])/lcurve['exptime']
        lcurve['mag'] = gxt.counts2mag(lcurve['cps'],band)
        lcurve['mag_bgsub'] = gxt.counts2mag(lcurve['cps_bgsub'],band)
        lcurve['mag_bgsub_cheese'] = gxt.counts2mag(
                                            lcurve['cps_bgsub_cheese'],band)
        lcurve['flux'] = gxt.counts2flux(lcurve['cps'],band)
        lcurve['flux_bgsub'] = gxt.counts2flux(lcurve['cps_bgsub'],band)
        lcurve['flux_bgsub_cheese'] = gxt.counts2flux(
                                            lcurve['cps_bgsub_cheese'],band)
        lcurve['detrad'] = mc.distance(lcurve['detxs'],lcurve['detys'],400,400)
    except ValueError:
        lcurve['cps']=[]
        lcurve['cps_bgsub']=[]
        lcurve['cps_bgsub_cheese']=[]
        lcurve['mag']=[]
        lcurve['mag_bgsub']=[]
        lcurve['mag_bgsub_cheese']=[]
        lcurve['flux']=[]
        lcurve['flux_bgsub']=[]
        lcurve['flux_bgsub_cheese']=[]
        lcurve['detrad']=[]
    if verbose:
        mc.print_inline("Done.")
        mc.print_inline("")
    return lcurve

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import distance [as 别名]
def construct_row(i,band,objid,mcat,data):
    # Note: mcat['skybg'] is in counts per second per square arcseconds
    #       where as gPhoton is reporting cps over the aperture area.
    return (objid, data['t0'][0], data['t1'][0],
            mcat[band]['expt'][i], data['exptime'][0],
            mcat['ra'][i], mcat['dec'][i],
            data['racent'][0], data['deccent'][0],
            mcat[band][4]['mag'][i], mcat[band][4]['err'][i],
            data['mag_bgsub_cheese'][0],
            data['mag_bgsub'][0], data['mag'][0],
            mc.distance(data['detxs'],data['detys'],400,400)[0],
            data['responses'][0], mcat[band]['skybg'][i],
            data['bg']['simple'][0], data['bg']['cheese'][0],
            data['bg']['eff_area'])