From e349bd4f0cfd698f0f9bb7d4055d62add061c2ac Mon Sep 17 00:00:00 2001
From: David Rafferty <rafferty@strw.leidenuniv.nl>
Date: Tue, 3 Jul 2012 08:29:34 +0000
Subject: [PATCH] Task #3403: Fixed bug in positional errors; added image
rotation to beam; updated documentation; added check for negative values in
rms image.
---
CEP/PyBDSM/doc/source/conf.py | 4 +-
CEP/PyBDSM/doc/source/whats_new.rst | 22 ++++-
CEP/PyBDSM/src/python/_version.py | 13 ++-
CEP/PyBDSM/src/python/functions.py | 7 +-
CEP/PyBDSM/src/python/gaul2srl.py | 14 +++-
CEP/PyBDSM/src/python/gausfit.py | 6 +-
CEP/PyBDSM/src/python/psf_vary.py | 4 +-
CEP/PyBDSM/src/python/pybdsm.py | 12 ++-
CEP/PyBDSM/src/python/readimage.py | 37 +++++---
CEP/PyBDSM/src/python/rmsimage.py | 125 +++++++++++++++++-----------
10 files changed, 172 insertions(+), 72 deletions(-)
diff --git a/CEP/PyBDSM/doc/source/conf.py b/CEP/PyBDSM/doc/source/conf.py
index 46d6a3618f3..8c90f2bf99a 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 a6094d39aa8..b6ea18e04f2 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 f6438f95f06..e08b867db4b 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 26e63dfdaf6..f1b668eb177 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 96cbf913692..438035dce9b 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 e21ca2b021b..6ce6110d367 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 8b622c12fb2..c772ab5a41b 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 b4535285248..60273a0da80 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 c9c4ce578ba..7585268146e 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 2b369bf6cef..60004381bf9 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')
--
GitLab