Adam Ginsburg's blog

• Which one is real?

  I've "successfully" combined 7m and 12m data, and I've also "successfully" processed 12m data. In both cases, I've self-calibrated the long-baseline 12m data.

  The question now is: why are these images so different? There are features in the 12m-only data that are a factor of ~2 brighter than in the 12m+7m data. That should not happen.

  This image shows the problem. 7m+12m left, 12m-only right:

  

  Note the sharp features on the left and the small tail (filament) toward the top right that appear in the 12m but not in the 7m+12m data.

  Which features are real? How do we decide?

  I really want to believe the 12m-only data, but I am very worried that it is artificially enhanced by sidelobes from the e8 filament:

  

  There is a ~20-sigma feature between e8 and e2 along the e8 filament axis, but it is completely absent in the 7m+12m data. That leads me to believe that it is not real, but at such high significance that is very worrisome. I'm also concerned by the loss of the compact source at the bottom of the e8 filament and to the left: these seem unlikely to be affected by the filament.

---

• CLEAN has extreme sensitivity to the threshold

  CLEAN diverges sometimes. The point at which divergence begins is very difficult to identify and unfortunately it looks like the images get mostly better as you slowly approach the toxic threshold:

  

  The precise clean command used, with only threshold changing, is:

  2016-03-10 07:49:28     INFO    tclean::::      tclean(vis="w51_contvis_selfcal_4.ms",selectdata=True,field="",spw="",timerange="",
  2016-03-10 07:49:28     INFO    tclean::::+             uvrange="",antenna="",scan="",observation="",intent="",
  2016-03-10 07:49:28     INFO    tclean::::+             datacolumn="corrected",imagename="selfcal_allspw_selfcal_4ampphase_mfs_tclean_deeper_4mJy",imsize=[3072, 3072],cell="0.05arcsec",phasecenter="J2000 19:23:41.629000 +14.30.42.38000",
  2016-03-10 07:49:28     INFO    tclean::::+             stokes="I",projection="SIN",startmodel="",specmode="mfs",reffreq="",
  2016-03-10 07:49:28     INFO    tclean::::+             nchan=-1,start="",width="",outframe="LSRK",veltype="radio",
  2016-03-10 07:49:28     INFO    tclean::::+             restfreq=[],interpolation="linear",gridder="mosaic",facets=1,wprojplanes=1,
  2016-03-10 07:49:28     INFO    tclean::::+             aterm=True,psterm=False,wbawp=True,conjbeams=True,cfcache="",
  2016-03-10 07:49:28     INFO    tclean::::+             computepastep=360.0,rotatepastep=360.0,pblimit=0.4,normtype="flatnoise",deconvolver="clark",
  2016-03-10 07:49:28     INFO    tclean::::+             scales=[],nterms=2,restoringbeam=[],outlierfile="",weighting="briggs",
  2016-03-10 07:49:28     INFO    tclean::::+             robust=-2.0,npixels=0,uvtaper=[],niter=100000,gain=0.1,
  2016-03-10 07:49:28     INFO    tclean::::+             threshold="4mJy",cycleniter=-1,cyclefactor=1.0,minpsffraction=0.05,maxpsffraction=0.8,
  2016-03-10 07:49:28     INFO    tclean::::+             interactive=False,mask="",overwrite=True,savemodel="modelcolumn",calcres=True,
  2016-03-10 07:49:28     INFO    tclean::::+             calcpsf=True,parallel=False)
  

  It looks like CLEAN needs a different type of threshold to finish on, perhaps something requiring each subsequent component to be some value less than the previous component. If each component is forced to be less than the previous one, the clean can never diverge, though clean would stop unexpectedly if subtracting a component increased the amplitude of some other part of the map (i.e., if the next brightest source is sitting in a negative bowl). Still, this might be a rarer occurrence than the divergence that I see all the time, which may be related to gridding (given the sharp grid pattern seen in panel 4 of the linked image).

---

