search

ratt-ru / pyxis

Python Extensions for Interferometry Scripting
3 stars 5 forks source link

pybdsm_search in lsm.py does not work in polarized mode #30

Open modhurita opened 9 years ago

modhurita commented 9 years ago

I am using a modified version of a section of the 3C147 pipeline to run pybdsm in polarized mode:

lsm.pybdsmsearch(image='img$i.fits',thresh_pix=THRESH_PIX[0],thresh_isl=THRESH_ISL[0],select="r.gt.30s",pol=True); lsm.tigger_convert("${lsm.PYBDSMOUTPUT} $DESTDIR/$LSMBASE$SUFFIX+pybdsm-apparent${i}.lsm.html -f");

Looking at the help for pybdsm_search (and looking at lsm.py), I gathered that all I needed to change to make pybdsm find polarization properties from an all-Stokes image was to set pol=True above.

Looking at the log file, it seems like it does compute polarization properties:

--> Opened 'img_0.fits' Image size .............................. : (512, 512) pixels Number of channels ...................... : 1 Number of Stokes parameters ............. : 4 Beam shape (major, minor, pos angle) .... : (0.0047, 0.00298, 80.8) degrees Frequency of image ...................... : 1424.500 MHz Number of blank pixels .................. : 0 (0.0%) Flux from sum of (non-blank) pixels ..... : 0.087 Jy --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (65, 22) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map (I) : (0.00005, 0.00028) Jy/beam Value of background mean (I) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (Q) : (0.00003, 0.00006) Jy/beam Value of background mean (Q) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (U) : (0.00003, 0.00006) Jy/beam Value of background mean (U) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (V) : (0.00003, 0.00005) Jy/beam Value of background mean (V) ............ : -0.0 Jy/beam --> Expected 5-sigma-clipped false detection rate < fdr_ratio --> Using sigma-clipping ('hard') thresholding Minimum number of pixels per island ..... : 10 Number of islands found ................. : 4 Fitting islands with Gaussians .......... : [=---] 1/4 Fitting islands with Gaussians .......... : [===] 3/4 Fitting islands with Gaussians .......... : [====] 4/4 Total number of Gaussians fit to image .. : 4 Total flux density in model ............. : 0.070 Jy --> Grouping Gaussians into sources Number of sources formed from Gaussians : 4

--> Checking PI image for new sources --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (51, 17) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map .... : (0.00002, 0.00005) Jy/beam Value of background mean ................ : 4e-05 Jy/beam Minimum number of pixels per island ..... : 10 Number of islands found ................. : 0 --> Grouping Gaussians into sources Number of sources formed from Gaussians : 0 New sources found in PI image ........... : 0 (4 total) Calculating polarisation properties .... : [=---] 1/4 Calculating polarisation properties .... : [===] 3/4 Calculating polarisation properties .... : [====] 4/4 2015/03/02 17:32:55 INFO: writing PyBDSM gaul catalog --> Wrote ASCII file '3C147-CD-LO-spw0-s3_pybdsm.lsm.html.gaul'

However, in the final LSM, Q,U,V are 0. These aren't 0 in the gaul file. It seems like Q,U,V are not getting written properly into the final LSM.

Location of the files, on elwood:

Gaul file: /home/mitra/PrimaryBeamReduction/3C147-CD-LO-spw0-s3_pybdsm.lsm.html.gaul Final LSM: /home/mitra/PrimaryBeamReduction/reduction-02Mar15/plots-3C147-CD-LO-spw0/3C147-GdB-spw0+pybdsm-apparent_0.lsm.html Log file: /home/mitra/PrimaryBeamReduction/reduction-02Mar15/log-3C147-CD-LO-spw0.txt

o-smirnov commented 9 years ago

Try it wish the fix I just pushed. Problem was a missing space on the command-line invocation, when adding polarization parameters...

modhurita commented 9 years ago

That doesn't fix it, and now even the gaul file doesn't have any polarization values - the last column is the S_Code one.

o-smirnov commented 9 years ago

My change only affects the format string passed to tigger-convert (in the old version, S_Code and q were stuck together, thus upsetting the order of the polarization fields). If the gaul file has no QUV, that's earlier on -- that means PyBDSM itself is invoked without the polarization flags. Are you still passing pol=True properly?

modhurita commented 9 years ago

