Python MCUtils类代码示例

def movie(band,skypos,tranges,skyrange,framesz=0,width=False,height=False,
	verbose=0,memlight=False,coadd=False,response=False,hdu=False,retries=20,
	detsize=1.1):
	"""Generate a movie (mov)."""
	# Not defining stepsz creates a single full depth image.
	print tranges
	if coadd or (len(tranges)==1 and not framesz) or (not len(tranges)):
		if verbose>2:
			print 'Coadding across '+str(tranges)
		mv = integrate_map(band,skypos,tranges,skyrange,width=width,
			height=height,verbose=verbose,memlight=memlight,hdu=hdu,
			retries=retries,response=response,detsize=detsize)
	else:
		for trange in tranges:
			stepsz = framesz if framesz else trange[1]-trange[0]
			steps = np.ceil((trange[1]-trange[0])/stepsz)
			for i,t0 in enumerate(np.arange(trange[0],trange[1],stepsz)):
				if verbose>1:
					mc.print_inline('Movie frame '+str(i+1)+' of '+
																str(int(steps)))
				t1 = trange[1] if i==steps else t0+stepsz

				img = integrate_map(band,skypos,[[t0,t1]],skyrange,
					width=width,height=height,verbose=verbose,
					memlight=memlight,hdu=hdu,retries=retries,
					response=response,detsize=detsize)
				if img.min() == 0 and img.max() == 0:
					if verbose>1:
						print 'No data in frame {i}. Skipping...'.format(i=i)
					continue
				try:
					mv.append(img)
				except:
					mv = [img]
	return np.array(mv)

def pullphotons(band, ra0, dec0, tranges, radius, events={}, verbose=0,
                tscale=1000., calpath='../cal/', chunksz=10e6):
    """Retrieve photons within an aperture from the database."""
    events = {'t':[],'ra':[],'dec':[],'xi':[],'eta':[]}
    if verbose:
        print "Retrieving photons at ["+str(ra0)+", "+str(dec0)+"] within a radius of "+str(radius)
    for trange in tranges:
        if verbose:
            mc.print_inline(" and between "+str(trange[0])+" and "+
                            str(trange[1])+".")
        stream = gQuery.getArray(
            gQuery.allphotons(band, ra0, dec0, trange[0], trange[1], radius),
                              verbose=verbose,retries=100)
        if not stream:
            continue
        events['t'] = events['t']+np.array(np.array(stream,
                                        dtype='float64')[:,0]/tscale).tolist()
        # The float64 precision _is_ significant for RA / Dec.
        events['ra'] = events['ra']+np.array(np.array(stream,
                                        dtype='float64')[:,1]).tolist()
        events['dec'] = events['dec']+np.array(np.array(stream,
                                        dtype='float64')[:,2]).tolist()
        events['xi'] = events['xi']+np.array(np.array(stream,
                                        dtype='float32')[:,3]).tolist()
        events['eta'] = events['eta']+np.array(np.array(stream,
                                        dtype='float32')[:,4]).tolist()
    events['t'] = np.array(events['t'],dtype='float64')
    events['ra'] = np.array(events['ra'],dtype='float64')
    events['dec'] = np.array(events['dec'],dtype='float64')
    events['xi'] = np.array(events['xi'],dtype='float32')
    events['eta'] = np.array(events['eta'],dtype='float32')
    events = hashresponse(band, events, calpath=calpath, verbose=verbose)
    return events

def integrate_map(band,skypos,tranges,skyrange,width=False,height=False,
				  verbose=0,memlight=False,hdu=False,retries=20,response=False,
				  detsize=1.1):
	""" Integrate an image over some number of time ranges. Use a reduced
	memory optimization (at the expense of more web queries) if requested.
	"""
	imsz = gxt.deg2pix(skypos,skyrange)
	img = np.zeros(imsz)
	for trange in tranges:
		# If memlight is requested, break the integration into
		#  smaller chunks.
		# FIXME: memlight gives slightly wrong answers right now
		# This is probably due to a quirk of SQL, per issue #140.
		# Deprecating memlight until this can be resolved.
		step = memlight if memlight else trange[1]-trange[0]
		for i in np.arange(trange[0],trange[1],step):
			t0,t1=i,i+step
			if verbose:
				mc.print_inline('Coadding '+str(t0)+' to '+str(t1))
			img += makemap(band,skypos,[t0,t1],skyrange,response=response,
							 verbose=verbose,detsize=detsize)
		if response: # This is an intensity map.
			img /= dbt.compute_exptime(band,trange,skypos=skypos,
	                             						verbose=verbose)

	return img

