Python MCUtils.print_inline方法代码示例

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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)

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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)

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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)

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
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['EQUINOX'],hdu.header['EPOCH'] = 2000., 2000.
	hdu.header['BAND'] = 1 if band=='NUV' else 2
	hdu.header['VERSION'] = 'v{v}'.format(v=__version__)

	return hdu

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
def countmap(band,skypos,tranges,skyrange,width=False,height=False,
			 verbose=0,memlight=False,hdu=False,retries=20):
	"""Create a count (cnt) map."""
	imsz = gxt.deg2pix(skypos,skyrange)
	count = np.zeros(imsz)
	for trange in tranges:
		# If memlight is requested, break the integration into
		#  smaller chunks.
		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))
			events = gQuery.getArray(gQuery.rect(band,skypos[0],skypos[1],t0,t1,
												 skyrange[0],skyrange[1]),
									 verbose=verbose,retries=retries)

			# Check that there is actually data here.
			if not events:
				if verbose>1:
					print "No data in "+str([t0,t1])
				continue

			times = np.array(events,dtype='float64')[:,0 ]/tscale
			coo   =	np.array(events,dtype='float64')[:,1:]

			# If there's no data, return a blank image.
			if len(coo)==0:
				if verbose:
					print 'No data in this frame: '+str([t0,t1])
				continue

			# Define World Coordinate System (WCS)
			wcs = define_wcs(skypos,skyrange,width=False,height=False)

			# Map the sky coordinates onto the focal plane
			foc = wcs.sip_pix2foc(wcs.wcs_world2pix(coo,1),1)

			# Bin the events into actual image pixels
			H,xedges,yedges=np.histogram2d(foc[:,1]-0.5,foc[:,0]-0.5,
								bins=imsz,range=([ [0,imsz[0]],[0,imsz[1]] ]))
			count += H

	return count

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [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 print_inline [as 别名]
def shift_and_add(filename,ix=None,verbose=1):
    """Performs a _shift and add_ style coadd on a sequence of images of a
    planet. Recenters the center of brightness to the center of the image.
    The file indicated by filename is assumed to be .ser format. The ix keyword
    takes a list of frame numbers to use.
    """
    if verbose>1:
        print 'Reading: {f}'.format(f=filename)
    header = ser.readheader(filename)
    cnt = header['FrameCount']
    img = np.zeros([header['ImageHeight'],header['ImageWidth']])
    for i in (ix if ix else range(cnt-1)):
        frame = ser.readframe(filename,i,header=header)
        w,h=find_cob(frame)
        vec = [header['ImageHeight']/2.-h,header['ImageWidth']/2.-w]
        if verbose:
            mc.print_inline('Shift and adding {i} by {vec}'.format(i=i,vec=vec))
        img+=scipy.ndimage.interpolation.shift(frame,vec)/cnt
    return img

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
def query_photons(band,ra0,dec0,tranges,radius,verbose=0):
    """Retrieve photons within an aperture from the database."""
    stream = []
    if verbose:
        print "Retrieving photons within {rad} degrees of [{r}, {d}]".format(
                                                        rad=radius,r=ra0,d=dec0)
    for trange in tranges:
        if verbose:
            mc.print_inline(" and between "+str(trange[0])+" and "+
                            str(trange[1])+".")
        thisstream = gQuery.getArray(
            gQuery.allphotons(band, ra0, dec0, trange[0], trange[1], radius),
                              verbose=verbose,retries=100)
        stream.extend(thisstream)
    stream = np.array(stream, 'f8').T
    colnames = ['t', 'ra', 'dec', 'xi', 'eta', 'x', 'y']
    dtypes = ['f8', 'f8', 'f8', 'f4', 'f4', 'f4', 'f4']
    cols = map(np.asarray, stream, dtypes)
    events = dict(zip(colnames, cols))
    events['t']/=tscale # Adjust the timestamp by tscale
    return events

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import print_inline [as 别名]
def movie_tbl(band,tranges,verbose=0,framesz=0,retries=20):
	"""Initialize a FITS table to contain movie frame information."""
	if verbose:
		mc.print_inline('Populating exposure time table.')
	tstarts,tstops,exptimes=[],[],[]
	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)):
			t1 = trange[1] if i==steps else t0+stepsz
			tstarts.append(t0)
			tstops.append(t1)
			exptimes.append(dbt.compute_exptime(band,[t0,t1],
							verbose=verbose,retries=retries))
	col1 = pyfits.Column(name='tstart',format='E',array=np.array(tstarts))
	col2 = pyfits.Column(name='tstop',format='E',array=np.array(tstops))
	col3 = pyfits.Column(name='exptime',format='E',array=np.array(exptimes))
	cols = pyfits.ColDefs([col1,col2,col3])
	tbl  = pyfits.BinTableHDU.from_columns(cols)

	return tbl