Yes, the log says polarisation_do=True, but it isn't computing polarization properties any more:

2015/03/02 21:16:42 INFO: running PyBDSM process_image(img_0.fits,polarisation_do=True,quiet=True,thresh_isl=15,thresh_pix=50) --> Opened 'img_0.fits' Image size .............................. : (512, 512) pixels Number of channels ...................... : 1 Number of Stokes parameters ............. : 4 Beam shape (major, minor, pos angle) .... : (0.0047, 0.00298, 80.8) degrees Frequency of image ...................... : 1424.500 MHz Number of blank pixels .................. : 0 (0.0%) Flux from sum of (non-blank) pixels ..... : 0.087 Jy --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (65, 22) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map .... : (0.00005, 0.00028) Jy/beam Value of background mean ................ : -0.0 Jy/beam --> Expected 5-sigma-clipped false detection rate < fdr_ratio --> Using sigma-clipping ('hard') thresholding Minimum number of pixels per island ..... : 10 Number of islands found ................. : 4 Fitting islands with Gaussians .......... : [=---] 1/4 Fitting islands with Gaussians .......... : [===] 3/4 Fitting islands with Gaussians .......... : [====] 4/4 Total number of Gaussians fit to image .. : 4 Total flux density in model ............. : 0.070 Jy --> Grouping Gaussians into sources Number of sources formed from Gaussians : 4 2015/03/02 21:16:46 INFO: writing PyBDSM gaul catalog --> Wrote ASCII file '3C147-CD-LO-spw0-s3_pybdsm.lsm.html.gaul'

o-smirnov commented 9 years ago

OK, try running pybdsm interactively and give it the same command. Something else must've changed, you can look at my commit, https://github.com/ska-sa/pyxis/commit/40a67c3b3d78d6714d80d7acc9f4df35581d355b, it has nothing to do with pybdsm invocation per se.

modhurita commented 9 years ago

Pybdsm in interactive mode works as expected:

Non-default input parameters in interactive mode:

PROCESS_IMAGE: Find and measure sources in an image. filename ....... 'img_0.fits' polarisation_do ....... True thresh_isl ............ 15.0 thresh_pix ............ 50.0

Pybdsm output:

--> Opened 'img_0.fits' Image size .............................. : (512, 512) pixels Number of channels ...................... : 1 Number of Stokes parameters ............. : 4 Beam shape (major, minor, pos angle) .... : (0.0047, 0.00298, 80.8) degrees Frequency of image ...................... : 1424.500 MHz Number of blank pixels .................. : 0 (0.0%) Flux from sum of (non-blank) pixels ..... : 0.094 Jy --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (65, 22) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map (I) : (0.00005, 0.00027) Jy/beam Value of background mean (I) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (Q) : (0.00003, 0.00006) Jy/beam Value of background mean (Q) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (U) : (0.00003, 0.00006) Jy/beam Value of background mean (U) ............ : -0.0 Jy/beam Derived rms_box (box size, step size) ... : (65, 22) pixels Min/max values of background rms map (V) : (0.00003, 0.00005) Jy/beam Value of background mean (V) ............ : -0.0 Jy/beam --> Expected 5-sigma-clipped false detection rate < fdr_ratio --> Using sigma-clipping ('hard') thresholding Minimum number of pixels per island ..... : 10 Number of islands found ................. : 4 Fitting islands with Gaussians .......... : [====] 4/4 Total number of Gaussians fit to image .. : 4 Total flux density in model ............. : 0.070 Jy --> Grouping Gaussians into sources Number of sources formed from Gaussians : 4 --> Checking PI image for new sources --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (51, 17) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map .... : (0.00002, 0.00005) Jy/beam Value of background mean ................ : 4e-05 Jy/beam Minimum number of pixels per island ..... : 10 Number of islands found ................. : 0 --> Grouping Gaussians into sources Number of sources formed from Gaussians : 0 New sources found in PI image ........... : 0 (4 total) Calculating polarisation properties .... : [====] 4/4

However, invoking it through pyxis causes it to disregard polarization:

In [10]: lsm.pybdsm_search(image='img_0.fits',thresh_pix=50,thresh_isl=15,pol=True) 2015/03/03 16:21:58 INFO: running PyBDSM process_image(img_0.fits,polarisation_do=True,quiet=True,thresh_isl=15,thresh_pix=50) --> Opened 'img_0.fits' Image size .............................. : (512, 512) pixels Number of channels ...................... : 1 Number of Stokes parameters ............. : 4 Beam shape (major, minor, pos angle) .... : (0.0047, 0.00298, 80.8) degrees Frequency of image ...................... : 1424.500 MHz Number of blank pixels .................. : 0 (0.0%) Flux from sum of (non-blank) pixels ..... : 0.094 Jy --> Calculating background rms and mean images Derived rms_box (box size, step size) ... : (65, 22) pixels --> Variation in rms image significant --> Using 2D map for background rms --> Variation in mean image not significant --> Using constant background mean Min/max values of background rms map .... : (0.00005, 0.00027) Jy/beam Value of background mean ................ : -0.0 Jy/beam --> Expected 5-sigma-clipped false detection rate < fdr_ratio --> Using sigma-clipping ('hard') thresholding Minimum number of pixels per island ..... : 10 Number of islands found ................. : 4 Fitting islands with Gaussians .......... : [====] 4/4 Total number of Gaussians fit to image .. : 4 Total flux density in model ............. : 0.070 Jy --> Grouping Gaussians into sources Number of sources formed from Gaussians : 4 2015/03/03 16:22:01 INFO: writing PyBDSM gaul catalog --> Wrote ASCII file '3C147-spw0-s1_pybdsm.lsm.html.gaul'

Yesterday I did a git-pull after you edited lsm.py. I was on the devel branch, now I am on master, and there were many changes that got pulled - not sure if that made a difference. Also, I had to comment out lofarsoft-trunk-related things from the .bashrc file before pybdsm in interactive mode worked properly for polarization (it was giving errors earlier).

modhurita commented 9 years ago

PyBDSM in polarized mode is again not working for me; I get this error:

ERROR: IndexError: list index out of range [lofar.bdsm.make_residimage]

Here is the entire output:

In [5]: lsm.pybdsm_search(image='VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits',pol=True)
2015/10/30 17:25:10 INFO: running PyBDSM process_image(VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits,polarisation_do=True,quiet=True)
--> Opened 'VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits'
Image size .............................. : (2048, 2048) pixels
Number of channels ...................... : 1
Number of Stokes parameters ............. : 4
Beam shape (major, minor, pos angle) .... : (4.87071e-03, 4.53302e-03, -88.1) degrees
Frequency of image ...................... : 1474.845 MHz
Number of blank pixels .................. : 0 (0.0%)
Flux from sum of (non-blank) pixels ..... : 1.049 Jy
--> Calculating background rms and mean images
Derived rms_box (box size, step size) ... : (205, 68) pixels
--> Variation in rms image significant
--> Using 2D map for background rms
--> Variation in mean image not significant
--> Using constant background mean
Min/max values of background rms map (I)  : (0.00009, 0.00014) Jy/beam
Value of background mean (I) ............ : -1e-05 Jy/beam
Derived rms_box (box size, step size) ... : (205, 68) pixels
Min/max values of background rms map (Q)  : (0.00000, 0.00000) Jy/beam
Value of background mean (Q) ............ : -1e-05 Jy/beam
Derived rms_box (box size, step size) ... : (205, 68) pixels
Min/max values of background rms map (U)  : (0.00000, 0.00000) Jy/beam
Value of background mean (U) ............ : -1e-05 Jy/beam
Derived rms_box (box size, step size) ... : (205, 68) pixels
Min/max values of background rms map (V)  : (0.00001, 0.00003) Jy/beam
Value of background mean (V) ............ : -1e-05 Jy/beam
--> Expected 5-sigma-clipped false detection rate < fdr_ratio
--> Using sigma-clipping ('hard') thresholding
Minimum number of pixels per island ..... : 15
Number of islands found ................. : 11
Total number of Gaussians fit to image .. : 38
Total flux density in model ............. : 1.381 Jy
--> Grouping Gaussians into sources
Number of sources formed from Gaussians   : 11