def fits_header(band,skypos,tranges,skyrange,width=False,height=False,
				verbose=0,hdu=False,retries=20):
	"""Populate a FITS header."""
	if verbose:
		mc.print_inline('Populating FITS header.')
	hdu = hdu if hdu else pyfits.PrimaryHDU()
	wcs = define_wcs(skypos,skyrange,width=width,height=height)
	hdu.header['CDELT1'],hdu.header['CDELT2'] = wcs.wcs.cdelt
	hdu.header['CTYPE1'],hdu.header['CTYPE2'] = wcs.wcs.ctype
	hdu.header['CRPIX1'],hdu.header['CRPIX2'] = wcs.wcs.crpix
	hdu.header['CRVAL1'],hdu.header['CRVAL2'] = wcs.wcs.crval
	#hdu.header['RA_CENT'],hdu.header['DEC_CENT'] = wcs.wcs.crval # Dupe.
	hdu.header['EQUINOX'],hdu.header['EPOCH'] = 2000., 2000.
	hdu.header['BAND'] = 1 if band=='NUV' else 2
	hdu.header['VERSION'] = 'v{v}'.format(v=__version__)
	# Do we want to set the following?
	#hdu.header['OW'] = 1
	#hdu.header['DIRECT'] = 1
	#hdu.header['GRISM'] = 0
	#hdu.header['OPAQUE'] = 0

	# Put the total exposure time into the primary header
	hdu.header['EXPTIME'] = 0.
	for trange in tranges:
		hdu.header['EXPTIME'] += dbt.compute_exptime(band,trange,
												verbose=verbose,retries=retries)

	if len(tranges)==1:
	# Put the time range into the primary header for a single frame image
		hdu.header['EXPSTART'],hdu.header['EXPEND'] = tranges[0]
		# These are the proper keywords for this:
		hdu.header['TIME-OBS'],hdu.header['TIME-END'] = tranges[0]

	return hdu

def movie(band,skypos,tranges,skyrange,framesz=0,width=False,height=False,
	verbose=0,memlight=False,coadd=False,response=False,hdu=False,retries=20):
	"""Generate a movie (mov)."""
	# Not defining stepsz creates a single full depth image.
	if coadd or (len(tranges)==1 and not framesz):
		if verbose>2:
			print 'Coadding across '+str(tranges)
		mv = integrate_map(band,skypos,tranges,skyrange,width=width,
			height=height,verbose=verbose,memlight=memlight,hdu=hdu,
			retries=retries,response=response)
		#rr.append(rrhr(band,skypos,tranges,skyrange,response=response,width=width,height=height,stepsz=1.,verbose=verbose,hdu=hdu,retries=retries)) if response else rr.append(np.ones(np.shape(mv)[1:]))
	else:
		for trange in tranges:
			stepsz = framesz if framesz else trange[1]-trange[0]
			steps = np.ceil((trange[1]-trange[0])/stepsz)
			for i,t0 in enumerate(np.arange(trange[0],trange[1],stepsz)):
				if verbose>1:
					mc.print_inline('Movie frame '+str(i+1)+' of '+
																str(int(steps)))
				t1 = trange[1] if i==steps else t0+stepsz

				img = integrate_map(band,skypos,[[t0,t1]],skyrange,
					width=width,height=height,verbose=verbose,
					memlight=memlight,hdu=hdu,retries=retries,
					response=response)
				try:
					mv.append(img)
				except:
					mv = [img]
					# # FIXME: This should not create an rr unless it's requested...
	# 			rr.append(rrhr(band,skypos,[[t0,t1]],skyrange,response=response,width=width,height=height,stepsz=1.,verbose=verbose,retries=retries)) if response else rr.append(np.ones(np.shape(mv)[1:]))

	return np.array(mv)

