Adam Ginsburg's blog

• ALMA Cycle 9 corrupted zip fix

  If you encounter a bug like this:

  

  "XML (syntax or validity) error" because you hit issue C1_032, the fix is really easy:

  import zipfile
  import os
  
  # not strictly necessary, but important if you want to avoid filling your directory with files
  os.mkdir('recovery')
  os.chdir('recovery')
  
  # replace with your AOT file name
  with zipfile.ZipFile('../corrupted.aot', 'r') as aot:
      # this will extract all the files into the current directory
      aot.extractall()
      files = aot.filelist
  
  with zipfile.ZipFile('../recovered.aot', 'w') as aot:
      for fh in files:
          aot.write(fh.filename)
  

  At least, this worked on the one case where I got bitten.

  This fix is based on Erik Rosolowsky's ALMA proposal tools

---

• Determing ALMA Array Configurations

  ALMA's named array configurations are not accessible through CASA. They are available from https://almascience.eso.org/observing/observing-configuration-schedule/prior-cycle-observing-and-configuration-schedule

  This snippet will let you grab the tables and create a mapping from date to array config name.

  However, it's a HUGE waste of time, because these data are stored directly in the MS!

  This is the right way:

  tb.open(vis+'/ASDM_EXECBLOCK')
  tb.getcol('configName')
  # array(['C43-3'], dtype='<U16')
  

  This is the hacky, reconstructed from the website, bad way:

  import requests
  from bs4 import BeautifulSoup
  from astropy.table import Table, vstack
  from astropy.io import ascii
  from astropy.time import Time
  
  url = "https://almascience.eso.org/observing/observing-configuration-schedule/prior-cycle-observing-and-configuration-schedule"
  
  response = requests.get(url)
  response.raise_for_status()
  soup = BeautifulSoup(response.text)
  tables = soup.findAll('table', class_="grid listing")
  
  def clean_lines(soup, badrows=['Long Baseline Campaign',
                                 'February Maintenance Period',
                                 'End of Cycle',
                                 'Engineering/Software',
                                 'Engineering/Software Time',
                                     ]):
      for tr in soup('tr'):
          if any(bad in str(tr) for bad in badrows):
               _=tr.extract()
      return soup
  
  tables = soup.findAll('table', class_="grid listing")
  tables = [ascii.read(str(clean_lines(tbl)), format='html') for tbl in tables]
  tables[3] = tables[3][:-1]
  tables[3]['Block'] = tables[3]['Block'].astype('int')
  tables[4]['maximumrecoverablescale2(")'] = tables[4]['maximumrecoverablescale2(")'].astype('str')
  
  stacked = vstack(tables)
  start_times = Time(stacked['Start date'])
  end_times = Time(stacked['End date'])
  

  Then, say you have a list of measurement sets mses, you can look up the array configuration for each date and field. The choice of fields[0] here only makes sense if your data had a single target; it's a bad choice if there are multiple targets in a scheduling block:

  import json
  from casatools import msmetadata
  msmd = msmetadata()
  
  results = {}
  
  for vis in mses:
      msmd.open(mses)
      obstime = Time(msmd.timerangeforobs(0)['begin']['m0']['value'], format='mjd')
      fieldnames = np.array(msmd.fieldnames())
      fields = np.unique(fieldnames[msmd.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')])
      msmd.close()
  
      array_config = stacked[(obstime > start_times) & (obstime < end_times)]['Approx\xa0Config.']
      if fields[0] in results:
          results[fields[0]][obstime.strftime('%Y-%m-%d')] = array_config[0]
      else:
          results[fields[0]] = {obstime.strftime('%Y-%m-%d'): array_config[0]}
  
  with open('array_configurations.json', 'w') as fh:
      json.dump(results, fh)
  

---

• FeO Nondetection in Orion

  I had a Cycle 3 ALMA program 2015.1.00262.S: Digging for rusty bullets at an explosion site that aimed to detect ferrous oxide, FeO in the shocked outflows of the Orion nebula.

  I haven't published anything on this data set because it resulted in a fairly boring non-detection. You can find the source code to retrieve and process the data from here. The nondetection is several orders of magnitude deeper than previous ones, but to demonstrate this definitively, we need to do some more modeling work, which my co-I's and I have decided isn't worth the time, since we're all flooded with very exciting new ALMA data that has interesting detections in it (see salts soon).

  In brief, we were looking at a location where we know molecules (H2, molecular hydrogen) and gas-phase iron (Fe II, as seen in in the [Fe II] 1.64 micron lines) are approximately coincident. The most likely (and broadly accepted) explanation is that dust particles are being destroyed (sputtered) in the high-velocity ($>30$ km/s, up to about $100-200$ km/s) shocks produced in an outflow. Since iron in the interstellar medium mostly lives in the cores of dust grains, we thought it might be released in a molecular form detectable by ALMA before being split into an atomic, ionized form.

  Despite having a good location and a sensitive observation, we didn't detect any. It might mean that iron is mostly in a different molecular form, or it might mean that when iron gets sputtered out of grains, it never goes into a molecular form at all, instead dissociating straight into atomic form.

  There are a few things that, with infinite time, I would do with this data set:

  1. Come up with a formal upper limit on the FeO abundance.
  2. Compare the lines we did detect to the outflows and see if we can better constrain their physical conditions
  3. Catalog the stars and publish the catalog for comparison to other data sets. (actually, this is done: this table contains Gaussian fits to each of the stars in the continuum image (source).)

  These figure show the apertures we used to extract spectra and search for the lines, then the spectra of different species for one of the apertures. The background image is the H2 (orange) + FeII (cyan) + K-band (mostly white) I made for my 2015 paper with John Bally.

  

---

• The extragalactic view of W51 at 220 and 290 GHz

  These are APEX data covering 5.5' x 5.5' centered on W51 Main: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/RYEANM These are calibrated to T_A* units. APEX telescope efficiencies can be found here: http://www.apex-telescope.org/telescope/efficiency/ and are approximately 0.75 at the observed frequencies.

  W51 is at a distance of 5.1 kpc, so this field of view is equivalent to 8.1 pc, or 1.7" at 1 Mpc or 0.33" at 5 Mpc.

  The linewidth is 12.7-15.2 km/s (H2CO 4-3, CO 2-1) FWHM averaged over this full field.

  

  Images of the spectra, with some lines identified, are shown here (open them in a new tab to see the full resolution). The baselines aren't great, but it's easy to fit a local linear baseline around any of the fairly narrow lines.

  

  The same images, but not zoomed in the y-direction. Note that 12CO 2-1 is very near the edge of the band; it peaks around 20 K:

  

  Scripts to obtain the data and process them are here: https://github.com/keflavich/W51_APEX_H2CO

  See the related project by Watanabe et al looking at W51 in the 3mm band: http://adsabs.harvard.edu/abs/2017ApJ...845..116W

---

• More CASA simulation & debugging

  (this post written from snowy-ish ALMA; just cloudy at the OSF though)

  (tl; dr version: sm.setvp doesn't work if you have TELESCOP="ALMA" in your header)

  I'm trying to repeat an experiment very similar to previous simulation work in order to quantify my completeness and largest angular scales by injecting synthetic signal into the real UV data.

  However, as you might expect, it doesn't work as it did last time, even when I (afaict) copy and paste the code directly.

  Here's the new error:

  2016-04-25 02:57:21 SEVERE  Simulator::createSkyEquation() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2200)        Caught exception: Transformations to/from frame "Undefined" are not possible.
  2016-04-25 02:57:21 SEVERE  Simulator::predict() (file /Users/rpmbuild/gradle/workdir/casasources/release-4_5/code/synthesis/MeasurementEquations/Simulator.cc, line 2118)  Failed to create SkyEquation
  Result of predict: False
  

  It's strange because there is nothing obviously undefined.

  I'm comparing 2 FITS files that should behave identically. The good header lacks these keywords that the bad header has:

  ['PV2_1',
   'PV2_2',
   'OBSERVER',
   'DATE-OBS',
   'OBSRA',
   'OBSDEC',
   'OBSGEO-X',
   'OBSGEO-Y',
   'OBSGEO-Z']
  

  while the good header has these and the bad lacks them:

  ['PC03_01',
   'PC04_01',
   'PC03_02',
   'PC04_02',
   'PC01_03',
   'PC02_03',
   'PC03_03',
   'PC04_03',
   'PC01_04',
   'PC02_04',
   'PC03_04',
   'PC04_04',
   'CREATOR',
   'INSTRUME',
   'BAND',
   'PROPOSAL',
   'PRTITLE',
   'OBSID001',
   'OBSID002',
   'HISTORY']
  

  The only significant difference is the PC values, unless CASA is doing something with the observatory location information.

  The only difference in common keywords is numerical, which seems very unlikely to cause this issue.

  I cannot discern any useful difference between these two CASA .image headers:

  2016-04-25 03:08:55 INFO imhead     ##########################################
  2016-04-25 03:08:55 INFO imhead     ##### Begin Task: imhead             #####
  2016-04-25 03:08:55 INFO imhead     imhead(imagename="../perseus_synth/perseus_250_to_w51.image",mode="summary",hdkey="",hdvalue="",verbose=True)
  2016-04-25 03:08:55 INFO ImageAnalysis
  2016-04-25 03:08:55 INFO ImageAnalysis      Image name       : perseus_250_to_w51.image
  2016-04-25 03:08:55 INFO ImageAnalysis      Object name      : perseus 04
  2016-04-25 03:08:55 INFO ImageAnalysis      Image type       : PagedImage
  2016-04-25 03:08:55 INFO ImageAnalysis      Image quantity   : Intensity
  2016-04-25 03:08:55 INFO ImageAnalysis      Pixel mask(s)    : mask0
  2016-04-25 03:08:55 INFO ImageAnalysis      Region(s)        : None
  2016-04-25 03:08:55 INFO ImageAnalysis      Image units      : Jy/pixel
  2016-04-25 03:08:55 INFO ImageAnalysis      Restoring Beam   : 0.465804 arcsec, 0.465804 arcsec, 0 deg
  2016-04-25 03:08:55 INFO ImageAnalysis
  2016-04-25 03:08:55 INFO ImageAnalysis      Direction reference : J2000
  2016-04-25 03:08:55 INFO ImageAnalysis      Spectral  reference : Undefined
  2016-04-25 03:08:55 INFO ImageAnalysis      Velocity  type      : RADIO
  2016-04-25 03:08:55 INFO ImageAnalysis      Rest frequency      : 2.33947e+11 Hz
  2016-04-25 03:08:55 INFO ImageAnalysis      Telescope           : Herschel
  2016-04-25 03:08:55 INFO ImageAnalysis      Observer            : UNKNOWN
  2016-04-25 03:08:55 INFO ImageAnalysis      Date observation    : UNKNOWN
  2016-04-25 03:08:55 INFO ImageAnalysis
  2016-04-25 03:08:55 INFO ImageAnalysis      Axis Coord Type      Name             Proj Shape Tile   Coord value at pixel    Coord incr Units
  2016-04-25 03:08:55 INFO ImageAnalysis      ------------------------------------------------------------------------------------------------
  2016-04-25 03:08:55 INFO ImageAnalysis      0    0     Direction Right Ascension   TAN  2370  158  19:23:41.765  1099.00 -1.552680e-01 arcsec
  2016-04-25 03:08:55 INFO ImageAnalysis      1    0     Direction Declination       TAN  2500  250 +14.30.45.850  1552.00  1.552680e-01 arcsec
  2016-04-25 03:08:55 INFO ImageAnalysis      2    2     Spectral  Frequency                 1    1   2.33947e+11     0.00  1.000000e+09 Hz
  2016-04-25 03:08:55 INFO ImageAnalysis                           Velocity                                     0     0.00 -1.281456e+03 km/s
  2016-04-25 03:08:55 INFO ImageAnalysis      3    1     Stokes    Stokes                    1    1             I
  2016-04-25 03:08:55 INFO imhead     ##### End Task: imhead               #####
  2016-04-25 03:08:55 INFO imhead     ##########################################
  2016-04-25 03:09:07 INFO imhead
  2016-04-25 03:09:07 INFO imhead     ##########################################
  2016-04-25 03:09:07 INFO imhead     ##### Begin Task: imhead             #####
  2016-04-25 03:09:07 INFO imhead     imhead(imagename="../simulations/simimage_firsttest.image",mode="summary",hdkey="",hdvalue="",verbose=False)
  2016-04-25 03:09:07 INFO ImageAnalysis
  2016-04-25 03:09:07 INFO ImageAnalysis      Image name       : simimage_firsttest.image
  2016-04-25 03:09:07 INFO ImageAnalysis      Object name      :
  2016-04-25 03:09:07 INFO ImageAnalysis      Image type       : PagedImage
  2016-04-25 03:09:07 INFO ImageAnalysis      Image quantity   : Intensity
  2016-04-25 03:09:07 INFO ImageAnalysis      Pixel mask(s)    : None
  2016-04-25 03:09:07 INFO ImageAnalysis      Region(s)        : None
  2016-04-25 03:09:07 INFO ImageAnalysis      Image units      : Jy/beam
  2016-04-25 03:09:07 INFO ImageAnalysis      Restoring Beam   : 0.21023 arcsec, 0.192666 arcsec, 79.8538 deg
  2016-04-25 03:09:07 INFO ImageAnalysis
  2016-04-25 03:09:07 INFO ImageAnalysis      Direction reference : J2000
  2016-04-25 03:09:07 INFO ImageAnalysis      Spectral  reference : Undefined
  2016-04-25 03:09:07 INFO ImageAnalysis      Velocity  type      : RADIO
  2016-04-25 03:09:07 INFO ImageAnalysis      Rest frequency      : 2.33947e+11 Hz
  2016-04-25 03:09:07 INFO ImageAnalysis      Pointing center     :  19:23:41.629000  +14.30.42.380000
  2016-04-25 03:09:07 INFO ImageAnalysis      Telescope           : ALMA
  2016-04-25 03:09:07 INFO ImageAnalysis      Observer            : keflavich
  2016-04-25 03:09:07 INFO ImageAnalysis      Date observation    : 2015/04/23/09:47:44
  2016-04-25 03:09:07 INFO ImageAnalysis      Telescope position: [2.22514e+06m, -5.44031e+06m, -2.48103e+06m] (ITRF)
  2016-04-25 03:09:07 INFO ImageAnalysis
  2016-04-25 03:09:07 INFO ImageAnalysis      Axis Coord Type      Name             Proj Shape Tile   Coord value at pixel    Coord incr Units
  2016-04-25 03:09:07 INFO ImageAnalysis      ------------------------------------------------------------------------------------------------
  2016-04-25 03:09:07 INFO ImageAnalysis      0    0     Direction Right Ascension   TAN  3072  192  19:23:41.629  1536.00 -5.000000e-02 arcsec
  2016-04-25 03:09:07 INFO ImageAnalysis      1    0     Direction Declination       TAN  3072  192 +14.30.42.380  1536.00  5.000000e-02 arcsec
  2016-04-25 03:09:07 INFO ImageAnalysis      2    1     Stokes    Stokes                    1    1             I
  2016-04-25 03:09:07 INFO ImageAnalysis      3    2     Spectral  Frequency                 1    1   2.33947e+11     0.00  1.000000e+09 Hz
  2016-04-25 03:09:07 INFO ImageAnalysis                           Velocity                                     0     0.00 -1.281456e+03 km/s
  2016-04-25 03:09:07 INFO imhead     ##### End Task: imhead               #####
  2016-04-25 03:09:07 INFO imhead     ##########################################
  

  if anything, the ALMA file includes extra data. Removing this data has no effect.

  The only really substantial difference left is the order of the axes, which is not what I put in for the ALMA data: it puts Stokes as 2 instead of 3 despite the fact that CTYPE4='STOKES' and CTYPE3='FREQ'.

  By doing a round trip FITS->image->FITS->image, I was able to fix the order, but that is not the underlying problem, apparently.

  Could the missing pixel mask be at fault? It doesn't jive at all with the error message, so I doubt it. OBSDATE is not the problem.

  Turns out, the error is that ALMA is the telescope. I added the header entry:

  ffile[0].header['TELESCOP'] = 'NotReal'
  

  and it Just Worked. It is approaching the end of my shift, so I think it is now acceptable to say "what kind of !#$@#$!@% up #@$$@% is that?!".

---

• 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:

  

---

• Self-calibration of Sgr B2 Mosaic

  The Sgr B2 ALMA data from project 2013.1.00269.S has imaging problems around the brightest sources, especially Sgr B2 M. Crystal Brogan & Todd Hunter gave some useful advice, essentially the following:

  • the line forest is causing trouble, so line-containing regions of the spectrum must be excluded when mapping
  • there is some nasty curvature in the spectral baseline that needs to be removed probably by flagging out some of the data
  • Self-calibration is needed to exceed a dynamic range ~100

  Todd developed a fitcontinuum.py script to address the first point, but I ended up doing this custom (script) because of the problem introduced by the second point.

  To exclude the lines, I had to use frequency-based selection, not channel based, because my spectral windows from different dates were offset by 10 MHz. The approach to use in CASA is ugly: you flag the line-containing data, then split (and in splitting, average) the continuum into another MS. From there, a phase-only selfcal is straigtforward (script).

  

---

Page 1 / 2 »