--> Checking PI image for new sources
--> Calculating background rms and mean images
Derived rms_box (box size, step size) ... : (204, 68) pixels
--> Variation in rms image significant
--> Using 2D map for background rms
--> Variation in mean image not significant
--> Using constant background mean
Min/max values of background rms map .... : (0.00000, 0.00000) Jy/beam
Value of background mean ................ : 0.0 Jy/beam
Minimum number of pixels per island ..... : 15
Number of islands found ................. : 146
Total number of Gaussians fit to image .. : 136
--> Grouping Gaussians into sources
Number of sources formed from Gaussians   : 130
New sources found in PI image ........... : 110 (121 total)
ERROR: IndexError: list index out of range [lofar.bdsm.make_residimage]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-5-55e9e6d330ee> in <module>()
----> 1 lsm.pybdsm_search(image='VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits',pol=True)

/home/mitra/MeqTrees/pyxis/Pyxides/lsm.pyc in pybdsm_search(image, output, pol, select, center, threshold, pbexp, **kw)
     63   info("running PyBDSM process_image($image,%s)"%",".join(sorted([ "%s=%s"%x for x in opts.iteritems() ])));
     64   from lofar import bdsm
---> 65   img = bdsm.process_image(image,**opts);
     66   info("writing PyBDSM gaul catalog");
     67   img.write_catalog(outfile=gaul,format='ascii',catalog_type='gaul',clobber=True);

/usr/lib/python2.7/dist-packages/lofar/bdsm/__init__.pyc in process_image(input, **kwargs)
    241     # Now process it. Any kwargs specified by the user will
    242     # override those read in from the parameter save file or dictionary.
--> 243     img.process(**kwargs)
    244     return img

/usr/lib/python2.7/dist-packages/lofar/bdsm/image.pyc in process(self, **kwargs)
    124         """Process Image object"""
    125         import interface
--> 126         success = interface.process(self, **kwargs)
    127         return success
    128 

/usr/lib/python2.7/dist-packages/lofar/bdsm/interface.pyc in process(img, **kwargs)
     51         img, op_chain = get_op_chain(img)
     52         if op_chain is not None:
---> 53             _run_op_list(img, op_chain)
     54             img._prev_opts = img.opts.to_dict()
     55         return True

/usr/lib/python2.7/dist-packages/lofar/bdsm/__init__.pyc in _run_op_list(img, chain)
    147                 answ = raw_input_no_history(prompt)
    148         op.__start_time = time()
--> 149         op(img)
    150         op.__stop_time = time()
    151         gc.collect()

/usr/lib/python2.7/dist-packages/lofar/bdsm/make_residimage.pyc in __call__(self, img)
     48               isl = img.islands[g.wisland_id]
     49             else:
---> 50               isl = img.islands[g.island_id]
     51             b = self.find_bbox(thresh*isl.rms, g)
     52 

IndexError: list index out of range

I am working in this directory: elwood: /home/mitra/MuellerMatrixCalibration.

o-smirnov commented 9 years ago

Hmmm, so what's changed from before? Is it perhaps that you now have a multi-channel cube, whereas before it was single channel? Could you please test whether it works with single-channel images? Also, will it work if you enable spectral mode?

modhurita commented 9 years ago

It is a single channel image - it is this image: /home/mitra/MuellerMatrixCalibration/VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits. It wasn't working for multichannel images, so I checked to see if it at least works for a single-channel, all-Stokes image - it doesn't.

o-smirnov commented 9 years ago

But this step works in the pipeline, so what's different?

On Fri, Oct 30, 2015 at 5:50 PM, Modhurita Mitra notifications@github.com wrote:

It is a single channel image - it is this image: /home/mitra/MuellerMatrixCalibration/VLAC_1GHZ_BW_50MHZ_CHANWIDTH_2_IQUV.MS.CORRECTED_DATA.channel.1ch.restored.fits. It wasn't working for multichannel images, so I checked to see if it at least works for a single-channel image - it doesn't.

— Reply to this email directly or view it on GitHub https://github.com/ska-sa/pyxis/issues/30#issuecomment-152563188.

modhurita commented 9 years ago

Which pipeline do you mean? The 3C147 pipeline doesn't perform source extraction on any full-polarization image.

modhurita commented 9 years ago

Also, can tigger-convert be modified so that one can choose to run it in either total intensity (dividing measured Stokes I by beam gain) or full polarization (multiplying measured Stokes vector by inverse of Mueller matrix) mode? Currently, the 3C147 pipeline is producing apparent LSMs which have only I, but the intrinsic LSMs produced by tigger-convert have non-zero I and Q but zero U and V.