• Simulated Observations of Perseus in W51

  The W51 and Sgr B2 ALMA mosaics contain ~50 and ~150 pointlike sources, respectively. However, they have pretty miserable high-noise regions such that determining completeness is a serious challenge.

  To determine completeness in a means comparable to other observations, I have therefore been using the Gould's Belt Survey Perseus data as an input to do synthetic observations.

  It turns out this is technically challenging because CASA does some strange things when importing files.

  I finally got a version working.

  Important details when reading in a FITS file:

  • The header must be fully populated with a 4-dimensional WCS. If stokes or frequency are missing, CASA will choke with uninformative errors such as:

    2016-02-26 11:25:24 SEVERE Simulator::createSkyEquation() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2200) Caught exception: (/Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc : 2266) Failed AlwaysAssert spectralIndex>=0 2016-02-26 11:25:24 SEVERE Simulator::predict() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2118) Failed to create SkyEquation

    and

    2016-02-26 11:27:36 SEVERE Simulator::predict() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2118) Caught exception: (/Users/rpmbuild/gradle/workdir/casasources/release-4_5/casacore/coordinates/Coordinates/CoordinateSystem.cc : 793) Failed AlwaysAssert which < nCoordinates() && coordinates_p[which]->type() == Coordinate::STOKES

  • Flux units should be in Jy/beam. Other flux units may be acceptable (e.g., MJy/sr, Jy/pixel), but their behavior is not entirely predictable. I did not get out the expected values when I used MJy/sr as input, but I did not pursue further why this happened. Using Jy/beam, and ensuring that the units were correctly set to Jy/beam in the file (by scaling them in the FITS image before importing), I got the expected results.

  • CASA apparently does not read beam information from headers. You need to include the beam information when doing importfits, e.g.:

    beam=["{0}deg".format(hdr['BMAJ']),

    "{0}deg".format(hdr['BMIN']), "{0}deg".format(hdr['BPA'])]

  • CASA will happily read a FITS file into a .image that cannot be used by other casa routines. You will not get an exception until you try a later step.

  Once the FITS image is read into a CASA .image, if it has the CASA-expected units and dimensions, you can use sm.predict to fourier transform the data and sample it onto the UV grid in your .ms file:

  sm.openfromms("file.ms")
  sm.setvp() # "set voltage pattern": not clear if this is needed
  success = sm.predict("myfile.image")
  sm.done()
  sm.close()
  

  This (seems to) populate the data column, not the model column as I had originally expected. Therefore, be warned! DO NOT DO THIS TO YOUR ORIGINAL DATA! It will probably overwrite it.

  Once your .ms file is populated, you can clean it as normal. Some examples:

  

---

• CASA Simulation

  More frustration with CASA:

  • If you try to load a FITS image into a CASA .image, beware that it is likely to have incorrect units. Round-tripping between .image and .fits doesn't work well, because CASA always assumes that CDELT is in radians even when it is explicitly set to be degrees.

  • EDIT: OK, this isn't true in all cases, but if you use imhead to put a keyword value, it will assume radians by default. That can be overcome by explicitly setting the units:

    imhead(perseus_casa_image, mode='put', hdkey='CDELT1', hdvalue={'value':hdr['CDELT1'], 'unit':'deg'})

    I was able to produce this awkward behavior without imhead as well, though I cannot repeat that now.

  The simulator gives this error message if I try to load an image:

  sm.setvp() sm.predict(perseus_casa_image)

  2016-02-25 09:26:34 SEVERE Simulator::predict() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2118) Caught exception: TiledShape: #elements in shape and tileShape differ

  My best guess is that CASA doesn't understand how to handle 2D images with 3D headers. However, the sensible workaround of expanding the image to 3D does not work.

  Re-importing after exporting gets yet another different image. After re-importing an image that was exported to FITS, I get a new error:

  Failed AlwaysAssert spectralIndex>=0 2016-02-25 09:37:35

  which is rather strange since the spectral index is a physical quantity.

  There is also the lovely message:

  2016-02-25 09:55:27 INFO importfits No usable restoring beam information found.

  for a header that contains these lines:

  BMAJ = 1.293900184840E-04 BMIN = 1.293900184840E-04 BPA = 0.000000000000E+00

  The image might finally be right, but now I can't find any way to convert it to a model. If I try to ft the model, I get:

  2016-02-25 10:54:54 SEVERE ft::imager::ft() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Imager.cc, line 4488) Exception: (/Users/rpmbuild/gradle/workdir/casasources/release-4_5/casacore/coordinates/Coordinates/CoordinateSystem.cc : 777) Failed AlwaysAssert which < nCoordinates() && coordinates_p[which]->type() == Coordinate::SPECTRAL

  which probably means that CASA doesn't know how to turn a single-band continuum image into a model for a multi-frequency .ms (which is all .ms's, including continuum, in the ALMA world).

