本文整理汇总了Python中admit.util.AdmitLogging.AdmitLogging.debug方法的典型用法代码示例。如果您正苦于以下问题:Python AdmitLogging.debug方法的具体用法?Python AdmitLogging.debug怎么用?Python AdmitLogging.debug使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类admit.util.AdmitLogging.AdmitLogging
的用法示例。
在下文中一共展示了AdmitLogging.debug方法的11个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: __init__
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def __init__(self, upper=True):
self.version = "27-apr-2016"
if have_ADMIT:
self.table = utils.admit_root() + "/etc/vlsr.tab"
self.cat = read_vlsr(self.table,upper)
logging.debug("VLSR: %s, found %d entries" % (self.table,len(self.cat)))
else:
logging.warning("VLSR: Warning, no ADMIT, empty catalogue")
self.cat = {}
示例2: peakstats
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def peakstats(image, freq, sigma, nsigma, minchan, maxgap, psample, peakfit = False):
""" Go through a cube and find peaks in the spectral dimension
It will gather a table of <peak>,<freq>,<sigma> which can be
optionally used for plotting
"""
if psample < 0: return
cutoff = nsigma * sigma
madata = casautil.getdata(image)
data = madata.data
shape = data.shape
logging.debug("peakstats: shape=%s cutoff=%g" % (str(shape),cutoff))
#print "DATA SHAPE:",shape
#print "cutoff=",cutoff
nx = shape[0]
ny = shape[1]
nz = shape[2]
chan = np.arange(nz)
# prepare the segment finder
# we now have an array data[nx,ny,nz]
sum = 0.0
pval = []
mval = []
wval = []
for x in range(0,nx,psample):
for y in range(0,ny,psample):
s0 = data[x,y,:]
spec = ma.masked_invalid(s0)
sum += spec.sum()
# using abs=True is a bit counter intuitive, but a patch to deal with the confusion in
# ADMITSegmentFinder w.r.t abs usage
asf = ADMITSegmentFinder(pmin=nsigma, minchan=minchan, maxgap=maxgap, freq=freq, spec=spec, abs=True)
#asf = ADMITSegmentFinder(pmin=nsigma, minchan=minchan, maxgap=maxgap, freq=freq, spec=spec, abs=False)
f = asf.line_segments(spec, nsigma*sigma)
for s in f:
if False:
for i in range(s[0],s[1]+1):
print "# ",x,y,i,spec[i]
## area preserving and peak are correlated, 18% difference
## fitgauss1Dm was about 5"
## with fitgauss1D was about 30", and still bad fits
par = utils.fitgauss1Dm(chan[s[0]:s[1]+1], spec[s[0]:s[1]+1], True) # peak from max
#par = utils.fitgauss1Dm(chan[s[0]:s[1]+1], spec[s[0]:s[1]+1], False) # peak from area preserving
if peakfit:
(par,cov) = utils.fitgauss1D (chan[s[0]:s[1]+1], spec[s[0]:s[1]+1],par)
#print "FIND: ",x,y,s,cutoff,0.0,0.0,par[0],par[1],par[2],s[1]-s[0]+1
pval.append(par[0])
mval.append(par[1])
wval.append(par[2])
#print "SUM:",sum
return (np.array(pval),np.array(mval),np.array(wval))
示例3: test_debug
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def test_debug(self):
msg = "unit_test_debug_message"
Alogging.debug(msg)
found = False
r = open(self.logfile, 'r')
for line in r.readlines():
if msg in line:
if(self.verbose):
print "\nFound message > ", line
found = True
r.close()
break
self.assertTrue(found)
示例4: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
#.........这里部分代码省略.........
caption = "Emission characteristics as a function of channel, as derived by CubeStats_AT "
caption += "(cyan: global rms,"
caption += " green: noise per channel,"
caption += " blue: peak value per channel,"
caption += " red: peak/noise per channel)."
self._summary["spectra"] = SummaryEntry([0, 0, str(specbox), 'Channel', imfile, thumbfile , caption, fin], "CubeStats_AT", self.id(True))
self._summary["chanrms"] = SummaryEntry([float(sigma0), fin], "CubeStats_AT", self.id(True))
# @todo Will imstat["max"][0] always be equal to s['datamax']? If not, why not?
if 'datamax' in s:
self._summary["dynrange"] = SummaryEntry([float(s['datamax']/sigma0), fin], "CubeStats_AT", self.id(True))
else:
self._summary["dynrange"] = SummaryEntry([float(imstat0["max"][0]/sigma0), fin], "CubeStats_AT", self.id(True))
self._summary["datamean"] = SummaryEntry([imstat0["mean"][0], fin], "CubeStats_AT", self.id(True))
title = bdp_name + "_1"
xlab = 'log(Peak,Noise,P/N)'
myplot.histogram([y1,y2,y3],title,bdp_name+"_1",xlab=xlab,thumbnail=True)
imfile = myplot.getFigure(figno=myplot.figno,relative=True)
thumbfile = myplot.getThumbnail(figno=myplot.figno,relative=True)
image1 = Image(images={bt.PNG:imfile},thumbnail=thumbfile,thumbnailtype=bt.PNG,description="CubeStats_1")
b2.addimage(image1,"im1")
# note that the 'y2' can have been clipped, which can throw off stats.robust()
# @todo should set a mask for those.
title = bdp_name + "_2"
xlab = 'log(Noise))'
n = len(y2)
ry2 = stats.robust(y2)
y2_mean = ry2.mean()
y2_std = ry2.std()
if n>9: logging.debug("NORMALTEST2: %s" % str(scipy.stats.normaltest(ry2)))
myplot.hisplot(y2,title,bdp_name+"_2",xlab=xlab,gauss=[y2_mean,y2_std],thumbnail=True)
title = bdp_name + "_3"
xlab = 'log(diff[Noise])'
n = len(y2)
# dy2 = y2[0:-2] - y2[1:-1]
dy2 = ma.masked_equal(y2[0:-2] - y2[1:-1],0.0).compressed()
rdy2 = stats.robust(dy2)
dy2_mean = rdy2.mean()
dy2_std = rdy2.std()
if n>9: logging.debug("NORMALTEST3: %s" % str(scipy.stats.normaltest(rdy2)))
myplot.hisplot(dy2,title,bdp_name+"_3",xlab=xlab,gauss=[dy2_mean,dy2_std],thumbnail=True)
title = bdp_name + "_4"
xlab = 'log(Signal/Noise))'
n = len(y3)
ry3 = stats.robust(y3)
y3_mean = ry3.mean()
y3_std = ry3.std()
if n>9: logging.debug("NORMALTEST4: %s" % str(scipy.stats.normaltest(ry3)))
myplot.hisplot(y3,title,bdp_name+"_4",xlab=xlab,gauss=[y3_mean,y3_std],thumbnail=True)
title = bdp_name + "_5"
xlab = 'log(diff[Signal/Noise)])'
n = len(y3)
dy3 = y3[0:-2] - y3[1:-1]
rdy3 = stats.robust(dy3)
dy3_mean = rdy3.mean()
dy3_std = rdy3.std()
if n>9: logging.debug("NORMALTEST5: %s" % str(scipy.stats.normaltest(rdy3)))
myplot.hisplot(dy3,title,bdp_name+"_5",xlab=xlab,gauss=[dy3_mean,dy3_std],thumbnail=True)
示例5: find
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def find(self):
""" Method that does the segment finding
Parameters
----------
None
Returns
-------
Tuple containing a list of the segments, the cutoff used, the
noise level, and a mean baseline.
"""
if self.abs:
self.spec = abs(self.spec)
temp = np.zeros(len(self.spec))
#self.see = ma.masked_array(temp, mask=self.spec.mask)
# parameters (some now from the function argument)
logging.debug("MIN/MAX " + str(self.spec.min()) + " / " +\
str(self.spec.max()))
n = len(self.spec) # data and freq assumed to be same size
if self.hanning:
h = np.hanning(5) / 2 # normalize the convolution array
h = np.array([0.25, 0.5, 0.25])
data2 = np.convolve(self.spec, h, 'same')
else:
data2 = self.spec
if len(data2) != len(self.freq):
raise Exception("ulines: data2 and freq not same array")
# do the work
dr = stats.robust(data2, self.f)
noise = dr.std()
logging.debug("ROBUST: (mean/median/noise) " + \
str(dr.mean()) + " / " + str(ma.median(dr)) + " / " + str(noise))
#print "\n\nD2",data2,"\n"
data3 = ma.masked_invalid(data2)
#print "\n\nD3\n",data3,"\n"
ddiff = data3[1:n] - data3[0:n-1]
logging.debug("DIFF: (mean/stdev) " + str(ddiff.mean()) +\
" / " + str(ddiff.std()))
#print "\n\n",ddiff,"\n",self.f,"\n"
ddr = stats.robust(ddiff, self.f)
logging.debug("RDIFF: (mean/median/stdev) " + \
str(ddr.mean()) + " / " + str(ma.median(ddr)) + " / " + \
str(ddr.std()))
#plt.show()
if self.bottom:
# first remind the classic
dmean1 = dr.mean()
dstd1 = dr.std()
logging.debug("CLASSIC MEAN/SIGMA: " + str(dmean1) + \
" / " + str(dstd1))
# see if we can find a better one?
# k should really depend on nchan, (like an nsigma), 2-3 should be ok for most.
k = 2.5
dmin = dr.min()
dmax = dr.min() + 2 * k * ddr.std() / 1.414214
logging.debug("DMIN/DMAX: " + str(dmin) + " / " + \
str(dmax))
dm = ma.masked_outside(dr, dmin, dmax)
dmean = max(0.0, dm.mean()) # ensure that the mean is positive or 0.0
dstd = dm.std()
if self.noise is not None:
cutoff = self.pmin * self.noise
elif self.nomean:
cutoff = self.pmin * dstd
else:
cutoff = dmean + self.pmin * dstd
logging.debug("BETTER MEAN/SIGMA: " + str(dmean) + \
" / " + str(dstd))
else:
# classic simple, but fails when robust didn't kill off (enough of) the signal
# sigma will be too high, cutoff too high and you could have no lines if there
# is one strong lines
dmean = dr.mean()
dstd = dr.std()
if self.noise is not None:
cutoff = self.pmin * self.noise
elif self.nomean:
cutoff = self.pmin * dstd
else:
cutoff = dmean + self.pmin * dstd
logging.debug("CLASSIC MEAN/SIGMA: " + str(dmean) + \
" / " + str(dstd))
logging.debug("SEGMENTS: f=%g pmin=%g maxgap=%d minchan=%d" % \
(self.f, self.pmin, self.maxgap, self.minchan))
#print "\nDATA\n\n",data2,"\n\n"
segments = self.line_segments(data2, cutoff)
#print "SEGMENTS",segments
nlines = len(segments)
logging.debug("Found %d segments above cutoff %f" % \
(nlines, cutoff))
segp = []
rmax = data2.max() + 0.1 # + 0.05*(data2.max()-data2.min())
segp.append([self.freq[0], self.freq[n - 1], cutoff, cutoff])
segp.append([self.freq[0], self.freq[n - 1], dmean, dmean])
for (l, s) in zip(range(nlines), segments):
ch0 = s[0]
#.........这里部分代码省略.........
示例6: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
#.........这里部分代码省略.........
if type(xpos)==int:
# for integers, boxes are allowed, even multiple
box = '%d,%d,%d,%d' % (xpos,ypos,xpos,ypos)
# convention for summary is (box)
cbox = '(%d,%d,%d,%d)' % (xpos,ypos,xpos,ypos)
# use extend here, not append, we want individual values in a list
sd.extend([xpos,ypos,cbox])
caption = "Average Spectrum at %s" % cbox
if False:
# this will fail on 3D cubes (see CAS-7648)
imval[i] = casa.imval(self.dir(fin),box=box)
else:
# work around that CAS-7648 bug
# another approach is the ia.getprofile(), see CubeStats, this will
# also integrate over regions, imval will not (!!!)
region = 'centerbox[[%dpix,%dpix],[1pix,1pix]]' % (xpos,ypos)
caption = "Average Spectrum at %s" % region
imval[i] = casa.imval(self.dir(fin),region=region)
elif type(xpos)==str:
# this is tricky, to stay under 1 pixel , or you get a 2x2 back.
region = 'centerbox[[%s,%s],[1pix,1pix]]' % (xpos,ypos)
caption = "Average Spectrum at %s" % region
sd.extend([xpos,ypos,region])
imval[i] = casa.imval(self.dir(fin),region=region)
else:
print "Data type: ",type(xpos)
raise Exception,"Data type for region not handled"
dt.tag("imval")
flux = imval[i]['data']
if len(flux.shape) > 1: # rare case if we step on a boundary between cells?
logging.warning("source %d has spectrum shape %s: averaging the spectra" % (i,repr(flux.shape)))
flux = np.average(flux,axis=0)
logging.debug('minmax: %f %f %d' % (flux.min(),flux.max(),len(flux)))
smax.append(flux.max())
if i==0: # for first point record few extra things
if len(imval[i]['coords'].shape) == 2: # normal case: 1 pixel
freqs = imval[i]['coords'].transpose()[2]/1e9 # convert to GHz @todo: input units ok?
elif len(imval[i]['coords'].shape) == 3: # rare case if > 1 point in imval()
freqs = imval[i]['coords'][0].transpose()[2]/1e9 # convert to GHz @todo: input units ok?
else:
logging.fatal("bad shape %s in freq return from imval - SHOULD NEVER HAPPEN" % imval[i]['coords'].shape)
chans = np.arange(len(freqs)) # channels 0..nchans-1
unit = imval[i]['unit']
restfreq = casa.imhead(self.dir(fin),mode="get",hdkey="restfreq")['value']/1e9 # in GHz
dt.tag("imhead")
vel = (1-freqs/restfreq)*utils.c # @todo : use a function (and what about relativistic?)
# construct the Table for CubeSpectrum_BDP
# @todo note data needs to be a tuple, later to be column_stack'd
labels = ["channel" ,"frequency" ,"flux" ]
units = ["number" ,"GHz" ,unit ]
data = (chans ,freqs ,flux )
if i==0:
# plane 0 : we are allowing a multiplane table, so the first plane is special
table = Table(columns=labels,units=units,data=np.column_stack(data),planes=["0"])
else:
# planes 1,2,3.... are stacked onto the previous one
table.addPlane(np.column_stack(data),"%d" % i)
# example plot , one per position for now
if use_vel:
x = vel
xlab = 'VLSR (km/s)'
else:
示例7: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def run(self):
""" Main program for OverlapIntegral
"""
dt = utils.Dtime("OverlapIntegral")
self._summary = {}
chans =self.getkey("chans")
cmap = self.getkey("cmap")
normalize = self.getkey("normalize")
doCross = True
doCross = False
myplot = APlot(pmode=self._plot_mode,ptype=self._plot_type,abspath=self.dir())
dt.tag("start")
n = len(self._bdp_in)
if n==0:
raise Exception,"Need at least 1 input Image_BDP "
logging.debug("Processing %d input maps" % n)
data = range(n) # array in which each element is placeholder for the data
mdata = range(n) # array to hold the max in each array
summarytable = admit.util.Table()
summarytable.columns = ["File name","Spectral Line ID"]
summarytable.description = "Images used in Overlap Integral"
for i in range(n):
bdpfile = self._bdp_in[i].getimagefile(bt.CASA)
if hasattr(self._bdp_in[i],"line"):
line = getattr(self._bdp_in[i],"line")
logging.info("Map %d: %s" % (i,line.uid))
lineid = line.uid
else:
lineid="no line"
data[i] = casautil.getdata(self.dir(bdpfile),chans)
mdata[i] = data[i].max()
logging.info("shape[%d] = %s with %d good data" % (i,data[i].shape,data[i].count()))
if i==0:
shape = data[i].shape
outfile = self.mkext("testOI","oi")
else:
if shape != data[i].shape:
raise Exception,"Shapes not the same, cannot overlap them"
# collect the file names and line identifications for the summary
summarytable.addRow([bdpfile,lineid])
logging.regression("OI: %s" % str(mdata))
if len(shape)>2 and shape[2] > 1:
raise Exception,"Cannot handle 3D cubes yet"
if doCross:
# debug: produce all cross-corr's of the N input maps (expensive!)
crossn(data, myplot)
dt.tag("crossn")
b1 = Image_BDP(outfile)
self.addoutput(b1)
b1.setkey("image", Image(images={bt.CASA:outfile}))
dt.tag("open")
useClone = True
# to create an output dataset, clone the first input, but using the chans=ch0~ch1
# e.g. using imsubimage(infile,outfile=,chans=
if len(chans) > 0:
# ia.regrid() doesn't have the chans=
taskinit.ia.open(self.dir(self._bdp_in[0].getimagefile(bt.CASA)))
taskinit.ia.regrid(outfile=self.dir(outfile))
taskinit.ia.close()
else:
# 2D for now
if not useClone:
logging.info("OVERLAP out=%s" % outfile)
taskinit.ia.fromimage(infile=self.dir(self._bdp_in[0].getimagefile(bt.CASA)),
outfile=self.dir(outfile), overwrite=True)
taskinit.ia.close()
dt.tag("fromimage")
if n==3:
# RGB
logging.info("RGB mode")
out = rgb1(data[0],data[1],data[2], normalize)
else:
# simple sum
out = data[0]
for i in range(1,n):
out = out + data[i]
if useClone:
casautil.putdata_raw(self.dir(outfile),out,clone=self.dir(self._bdp_in[0].getimagefile(bt.CASA)))
else:
taskinit.ia.open(self.dir(outfile))
s1 = taskinit.ia.shape()
s0 = [0,0,0,0]
r1 = taskinit.rg.box(blc=s0,trc=s1)
pixeldata = out.data
pixelmask = ~out.mask
taskinit.ia.putregion(pixels=pixeldata, pixelmask=pixelmask, region=r1)
taskinit.ia.close()
title = "OverlapIntegral"
#.........这里部分代码省略.........
示例8: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def run(self):
#
self._summary = {} # prepare to make a summary here
dt = utils.Dtime("Ingest") # timer for debugging
do_cbeam = True # enforce a common beam
#
pb = self.getkey('pb')
do_pb = len(pb) > 0
use_pb = self.getkey("usepb")
#
create_mask = self.getkey('mask') # create a new mask ?
box = self.getkey("box") # corners in Z, XY or XYZ
edge = self.getkey("edge") # number of edge channels to remove
restfreq = self.getkey("restfreq") # < 0 means not activated
# smooth= could become deprecated, and/or include a decimation option to make it useful
# again, Smooth_AT() does this also , at the cost of an extra cube to store
smooth = self.getkey("smooth") #
#
vlsr = self.getkey("vlsr") # see also LineID, where this could be given again
# first place a fits file in the admit project directory (symlink)
# this is a bit involved, depending on if an absolute or relative path was
# give to Ingest_AT(file=)
fitsfile = self.getkey('file')
if fitsfile[0] != os.sep:
fitsfile = os.path.abspath(os.getcwd() + os.sep + fitsfile)
logging.debug('FILE=%s' % fitsfile)
if fitsfile[0] != os.sep:
raise Exception,"Bad file=%s, expected absolute name",fitsfile
# now determine if it could have been a CASA (or MIRIAD) image already
# which we'll assume if it's a directory; this is natively supported by CASA
# but there are tools where if you pass it a FITS or MIRIAD
# MIRIAD is not recommended for serious work, especially big files, since there
# is a performance penalty due to tiling.
file_is_casa = casautil.iscasa(fitsfile)
loc = fitsfile.rfind(os.sep) # find the '/'
ffile0 = fitsfile[loc+1:] # basename.fits
basename = self.getkey('basename') # (new) basename allowed (allow no dots?)
if len(basename) == 0:
basename = ffile0[:ffile0.rfind('.')] # basename
logging.info("basename=%s" % basename)
target = self.dir(ffile0)
if not os.path.exists(target) :
cmd = 'ln -s "%s" "%s"' % (fitsfile, target)
logging.debug("CMD: %s" % cmd)
os.system(cmd)
readonly = False
if file_is_casa:
logging.debug("Assuming input %s is a CASA (or MIRIAD) image" % ffile0)
bdpfile = self.mkext(basename,"im")
if bdpfile == ffile0:
logging.warning("No selections allowed on CASA image, since no alias was given")
readonly = True
b1 = SpwCube_BDP(bdpfile)
self.addoutput(b1)
b1.setkey("image", Image(images={bt.CASA:bdpfile}))
# @todo b2 and PB?
else:
# construct the output name and construct the BDP based on the CASA image name
# this also takes care of the behind the scenes alias= substitution
bdpfile = self.mkext(basename,"im")
if bdpfile == basename:
raise Exception,"basename and bdpfile are the same, Ingest_AT needs a fix for this"
b1 = SpwCube_BDP(bdpfile)
self.addoutput(b1)
if do_pb:
print "doing the PB"
bdpfile2 = self.mkext(basename,"pb")
b2 = Image_BDP(bdpfile2)
self.addoutput(b2)
# @todo we should also set readonly=True if no box, no mask etc. and still an alias
# that way it will speed up and not make a copy of the image ?
# fni and fno are full (abspath) filenames, ready for CASA
# fni is the same as fitsfile
fni = self.dir(ffile0)
fno = self.dir(bdpfile)
if do_pb: fno2 = self.dir(bdpfile2)
dt.tag("start")
if file_is_casa:
taskinit.ia.open(fni)
else:
if do_pb and use_pb:
# @todo this needs a fix for the path for pb, only works if abs path is given
# impbcor(im.fits,pb.fits,out.im,overwrite=True,mode='m')
if False:
# this may seem like a nice shortcut, to have the fits->casa conversion be done
# internally in impbcor, but it's a terrible performance for big cubes. (tiling?)
# we keep this code here, perhaps at some future time (mpi?) this performs better
# @todo fno2
impbcor(fni,pb,fno,overwrite=True,mode='m')
dt.tag("impbcor-1")
#.........这里部分代码省略.........
示例9: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
#.........这里部分代码省略.........
# the cube apparently must be closed afterwards and
# then reopened if spatial smoothing is to be done.
if velres['value'] > 0:
# handle the different units allowed. CASA doesn't
# like lowercase for hz units...
if not hertz_input:
freq_res = str(velres['value']*1000.0/CC *rest_freq )+'Hz'
else:
freq_res = str(velres['value'])
# try to convert velres to km/s for debug purposes
velres['value'] = velres['value']/rest_freq*CC / 1000.0
if(velres['unit'] == 'khz'):
velres['value'] = velres['value']*1000.0
velres['unit'] = 'kHz'
elif(velres['unit']=='mhz'):
velres['value'] = velres['value']*1E6
velres['unit'] = 'MHz'
elif(velres['unit']=='ghz'):
velres['value'] = velres['value']*1E9
velres['unit'] = 'GHz'
freq_res = freq_res + velres['unit']
# NB: there is apparently a bug in CASA. only smoothing along the frequency
# axis does not work. sepconvolve gives a unit error (says axis unit is radian rather
# than Hz). MUST smooth in 2+ dimensions if you want this to work.
if(velres['value'] < vel_scale):
raise Exception,"Desired velocity resolution %g less than pixel scale %g" % (velres['value'],vel_scale)
image_tmp = self.dir('tmp.smooth')
im2=taskinit.ia.sepconvolve(outfile=image_tmp,axes=[0,1,2], types=["boxcar","boxcar","gauss"],\
widths=['1pix','1pix',freq_res], overwrite=True)
im2.done()
logging.debug("sepconvolve to %s" % image_out)
# for some reason, doing this in memory does not seem to work, so outfile must be specified.
logging.info("Smoothing cube to a velocity resolution of %s km/s" % str(velres['value']))
logging.info("Smoothing cube to a frequency resolution of %s" % freq_res)
taskinit.ia.close()
taskinit.ia.open(image_tmp)
dt.tag("sepconvolve")
else:
image_tmp = image_out
# now do the spatial smoothing
convolve_to_min_beam = True # default is to convolve to a min enclosing beam
if bmaj > 0 and bmin > 0:
# form qa objects out of these so that casa can understand
bmaj = taskinit.qa.quantity(bmaj,'arcsec')
bmin = taskinit.qa.quantity(bmin,'arcsec')
bpa = taskinit.qa.quantity(bpa,'deg')
target_res={}
target_res['major'] = bmaj
target_res['minor'] = bmin
target_res['positionangle'] = bpa
# throw an exception if cannot be convolved
try:
# for whatever reason, if you give convolve2d a beam parameter,
# it complains ...
im2=taskinit.ia.convolve2d(outfile=image_out,major = bmaj,\
minor = bmin, pa = bpa,\
示例10: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
#.........这里部分代码省略.........
#args2["mask"] = ""
args2["point"] = False
args2["width"] = 5
args2["negfind"] = False
# set-up for SourceList_BDP
slbdp = SourceList_BDP(slbase)
# connect to casa image and call casa ia.findsources tool
taskinit.ia.open(self.dir(infile))
# findsources() cannot deal with 'Jy/beam.km/s' ???
# so for the duration of findsources() we patch it
bunit = taskinit.ia.brightnessunit()
if bpatch and bunit != 'Jy/beam':
logging.warning("Temporarely patching your %s units to Jy/beam for ia.findsources()" % bunit)
taskinit.ia.setbrightnessunit('Jy/beam')
else:
bpatch = False
atab = taskinit.ia.findsources(**args2)
if bpatch:
taskinit.ia.setbrightnessunit(bunit)
taskargs = "nsigma=%4.1f sigma=%g region=%s robust=%s snmax=%5.1f" % (nsigma,sigma,str(region),str(robust),snmax)
dt.tag("findsources")
nsources = atab["nelements"]
xtab = []
ytab = []
logscale = False
sumflux = 0.0
if nsources > 0:
# @TODO: Why are Xpix, YPix not stored in the table?
# -> PJT: I left them out since they are connected to an image which may not be available here
# but we should store the frequency of the observation here for later bandmerging
logging.debug("%s" % str(atab['component0']['shape']))
logging.info("Right Ascen. Declination X(pix) Y(pix) Peak Flux Major Minor PA SNR")
funits = atab['component0']['flux']['unit']
if atab['component0']['shape'].has_key('majoraxis'):
sunits = atab['component0']['shape']['majoraxis']['unit']
aunits = atab['component0']['shape']['positionangle']['unit']
else:
sunits = "n/a"
aunits = "n/a"
punits = taskinit.ia.summary()['unit']
logging.info(" %s %s %s %s %s" % (punits,funits,sunits,sunits,aunits))
#
# @todo future improvement is to look at image coordinates and control output appropriately
#
if ds9:
# @todo variable name
regname = self.mkext(infile,'ds9.reg')
fp9 = open(self.dir(regname),"w!")
for i in range(nsources):
c = "component%d" % i
name = "%d" % (i+1)
r = atab[c]['shape']['direction']['m0']['value']
d = atab[c]['shape']['direction']['m1']['value']
pixel = taskinit.ia.topixel([r,d])
xpos = pixel['numeric'][0]
ypos = pixel['numeric'][1]
rd = taskinit.ia.toworld([xpos,ypos],'s')
ra = rd['string'][0][:12]
dec = rd['string'][1][:12]
flux = atab[c]['flux']['value'][0]
sumflux = sumflux + flux
if atab[c]['shape'].has_key('majoraxis'):
smajor = atab[c]['shape']['majoraxis']['value']
示例11: run
# 需要导入模块: from admit.util.AdmitLogging import AdmitLogging [as 别名]
# 或者: from admit.util.AdmitLogging.AdmitLogging import debug [as 别名]
def run(self):
""" The run method, calculates the moments and creates the BDP(s)
Parameters
----------
None
Returns
-------
None
"""
self._summary = {}
momentsummary = []
dt = utils.Dtime("Moment")
# variable to track if we are using a single cutoff for all moment maps
allsame = False
moments = self.getkey("moments")
numsigma = self.getkey("numsigma")
mom0clip = self.getkey("mom0clip")
# determine if there is only 1 cutoff or if there is a cutoff for each moment
if len(moments) != len(numsigma):
if len(numsigma) != 1:
raise Exception("Length of numsigma and moment lists do not match. They must be the same length or the length of the cutoff list must be 1.")
allsame = True
# default moment file extensions, this is information copied from casa.immoments()
momentFileExtensions = {-1: ".average",
0: ".integrated",
1: ".weighted_coord",
2: ".weighted_dispersion_coord",
3: ".median",
4: "",
5: ".standard_deviation",
6: ".rms",
7: ".abs_mean_dev",
8: ".maximum",
9: ".maximum_coord",
10: ".minimum",
11: ".minimum_coord",
}
logging.debug("MOMENT: %s %s %s" % (str(moments), str(numsigma), str(allsame)))
# get the input casa image from bdp[0]
# also get the channels the line actually covers (if any)
bdpin = self._bdp_in[0]
infile = bdpin.getimagefile(bt.CASA)
chans = self.getkey("chans")
# the basename of the moments, we will append _0, _1, etc.
basename = self.mkext(infile, "mom")
fluxname = self.mkext(infile, "flux")
# beamarea = nppb(self.dir(infile))
beamarea = 1.0 # until we have it from the MOM0 map
sigma0 = self.getkey("sigma")
sigma = sigma0
dt.tag("open")
# if no CubseStats BDP was given and no sigma was specified, find a
# noise level via casa.imstat()
if self._bdp_in[1] is None and sigma <= 0.0:
raise Exception("A sigma or a CubeStats_BDP must be input to calculate the cutoff")
elif self._bdp_in[1] is not None:
sigma = self._bdp_in[1].get("sigma")
# immoments is a bit peculiar. If you give one moment, it will use
# exactly the outfile you picked for multiple moments, it will pick
# extensions such as .integrated [0], .weighted_coord [1] etc.
# we loop over the moments and will use the numeric extension instead.
# Might be laborious loop for big input cubes
#
# arguments for immoments
args = {"imagename" : self.dir(infile),
"moments" : moments,
"outfile" : self.dir(basename)}
# set the channels if given
if chans != "":
args["chans"] = chans
# error check the mom0clip input
if mom0clip > 0.0 and not 0 in moments:
logging.warning("mom0clip given, but no moment0 map was requested. One will be generated anyway.")
# add moment0 to the list of computed moments, but it has to be first
moments.insert(0,0)
if not allsame:
numsigma.insert(0, 2.0*sigma)
if allsame:
# this is only executed now if len(moments) > 1 and len(cutoff)==1
args["excludepix"] = [-numsigma[0] * sigma, numsigma[0] * sigma]
casa.immoments(**args)
dt.tag("immoments-all")
else:
# this is execute if len(moments)==len(cutoff) , even when len=1
for i in range(len(moments)):
args["excludepix"] = [-numsigma[i] * sigma, numsigma[i] * sigma]
args["moments"] = moments[i]
args["outfile"] = self.dir(basename + momentFileExtensions[moments[i]])
casa.immoments(**args)
#.........这里部分代码省略.........