def movie(band,skypos,tranges,skyrange,framesz=0,width=False,height=False,
		  verbose=0,memlight=False,coadd=False,
		  response=False,hdu=False,retries=20):
	"""Generate a movie (mov)."""
	# Not defining stepsz effectively creates a count map.
	mv = []
	rr = []
	if coadd:
		if verbose>2:
			print 'Coadding across '+str(tranges)
		mv.append(countmap(band,skypos,tranges,skyrange,width=width,
				  height=height,verbose=verbose,memlight=memlight,
				  hdu=hdu,retries=retries))
		rr.append(rrhr(band,skypos,tranges,skyrange,response=response,width=width,height=height,stepsz=1.,verbose=verbose,hdu=hdu,retries=retries)) if response else rr.append(np.ones(np.shape(mv)[1:]))
	else:
		for trange in tranges:
			stepsz = framesz if framesz else trange[1]-trange[0]
			steps = np.ceil((trange[1]-trange[0])/stepsz)
			for i,t0 in enumerate(np.arange(trange[0],trange[1],stepsz)):
				if verbose>1:
					mc.print_inline('Movie frame '+str(i+1)+' of '+
																str(int(steps)))
				t1 = trange[1] if i==steps else t0+stepsz
				mv.append(countmap(band,skypos,[[t0,t1]],skyrange,width=width,
								   height=height,verbose=verbose,
								   memlight=memlight,hdu=hdu,retries=retries))
	# FIXME: This should not create an rr unless it's requested...
				rr.append(rrhr(band,skypos,[[t0,t1]],skyrange,response=response,width=width,height=height,stepsz=1.,verbose=verbose,retries=retries)) if response else rr.append(np.ones(np.shape(mv)[1:]))

	return np.array(mv),np.array(rr)

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)

def error(data,band,radius,annulus):
    N_a = 1
    N_b0 = (mc.area(annulus[1])-mc.area(annulus[0]))/mc.area(radius)
    N_b = data[band]['bg_eff_area']/mc.area(radius)
    B0 = data[band]['bg']
    B = data[band]['bg_cheese']
    S = gt.mag2counts(data[band]['mag'],band)*data[band]['t_eff']
    s2 = {'bg_cheese_err':(S-B)+(N_a+(N_a**2.)/N_b),
          'bg_err':(S-B0)+(N_a+(N_a**2.)/N_b0)}
    return s2

def define_wcs(skypos,skyrange,width=False,height=False,verbose=0,
			   pixsz=0.000416666666666667):
	"""Define the world coordinate system (WCS)."""
	if verbose:
		mc.print_inline('Defining World Coordinate System (WCS).')
	wcs = pywcs.WCS(naxis=2) # NAXIS = 2
	imsz = gxt.deg2pix(skypos,skyrange)
	wcs.wcs.cdelt = np.array([-pixsz,pixsz])
	wcs.wcs.ctype = ['RA---TAN','DEC--TAN']
	wcs.wcs.crpix = [(imsz[1]/2.)+0.5,(imsz[0]/2.)+0.5]
	wcs.wcs.crval = skypos
	return wcs