---

• Retrieving and Selecting Data in CASA

  I learned a few important things at the 7m+12m workshop in Garching:

  1. clean treats input models in Jy/beam and Jy/pixel qualitatively (not just quantitatively) differently. Jy/pixel images are taken to be model images exactly, where Jy/beam are... somehow interpreted differently, maybe as clean component assemblies. I don't know the details, but this is a major caution.

  2. The tbtool is not very good for extracting data for plotting. The flags are not guaranteed to be in a sane shape, which will lead to crashes if you try tb.getcolumn('FLAG'). ms.msselect, however, allows you access to most of the normal selection routines used in clean, split, etc. It is tricky, though, because ms.select does not allow you access to the same selection. Also, while spw-based selection normally allows you to select channels, in ms.msselect it does not: instead you must use ms.selectchannel. See this example.

     Details about the selection are at this URL, but the documentation is incomplete and, in places, incorrect:

     http://www.aoc.nrao.edu/~sbhatnag/misc/msselection/FrequencySelection.html

     (for example, chanspec is not valid in msselect, or at least it doesn't do anything)

---

• A tclean vs clean comparison

  Comparing tclean to clean with the following parameters:

  clean(vis=vis0, imagename=myimagebase, field="", spw='',
        mode='mfs', outframe='LSRK', interpolation='linear',
        imagermode='mosaic', multiscale=multiscale, interactive=False,
        niter=1000, threshold='100mJy', imsize=imsize, minpb=0.5, cell=cell,
        phasecenter=phasecenter, weighting='briggs', usescratch=True,
        pbcor=True, robust=-2.0)
  
  tclean(vis=vis0, imagename=myimagebase, field="", spw='',
         outframe='LSRK', interpolation='linear', gridder='mosaic',
         scales=multiscale, interactive=False, niter=1000,
         threshold='100mJy', imsize=imsize, mask='auto-pb', specmode='mfs',
         cell=cell, phasecenter=phasecenter, weighting='briggs', robust=-2.0)
  

  Note the very bad cleaning instability creeping in in the clean image:

  

---

• CASA Selfcal: Good & Awful

  Self-calibration is necessary for some images to get a decent dynamic range or to see anything at all in the images. So, I did some tests using self-cal on the central field of Sgr B2 (source) and it worked great:

  

  That image is actually from self-calibrating the whole field, which also worked great. Except, it deleted all the flux outside of the bright regions:

  

  This is awful. It is just deleting flux! I had previously thought that the applyonly option would only apply solutions that it actually found, but that appears not to bt true. Instead, it seems that even phase solutions with no signal are being applied, which effectively removes everything from the image. I could have understood if this was happening to amplitude self-calibration, but I can't understand how it is happening with phase self-cal.

---

• CASA imaging of huge cube

  Imaging full cubes, i.e. every channel and every pixel, from ALMA data proved too much for the first few machines I tried it on, so I then tried splitting the cubes into chunks (source).

  For the most part, this works well, but there are some serious problems:

  

  Note the range from 233.8-234 GHz. It has a continuum level way above the rest. This is utter nonsense, and turns out to be due to some imaging error. The first image below shows a "good" chunk of the cube, the second image shows a "bad" chunk.

  

  It is clear that the second image was restore using a too-large beam, and this is verified by examining the header:

  In [8]: badcube.beam
  Out[8]: Beam: BMAJ=1.15963029861 arcsec BMIN=0.96223282814 arcsec BPA=-71.0310058594 deg
  
  In [9]: goodcube.beam
  Out[9]: Beam: BMAJ=0.408393889665 arcsec BMIN=0.229679107666 arcsec BPA=45.966835022 deg
  

  I don't yet know what is causing this error. When I try re-doing the clean with tclean, I get the following message, which is a hint:

  2016-01-25 09:13:04     WARN    task_tclean::SIImageStore::getPSFGaussian (file
  /var/rpmbuild/BUILD/casa-prerelease/casa-prerelease-4.5.0/code/synthesis/ImagerObjects/SIImageStore.cc,
  line 1262)      PSF is blank for[C139:P0] [C140:P0] [C141:P0] [C142:P0]
  [C143:P0] [C144:P0] [C145:P0] [C146:P0] [C147:P0] [C148:P0] [C149:P0] [C150:P0]
  [C151:P0] [C152:P0] [C153:P0] [C154:P0] [C155:P0] [C156:P0] [C157:P0] [C158:P0]
  [C159:P0] [C160:P0] [C161:P0] [C162:P0] [C163:P0] [C164:P0] [C165:P0] [C166:P0]
  [C167:P0] [C168:P0] [C169:P0] [C170:P0] [C171:P0] [C172:P0] [C173:P0] [C174:P0]
  [C175:P0] [C176:P0] [C177:P0] [C178:P0] [C179:P0] [C180:P0] [C181:P0] [C182:P0]
  [C183:P0] [C184:P0] [C185:P0] [C186:P0]
  

  Apparently tclean solves this problem! Instead of using a single beam for all channels, it creates a CASAMBM table in the FITS output and uses different beams at each channel. There must be genuinely bad data (probably an atmospheric absorption line) at the specified frequencies. At least now, that will come up more naturally, rather than spiking the data.

---

• Comparing tclean parameters

  Comparing clean methods with identical data & parameters (as far as possible): code.

  An example using natural weighting:

  

  None of the methods are acceptable. Multiscale is the worst, but hogbom is extremely bad: there are major ringing artifacts (stripes parallel to the main filament).

  Disturbingly, even though I have given multiple scales, there are clearly no non-pointlike scales in any of the model images, including mem, so I think there's a major problem with tclean:

  

  I'm re-doing it with uniform weighting, because that's the case I was actually interested in.

  

  The uniform tests come up with a clear and obvious answer: use clark clean, DO NOT use any other methods as they are all unstable.

---

• Self-cal with EVLA data

  I'm trying to self-calibrate my A+C array C-band data on W51, since we're in a clearly dynamic-range-limited regime. I've never had any luck with self-cal, though. The images always get worse, and I still don't really understand why.

  This is a clue, though:

  

  For many of the antennae, especially those in C array, self-calibration results in phase plots that look nice, with corrections confined to +/-30 degrees as you would hope and expect.

  These plots show that some antennae are going through dozens or hundreds of phase wraps. That is obviously unacceptable and incorrect: the original images actually looked pretty good, so it is obvious (to me, at least) that there should be no phase wrapping in the self-cal. I'm surprised at this, since each of these points should have S/N > 3 and there is a lot of consistency in these large phase shifts.

  My suspicion is still that the clean process has failed. I have been trying multiscale clean, since the emission is completely dominated by extended components. However, even if I set smallscalebias to different values, there are essentially no pointlike components being added to the clean model. This is wrong! There are genuine point sources. So I still suspect that this is a modeling problem and the phase wraps are being caused by bad models.

  I'm trying again with smaller pixels in the hope that a better match between the beam size and the pixel size will produce a better model. The only other thing I can think to try is increasing the solint, but given the high S/N already present, I doubt that's an important factor.

  Smaller pixels didn't help, but forcing clean to include point-source components by performing a second round of clean with multiscale=[] worked! The solutions are MUCH more reasonable (no more phase wrap every integration), and selfcal no longer diverges. I also increased the solution interval to 30s since it seems that the phase change timescale is significantly longer than an integration. I increased the S/N cutoff to 5 too.

  There are still some "noisy" phases, even in the 3rd iteration:

  

  I think more iterations may be necessary, but I will also increase the S/N threshold to 6 and solint to 60s to see if that helps.

  RMS improvements are evident but not spectacular:

  SPW4 Iteration 0: RMS=1.003e-04
  SPW4 Iteration 1: RMS=6.009e-05
  SPW4 Iteration 2: RMS=6.698e-05
  SPW4 Iteration 3: RMS=7.061e-05
  

  The bump in the last one is because I forgot to add a point clean for the last iteration.

  Perhaps not unfortunately, it looks like the improvements in image quality are tiny when self-cal is added. The only real hope now is that self-cal on all bands will result in an artifact-free final combined image.

  On the "noisy" phases, now that it looks like there are decent solutions, it may actually be worth keeping the flags generated by applycal. <10% of the data is being flagged out, which is negligible, and it probably is genuinely bad (or signal-free, anyway) in that case.

  EDIT: Looks better... or at least different... now:

  SPW4 Iteration 0: RMS=1.003e-04
  SPW4 Iteration 1: RMS=5.962e-05
  SPW4 Iteration 2: RMS=6.817e-05
  SPW4 Iteration 3: RMS=6.775e-05
  

---

« Page 5 / 51 »