diff --git a/CEP/PyBDSM/doc/source/conf.py b/CEP/PyBDSM/doc/source/conf.py
index 46d6a3618f3c1ca0423a7cc76f15dd9e8557d8a9..8c90f2bf99ade4d355f6f9e36c200246387afa7e 100644
--- a/CEP/PyBDSM/doc/source/conf.py
+++ b/CEP/PyBDSM/doc/source/conf.py
@@ -48,9 +48,9 @@ copyright = u'2012, David Rafferty and Niruj Mohan'
# built documents.
#
# The short X.Y version.
-version = '1.2'
+version = '1.3'
# The full version, including alpha/beta/rc tags.
-release = '1.2'
+release = '1.3'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
diff --git a/CEP/PyBDSM/doc/source/whats_new.rst b/CEP/PyBDSM/doc/source/whats_new.rst
index a6094d39aa87175b9b7e256110366829ec39b2f7..b6ea18e04f2be6c0a0d4148feaeb5508adfa1eb7 100644
--- a/CEP/PyBDSM/doc/source/whats_new.rst
+++ b/CEP/PyBDSM/doc/source/whats_new.rst
@@ -4,6 +4,26 @@
What's New
**********
+Version 1.3 (2012/07/03):
+
+ * Fixed a bug in the calculation of positional errors for Gaussians.
+
+ * Adjusted ``rms_box`` algorithm to check for negative rms values (due to interpolation with cubic spline). If negative values are found, either the box size is increased or the interpolation is done with ``order=1`` (bilinear) instead.
+
+ * Output now includes the residual image produced using only wavelet Gaussians (if any) when ``atrous_do=True`` and ``output_all=True``.
+
+ * Improved organization of files when ``output_all=True``.
+
+ * Added logging of simple statistics (mean, std. dev, skew, and kurtosis) of the residual images.
+
+ * Included image rotation (if any) in beam definition. Rotation angle can vary across the image (defined by image WCS).
+
+ * Added Sagecal output format for Gaussian catalogs.
+
+ * Added check for newer versions of the PyBDSM software ``tar.gz`` file available on ftp.strw.leidenuniv.nl.
+
+ * Added total island flux (from sum of pixels) to ``gaul`` and ``srl`` catalogs.
+
Version 1.2 (2012/06/06):
* Added option to output flux densities for every channel found by the spectral index module.
@@ -12,7 +32,7 @@ Version 1.2 (2012/06/06):
* Implemented an adaptive scaling scheme for the ``rms_box`` parameter that shrinks the box size near bright sources and expands it far from them (enabled with the ``adaptive_rms_box`` option when ``rms_box`` is None). This scheme generally results in improved rms and mean maps when both strong artifacts and extended sources are present.
- * Improved speed of Gaussian fitting to wavelet images.
+ * Improved speed of Gaussian fitting to wavelet images.
* Added option to calculate fluxes within a specified aperture radius in pixels (set with the ``aperture`` option). Aperture fluxes, if measured, are output in the ``srl`` format catalogs.
diff --git a/CEP/PyBDSM/src/python/_version.py b/CEP/PyBDSM/src/python/_version.py
index f6438f95f06bba93e1af4fb2a43aa7dc206c325f..e08b867db4bd7470620e04fdfe89f4ba7e5feefa 100644
--- a/CEP/PyBDSM/src/python/_version.py
+++ b/CEP/PyBDSM/src/python/_version.py
@@ -9,7 +9,7 @@ adding to the changelog will naturally do this.
"""
# Version number
-__version__ = '1.2'
+__version__ = '1.3'
# Store svn Revision number. For this to work, one also needs to do:
#
@@ -27,6 +27,16 @@ def changelog():
PyBDSM Changelog.
-----------------------------------------------------------------------
+ 2012/07/03 - Version 1.3
+
+ 2012/07/03 - Fixed a bug in calculation of the positional errors of
+ Gaussians.
+
+ 2012/07/02 - Adjusted rms_box algorithm to check for negative rms
+ values (due to interpolation with cubic spline). If negative
+ values are found, either the box size is increased or the
+ interpolation is done with order=1 (bilinear) instead.
+
2012/06/28 - Output now includes the residual image produced by
using only wavelet Gaussians (if any) when atrous_do=True and
output_all=True. Improved organization of files when
@@ -34,6 +44,7 @@ def changelog():
std. dev, skew, and kurtosis) of the residual images.
2012/06/22 - Included image rotation (if any) in beam definition.
+ Rotation angle can vary across the image (defined by image WCS).
2012/06/19 - Changed exception handling to raise exceptions when
the interactive shell is not being used. Fixed bug that
diff --git a/CEP/PyBDSM/src/python/functions.py b/CEP/PyBDSM/src/python/functions.py
index 26e63dfdaf6d447a55cb092669520e5bda76b9a3..f1b668eb1773474bbe43d731fde1bbe071081adc 100755
--- a/CEP/PyBDSM/src/python/functions.py
+++ b/CEP/PyBDSM/src/python/functions.py
@@ -717,7 +717,7 @@ def get_errors(img, p, stdav, bm_pix=None):
pp = p[i*7:i*7+7]
### Now do error analysis as in Condon (and fBDSM)
size = pp[3:6]
- size = corrected_size(size)
+ size = corrected_size(size) # angle is now degrees CCW from +y-axis
sq2 = sqrt(2.0)
if bm_pix == None:
bm_pix = N.array([img.pixel_beam[0]*fwsig, img.pixel_beam[1]*fwsig, img.pixel_beam[2]])
@@ -729,10 +729,11 @@ def get_errors(img, p, stdav, bm_pix=None):
d2 = (size[0]*size[0]-size[1]*size[1])/(size[0]*size[0])
try:
e_peak = pp[0]*sq2/(dumrr3*pow(dumrr1,0.75)*pow(dumrr2,0.75))
- e_x0=size[0]/fwsig*sq2/(d1*dumrr3*pow(dumrr1,1.25)*pow(dumrr2,0.25))
- e_y0=size[1]/fwsig*sq2/(d1*dumrr3*pow(dumrr1,0.25)*pow(dumrr2,1.25))
e_maj=size[0]*sq2/(dumrr3*pow(dumrr1,1.25)*pow(dumrr2,0.25))
e_min=size[1]*sq2/(dumrr3*pow(dumrr1,0.25)*pow(dumrr2,1.25)) # in fw
+ pa_rad = size[2]*pi/180.0
+ e_x0 = sqrt( (e_maj*N.sin(pa_rad))**2 + (e_min*N.cos(pa_rad))**2 ) / d1
+ e_y0 = sqrt( (e_maj*N.cos(pa_rad))**2 + (e_min*N.sin(pa_rad))**2 ) / d1
e_pa=2.0/(d2*dumrr3*pow(dumrr1,0.25)*pow(dumrr2,1.25))
e_pa=e_pa*180.0/pi
e_tot=pp[0]*sqrt(e_peak*e_peak/(pp[0]*pp[0])+(0.25/dumr/dumr)*(e_maj*e_maj/(size[0]*size[0])+e_min*e_min/(size[1]*size[1])))
diff --git a/CEP/PyBDSM/src/python/gaul2srl.py b/CEP/PyBDSM/src/python/gaul2srl.py
index 96cbf91369240f38c2ec1211037ae5b08b41f307..438035dce9b18b81a7dfe22c862ca0674c73aacc 100644
--- a/CEP/PyBDSM/src/python/gaul2srl.py
+++ b/CEP/PyBDSM/src/python/gaul2srl.py
@@ -336,7 +336,7 @@ class Op_gaul2srl(Op):
gaus_c = [mompara[3], mompara[4], mompara[5]]
gaus_bm = [bm_pix[0], bm_pix[1], bm_pix[2]]
gaus_dc, err = func.deconv2(gaus_bm, gaus_c)
- deconv_size_sky = [img.pix2beam(gaus_dc), [0.0, 0.0, 0.0]]
+ deconv_size_sky = [img.pix2beam(gaus_dc, [mompara[1]+delc[0], mompara[2]+delc[1]]), [0.0, 0.0, 0.0]]
# update all objects etc
tot = 0.0
@@ -346,7 +346,7 @@ class Op_gaul2srl(Op):
totE_sq += g.total_fluxE**2
totE = sqrt(totE_sq)
size_pix = [mompara[3], mompara[4], mompara[5]]
- size_sky = img.pix2beam(size_pix)
+ size_sky = img.pix2beam(size_pix, [mompara[1]+delc[0], mompara[2]+delc[1]])
# Estimate errors using Monte Carlo technique
nMC = 20
@@ -376,10 +376,20 @@ class Op_gaul2srl(Op):
mompara5_MC[i] = mompara[5]
mompara0E = N.std(mompara0_MC)
mompara1E = N.std(mompara1_MC)
+ if mompara1E > 2.0*mompara[1]:
+ mompara1E = 2.0*mompara[1] # Don't let errors get too large
mompara2E = N.std(mompara2_MC)
+ if mompara2E > 2.0*mompara[2]:
+ mompara2E = 2.0*mompara[2] # Don't let errors get too large
mompara3E = N.std(mompara3_MC)
+ if mompara3E > 2.0*mompara[3]:
+ mompara3E = 2.0*mompara[3] # Don't let errors get too large
mompara4E = N.std(mompara4_MC)
+ if mompara4E > 2.0*mompara[4]:
+ mompara4E = 2.0*mompara[4] # Don't let errors get too large
mompara5E = N.std(mompara5_MC)
+ if mompara5E > 2.0*mompara[5]:
+ mompara5E = 2.0*mompara[5] # Don't let errors get too large
size_skyE = [mompara3E*sqrt(cdeltsq), mompara4E*sqrt(cdeltsq), mompara5E]
sraE, sdecE = (mompara1E*sqrt(cdeltsq), mompara2E*sqrt(cdeltsq))
diff --git a/CEP/PyBDSM/src/python/gausfit.py b/CEP/PyBDSM/src/python/gausfit.py
index e21ca2b021b3cbcea9c3783fad6e46578b731aee..6ce6110d3675934ab97531df462ca6bf38b56165 100644
--- a/CEP/PyBDSM/src/python/gausfit.py
+++ b/CEP/PyBDSM/src/python/gausfit.py
@@ -862,7 +862,7 @@ class Gaussian(object):
size = img.pixel_beam
size = func.corrected_size(size) # gives fwhm and P.A.
self.size_pix = size # FWHM in pixels and P.A. CCW from +y axis
- self.size_sky = img.pix2beam(size) # FWHM in degrees and P.A. CCW from North
+ self.size_sky = img.pix2beam(size, self.centre_pix) # FWHM in degrees and P.A. CCW from North
# Check if this is a wavelet image. If so, use orig_pixel_beam
# for flux calculation, as pixel_beam has been altered to match
@@ -888,7 +888,7 @@ class Gaussian(object):
self.centre_pixE = errors[1:3]
self.centre_skyE = img.pix2coord(errors[1:3])
self.size_pixE = errors[3:6]
- self.size_skyE = img.pix2beam(errors[3:6])
+ self.size_skyE = img.pix2beam(errors[3:6], self.centre_pix)
self.rms = img.islands[isl_idx].rms
self.mean = img.islands[isl_idx].mean
self.total_flux_isl = img.islands[isl_idx].total_flux
@@ -902,7 +902,7 @@ class Gaussian(object):
if flag == 0:
#gaus_dc = func.deconv(bm_pix, size)
gaus_dc, err = func.deconv2(bm_pix, size)
- self.deconv_size_sky = img.pix2beam(gaus_dc)
+ self.deconv_size_sky = img.pix2beam(gaus_dc, self.centre_pix)
self.deconv_size_skyE = [0., 0., 0.]
else:
self.deconv_size_sky = [0., 0., 0.]
diff --git a/CEP/PyBDSM/src/python/psf_vary.py b/CEP/PyBDSM/src/python/psf_vary.py
index 8b622c12fb2db89d591309eb8691fb40cf45fc09..c772ab5a41b37b147f4a721b669ef701d125c901 100644
--- a/CEP/PyBDSM/src/python/psf_vary.py
+++ b/CEP/PyBDSM/src/python/psf_vary.py
@@ -268,12 +268,12 @@ class Op_psf_vary(Op):
gaus_c = img.beam2pix(src.size_sky)
gaus_bm = [psf_maj_int[src_pos_int]*fwsig, psf_min_int[src_pos_int]*fwsig, psf_pa_int[src_pos_int]*fwsig]
gaus_dc, err = func.deconv2(gaus_bm, gaus_c)
- src.deconv_size_sky = img.pix2beam(gaus_dc)
+ src.deconv_size_sky = img.pix2beam(gaus_dc, src_pos)
src.deconv_size_skyE = [0.0, 0.0, 0.0]
for g in src.gaussians:
gaus_c = img.beam2pix(g.size_sky)
gaus_dc, err = func.deconv2(gaus_bm, gaus_c)
- g.deconv_size_sky = img.pix2beam(gaus_dc)
+ g.deconv_size_sky = img.pix2beam(gaus_dc, g.centre_pix)
g.deconv_size_skyE = [0.0, 0.0, 0.0]
if img.opts.quiet == False:
bar2.spin()
diff --git a/CEP/PyBDSM/src/python/pybdsm.py b/CEP/PyBDSM/src/python/pybdsm.py
index b45352852482b5bef7ec6e8a23f8ff87cad5ad7c..60273a0da800f0eee8e9c47ca3424bacf4d94e2c 100644
--- a/CEP/PyBDSM/src/python/pybdsm.py
+++ b/CEP/PyBDSM/src/python/pybdsm.py
@@ -677,7 +677,17 @@ if aps_local_val == None:
file_list = []
file_list = f.nlst('pub/rafferty/PyBDSM')
f.close()
- if 'pub/rafferty/PyBDSM/PyBDSM-' + __version__ + '.tar.gz' not in file_list:
+ ftp_version = []
+ for file in file_list:
+ if 'tar.gz' in file:
+ ver_start_indx = file.find('-') + 1
+ ftp_version.append(float(file[ver_start_indx:ver_start_indx+3]))
+ if ftp_version == []:
+ # No files found, continue without message
+ pass
+ elif float(__version__) < max(ftp_version):
+ print __version__
+ print max(ftp_version)
print '\n' + '*' * 72
print "There appears to be a newer version of PyBDSM available at:"
print " ftp://ftp.strw.leidenuniv.nl/pub/rafferty/PyBDSM/"
diff --git a/CEP/PyBDSM/src/python/readimage.py b/CEP/PyBDSM/src/python/readimage.py
index c9c4ce578ba95b7e8bc42fee44c7932df21bc6f0..7585268146ea448336d357c76634411d4ce306e0 100644
--- a/CEP/PyBDSM/src/python/readimage.py
+++ b/CEP/PyBDSM/src/python/readimage.py
@@ -262,17 +262,19 @@ class Op_readimage(Op):
hdr = img.header
cdelt1, cdelt2 = img.wcs_obj.acdelt[0:2]
- brot = self.get_rot(img.wcs_obj) # beam rotation delta CCW (in degrees) between N and +y axis of image
### define beam conversion routines:
- def beam2pix(x):
+ def beam2pix(x, location=None):
""" Converts beam in deg to pixels.
+ location specifies the location in pixels (x, y) for which beam is desired
Input beam angle should be degrees CCW from North.
The output beam angle is degrees CCW from the +y axis of the image.
"""
bmaj, bmin, bpa = x
- s1 = abs(bmaj/cdelt1)
+ brot = self.get_rot(img, location) # beam rotation delta CCW (in degrees) between N and +y axis of image
+
+ s1 = abs(bmaj/cdelt1)
s2 = abs(bmin/cdelt2)
th = bpa + brot
return (s1, s2, th)
@@ -283,15 +285,17 @@ class Op_readimage(Op):
s2 = abs(y*cdelt2)
return (s1, s2)
- def pix2beam(x):
+ def pix2beam(x, location=None):
""" Converts beam in pixels to deg.
+ location specifies the location in pixels (x, y) for which beam is desired
Input beam angle should be degrees CCW from the +y axis of the image.
The output beam angle is degrees CCW from North.
"""
s1, s2, th = x
bmaj = abs(s1*cdelt1)
bmin = abs(s2*cdelt2)
+ brot = self.get_rot(img, location) # beam rotation delta CCW (in degrees) between N and +y axis of image
bpa = th - brot
if bmaj < bmin:
bmaj, bmin = bmin, bmaj
@@ -348,7 +352,7 @@ class Op_readimage(Op):
found = True
if not found: raise RuntimeError("No beam information found in image header.")
- ### convert beam into pixels
+ ### convert beam into pixels (at image center)
pbeam = beam2pix(beam)
pbeam = (pbeam[0]/fwsig, pbeam[1]/fwsig, pbeam[2]) # IN SIGMA UNITS
@@ -460,9 +464,22 @@ class Op_readimage(Op):
if sys[:3] == 'FK4': year = 1950.0
return year, code
- def get_rot(self, wcs_obj):
- """Returns CCW rotation angle (in degrees) between N and +y axis of image"""
- if wcs_obj.crota == []:
- return 0.0
+ def get_rot(self, img, location=None):
+ """Returns CCW rotation angle (in degrees) between N and +y axis of image
+
+ location specifies the location in pixels (x, y) for which beam is desired
+ """
+ if location == None:
+ x1 = int(img.image.shape[2]/2.0)
+ y1 = int(img.image.shape[3]/2.0)
else:
- return N.mean(wcs_obj.crota) * 180.0 / N.pi
+ x1, y1 = location
+ delta_x = 0
+ delta_y = 10
+ w1 = img.pix2sky((x1, y1))
+ w2 = img.pix2sky((x1+delta_x, y1+delta_y))
+
+ rot_ang_rad = N.arctan2((w2[0] - w1[0]) , (w2[1] - w1[1]))
+ return rot_ang_rad*180.0/N.pi
+
+
\ No newline at end of file
diff --git a/CEP/PyBDSM/src/python/rmsimage.py b/CEP/PyBDSM/src/python/rmsimage.py
index 2b369bf6cef17008bec4f58a89df6444ebee98a4..60004381bf96a9c2b401a55323aeb0dbfb9e3da0 100644
--- a/CEP/PyBDSM/src/python/rmsimage.py
+++ b/CEP/PyBDSM/src/python/rmsimage.py
@@ -187,50 +187,11 @@ class Op_rmsimage(Op):
img.rms_box2 = (bsize2, bstep2)
else:
img.rms_box2 = opts.rms_box
- if (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']):
- if do_adapt:
- if opts.rms_box_bright is None:
- mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels (small scale)')
- else:
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels (small scale)')
- if opts.rms_box is None:
- mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
- '(' + str(img.rms_box2[0]) + ', ' +
- str(img.rms_box2[1]) + ') pixels (large scale)')
- else:
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box2[0]) + ', ' +
- str(img.rms_box2[1]) + ') pixels (large scale)')
- mylogger.userinfo(mylog, 'Number of sources using small scale', str(len(isl_pos)))
- else:
- mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels')
- else:
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels')
else:
img.rms_box = opts.rms_box_bright
img.rms_box2 = opts.rms_box
- if (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']):
- if do_adapt:
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels (small scale)')
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box2[0]) + ', ' +
- str(img.rms_box2[1]) + ') pixels (large scale)')
- mylogger.userinfo(mylog, 'Number of sources using small scale', str(len(isl_pos)))
- else:
- mylogger.userinfo(mylog, 'Using user-specified rms_box',
- '(' + str(img.rms_box[0]) + ', ' +
- str(img.rms_box[1]) + ') pixels')
-
+ self.output_rmsbox_size(img)
+
map_opts = (opts.kappa_clip, img.rms_box, opts.spline_rank)
for ipol, pol in enumerate(pols):
data = ch0_images[ipol]
@@ -244,6 +205,7 @@ class Op_rmsimage(Op):
## calculate rms/mean maps if needed
if ((opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const'])) and img.rms_box[0] > min(img.ch0.shape)/4.0:
# rms box is too large - just use constant rms and mean
+ self.output_rmsbox_size(img)
mylog.warning('Size of rms_box larger than 1/4 of image size')
mylogger.userinfo(mylog, 'Using constant background rms and mean')
img.use_rms_map = False
@@ -251,19 +213,56 @@ class Op_rmsimage(Op):
else:
if (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']):
if len(data.shape) == 2: ## 2d case
- self.map_2d(data, mean, rms, mask, *map_opts, do_adapt=do_adapt,
- bright_pt_coords=isl_pos, rms_box2=img.rms_box2)
+ rms_ok = False
+ while not rms_ok:
+ self.map_2d(data, mean, rms, mask, *map_opts, do_adapt=do_adapt,
+ bright_pt_coords=isl_pos, rms_box2=img.rms_box2,
+ logname="PyBDSM."+img.log)
+ if N.any(rms < 0.0):
+ rms_ok = False
+ if (opts.rms_box_bright is None and do_adapt) or (opts.rms_box is None and not do_adapt):
+ # Increase box by 20%
+ img.rms_box = (int(img.rms_box[0]*1.2),int(img.rms_box[1]*1.2))
+ map_opts = (opts.kappa_clip, img.rms_box, opts.spline_rank)
+ else:
+ # User has specified box size, use order=1 to prevent negatives
+ if opts.spline_rank > 1:
+ mylog.warning('Negative values found in rms map interpolated with spline_rank = %i' % opts.spline_rank)
+ mylog.warning('Using spline_rank = 1 (bilinear interpolation) instead')
+ map_opts = (opts.kappa_clip, img.rms_box, 1)
+ else:
+ rms_ok = True
+
elif len(data.shape) == 3: ## 3d case
if not isinstance(mask, N.ndarray):
mask = N.zeros(data.shape[0], dtype=bool)
for i in range(data.shape[0]):
## iterate each plane
- self.map_2d(data[i], mean[i], rms[i], mask[i], *map_opts,
- do_adapt=do_adapt, bright_pt_coords=isl_pos,
- rms_box2=img.rms_box2)
+ rms_ok = False
+ while not rms_ok:
+ self.map_2d(data[i], mean[i], rms[i], mask[i], *map_opts,
+ do_adapt=do_adapt, bright_pt_coords=isl_pos,
+ rms_box2=img.rms_box2, logname="PyBDSM."+img.log)
+ if N.any(rms[i] < 0.0):
+ rms_ok = False
+ if (opts.rms_box_bright is None and do_adapt) or (opts.rms_box is None and not do_adapt):
+ # Increase box by 20%
+ img.rms_box = (int(img.rms_box[0]*1.2),int(img.rms_box[1]*1.2))
+ map_opts = (opts.kappa_clip, img.rms_box, opts.spline_rank)
+ else:
+ # User has specified box size, use order=1 to prevent negatives
+ if opts.spline_rank > 1:
+ mylog.warning('Negative values found in rms map interpolated with spline_rank = %i' % opts.spline_rank)
+ mylog.warning('Using spline_rank = 1 (bilinear interpolation) instead')
+ map_opts = (opts.kappa_clip, img.rms_box, 1)
+ else:
+ rms_ok = True
else:
mylog.critical('Image shape not handleable' + pol_txt)
raise RuntimeError("Can't handle array of this shape" + pol_txt)
+ self.output_rmsbox_size(img)
+ if do_adapt:
+ mylogger.userinfo(mylog, 'Number of sources using small scale', str(len(isl_pos)))
mylog.info('Background rms and mean images computed' + pol_txt)
## check if variation of rms/mean maps is significant enough:
@@ -425,7 +424,7 @@ class Op_rmsimage(Op):
def map_2d(self, arr, out_mean, out_rms, mask=False,
kappa=3, box=None, interp=1, do_adapt=False,
- bright_pt_coords=None, rms_box2=None):
+ bright_pt_coords=None, rms_box2=None, logname=''):
"""Calculate mean&rms maps and store them into provided arrays
Parameters:
@@ -434,6 +433,7 @@ class Op_rmsimage(Op):
mask: mask
kappa: clipping value for rms/mean calculations
box: tuple of (box_size, box_step) for calculating map
+ rms_box2 = large-scale box size
interp: order of interpolating spline used to interpolate
calculated map
do_adapt: use adaptive binning
@@ -499,6 +499,9 @@ class Op_rmsimage(Op):
mean_map = mean_map1
# Interpolate to image coords
+ # Check whether default (cubic) or user-specified order causes problems.
+ # If so, use order=1.
+ mylog = mylogger.logging.getLogger(logname+"Rmsimage")
nd.map_coordinates(rms_map, ax[-1::-1], order=interp, output=out_rms)
nd.map_coordinates(mean_map, ax[-1::-1], order=interp, output=out_mean)
@@ -842,3 +845,31 @@ class Op_rmsimage(Op):
src_center = [xindx, yindx]
return [slice(xlow, xhigh, None), slice(ylow, yhigh, None)], src_center
+
+ def output_rmsbox_size(self, img):
+ """Prints rms/mean box size"""
+ opts = img.opts
+ do_adapt = opts.adaptive_rms_box
+ mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"RMSimage")
+ if (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']):
+ if do_adapt:
+ if opts.rms_box_bright is None:
+ mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
+ '(' + str(img.rms_box[0]) + ', ' +
+ str(img.rms_box[1]) + ') pixels (small scale)')
+ else:
+ mylogger.userinfo(mylog, 'Using user-specified rms_box',
+ '(' + str(img.rms_box[0]) + ', ' +
+ str(img.rms_box[1]) + ') pixels (small scale)')
+ if opts.rms_box is None:
+ mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
+ '(' + str(img.rms_box2[0]) + ', ' +
+ str(img.rms_box2[1]) + ') pixels (large scale)')
+ else:
+ mylogger.userinfo(mylog, 'Using user-specified rms_box',
+ '(' + str(img.rms_box2[0]) + ', ' +
+ str(img.rms_box2[1]) + ') pixels (large scale)')
+ else:
+ mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)',
+ '(' + str(img.rms_box[0]) + ', ' +
+ str(img.rms_box[1]) + ') pixels')