def rrhr(band,skypos,tranges,skyrange,width=False,height=False,stepsz=1.,
		 verbose=0,response=True,hdu=False,retries=20):
	"""Generate a high resolution relative response (rrhr) map."""
	imsz = gxt.deg2pix(skypos,skyrange)
	# TODO the if width / height

	flat, flatinfo = cal.flat(band)
	npixx,npixy 	= flat.shape
	fltsz 		= flat.shape
	pixsz = flatinfo['CDELT2']
	detsize = 1.25

	# Rotate the flat into the correct orientation to start.
	flat   = np.flipud(np.rot90(flat))

	# NOTE: This upsample interpolation is done _last_ in the canonical
	#	pipeline as part of the poissonbg.c routine.
	# 	The interpolation function is "congrid" in the same file.
	# TODO: Should this be first order interpolation? (i.e. bilinear)
	hrflat = scipy.ndimage.interpolation.zoom(flat,4.,order=0,prefilter=False)
	img = np.zeros(hrflat.shape)[
				hrflat.shape[0]/2.-imsz[0]/2.:hrflat.shape[0]/2.+imsz[0]/2.,
				hrflat.shape[1]/2.-imsz[1]/2.:hrflat.shape[1]/2+imsz[1]/2.]

	for trange in tranges:
		t0,t1=trange
		entries = gQuery.getArray(gQuery.aspect(t0,t1),retries=retries)
		n = len(entries)

		asptime = np.float64(np.array(entries)[:,2])/tscale
		aspra   = np.float32(np.array(entries)[:,3])
		aspdec  = np.float32(np.array(entries)[:,4])
		asptwist= np.float32(np.array(entries)[:,5])
		aspflags= np.float32(np.array(entries)[:,6])
		asptwist= np.float32(np.array(entries)[:,9])
		aspra0  = np.zeros(n)+skypos[0]
		aspdec0 = np.zeros(n)+skypos[1]

		xi_vec, eta_vec = gnomonic.gnomfwd_simple(
							aspra,aspdec,aspra0,aspdec0,-asptwist,1.0/36000.,0.)

		col = 4.*( ((( xi_vec/36000.)/(detsize/2.)*(detsize/(fltsz[0]*pixsz)) + 1.)/2. * fltsz[0]) - (fltsz[0]/2.) )
		row = 4.*( (((eta_vec/36000.)/(detsize/2.)*(detsize/(fltsz[1]*pixsz)) + 1.)/2. * fltsz[1]) - (fltsz[1]/2.) )

		vectors = mc.rotvec(np.array([col,row]),-asptwist)

		for i in range(n):
			if verbose>1:
				mc.print_inline('Stamping '+str(asptime[i]))
				# FIXME: Clean this mess up a little just for clarity.
	        	img += scipy.ndimage.interpolation.shift(scipy.ndimage.interpolation.rotate(hrflat,-asptwist[i],reshape=False,order=0,prefilter=False),[vectors[1,i],vectors[0,i]],order=0,prefilter=False)[hrflat.shape[0]/2.-imsz[0]/2.:hrflat.shape[0]/2.+imsz[0]/2.,hrflat.shape[1]/2.-imsz[1]/2.:hrflat.shape[1]/2+imsz[1]/2.]*dbt.compute_exptime(band,[asptime[i],asptime[i]+1],verbose=verbose,retries=retries)*gxt.compute_flat_scale(asptime[i]+0.5,band,verbose=0)

	return img

def hashresponse(band,events,verbose=0):
    """Given detector xi, eta, return the response at each position."""
    # Hash out the response correction
    if verbose:
        mc.print_inline("Applying the response correction.")
    flat, _ = cal.flat(band)
    events['col'], events['row'] = xieta2colrow(
                                            events['xi'], events['eta'], band)
    events['flat'] = flat[np.array(events['col'], dtype='int16'),
                          np.array(events['row'], dtype='int16')]
    events['scale'] = gxt.compute_flat_scale(events['t'], band)
    # TODO: Separately do the binlinearly interpolated response
    events['response'] = (events['flat']*events['scale'])
    return events

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

def bg_mask_sources(band,ra0,dec0,ras,decs,responses,sources,maskradius=1.5):
    # At present, masks to 1.5 sigma where FWHM = 2.3548*sigma
    for i in range(len(sources['ra'])):
        ix = np.where(mc.angularSeparation(sources['ra'][i], sources['dec'][i],
                ras,decs)>=(maskradius/2.3548)*np.median(sources['fwhm'][i,:]))
        ras, decs, responses = ras[ix], decs[ix], responses[ix]
    return ras,decs,responses

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

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"],
    )