Open in Colab: https://colab.research.google.com/github/casangi/casadocs/blob/2348961/docs/notebooks/simulation.ipynb


Simulations

The capability of simulating observations and data sets from VLA, ALMA, and other existing and future observatories is an important use-case for CASA. This not only allows the user to get an idea of the capabilities of these instruments for doing science, but also provides benchmarks for the performance and utility of the software to process “realistic” data sets (with atmospheric and instrumental effects). Simulations can also be used to tune parameters of the data reduction and therefore help to optimize the process. CASA can calculate visibilities (create a MeasurementSet) for any interferometric array, and calculate and apply calibration tables representing some of the most important corrupting effects.

Tasks available for simulating observations are:

  • simobserve - simulate and create custom synthetic MeasurementSets for an interferometric or total power observation

  • simanalyze - image and analyze simulated data set, including diagnostic images and plots

  • simalma - simulate an ALMA observation including multiple configurations of the 12-m interferometric array, the 7-m ACA, and total power measurements by streamlining the capabilities of both simobserve and simanalyze

Inside the Toolkit: The simulator methods are in the simulator tool sm. Many of the other CASA tools are helpful when constructing and analyzing simulations. Following general CASA practice, the greatest flexibility and functionality is available in the Toolkit, and the most commonly used procedures are bundled for convenience into the tasks.

Utility functions: The simutil python class contains numerous utility methods which can be used to facilitate simulations, especially when using the Toolkit.

Simulating interferometric observations using the simobserve and simanalyze tasks proceeds in the following steps:

  1. Make a model image or component list. The model is a representation of the sky brightness distribution that you would like to simulate observing (details on model specification in the simobserve documentation).

  2. Use the simobserve task to create a MeasurementSet (uv data) that would be measured by a telescope observing the specified input model of sky brightness. simobserve can also introduce corruption modeling thermal noise or atmospheric effects. The task can be called multiple times to generate e.g., observations from the same model using different array configurations. simobserve can also simulate total power observations, which can be combined with interferometric data in simanalyze.

  3. Image (grid, invert, and deconvolve) the simulated observation(s) with the simanalyze task. simanalyze can also compare the simulated image with your input (convolved with the output clean beam) and then calculate a “fidelity image” that indicates how well the simulated output matches the convolved input image. Alternately, you can create an image yourself with the tclean task, and then use simanalyze to compare that to the sky model input.

ALMA users, especially those less familiar with CASA, are encouraged to use the simalma task, which provides additional information on the multiple simobserve and simanalyze calls required to simulate an ALMA observation which may consist of 12m interferometric, 7m interferometric, and 12m total power data, and attempts to provide useful feedback on those different observation components, to help the user better understand the observing considerations.

Note that simobserve and simalma currently only handle linear polarization. For implementing full Stokes behavior, please see these examples.

More simulation examples can be found in the CASA guides http://casaguides.nrao.edu, under “Simulating Observations in CASA”, and in the Notebook examples of Simulation in CASA in CASA Docs. It is possible to run the steps independently and optionally, as long as you follow the simobserve and simanalyze conventions about filenames.

Tip: A list of antenna configuration files known to CASA is linked from the simulation CASA guides. On Unix, Linux, and Mac computers, you can usually also find this list yourself by typing, for instance, “locate alma.cycle4.1.cfg” and looking at other files in that directory.


ALMA simulations

The task simalma simulates an ALMA observation by ALMA 12-m, ACA-7m and total power arrays. It takes an input model image or a list of components, plus configurations of ALMA antennas (locations and sizes), and simulates a particular ALMA observation (specified by mosaic setup and observing cycles and times). The outputs are MeasurementSets. The task optionally generates images from the MeasurementSets.

Technically speaking, simalma internally calls simobserve and simanalyze as many times as necessary to simulate and analyze an ALMA observation. Some of the simulation and imaging parameters are automatically set to values typical of ALMA observations. Thus, it has a simpler task interface compared to simobserve plus simanalyze at the cost of limited flexibility. If the user wants to have more control of simulation setup, it is available by manually running simobserve and simanalyze multiple times or by using the simulator (sm) tool.

WARNING: The task simalma is designed to only be invoked once for a simulation setup. It always sets up a skymodel and pointings. That means that simalma is not supposed to be run multiple times for a project, unlike simobserve and simanalyze. The task simalma may ignore or overwrite the old results when it is run more than once with the same project name.

There are options in simalma to simulate observation of ACA 7-m and total power arrays, to apply thermal noise, and/or to generate images from simulated MeasurementSets. One inputs a vector of configurations, and a corresponding vector of total times to observe each component. Thermal noise is added to visibilities when pwv > 0 . The ATM atmospheric model is constructed from the characteristics of the ALMA site and a user defined Precipitable Water Vapour (pwv) value. Set pwv = 0 to omit the thermal noise. Finally, when image = True, synthesized images are generated from the simulated MeasurementSets.

Antenna Configuration

The configurations of the ALMA 12-m and 7-m arrays are defined by the antennalist parameter, which can be a vector. Each element of the vector can be either the name of an antenna configuration file or a desired resolution, e.g., ‘alma;cycle1;5arcsec’. Some examples:

  • antennalist = [‘alma.cycle2.5.cfg’,’aca.cycle2.i.cfg’]; totaltime = [‘20min’,’2h’]’: Will observe the 12-m array in configuration C32-5 for 20 minutes and the ACA 7-m array for 2 hours.

  • antennalist = [‘alma;cycle2;0.5arcsec’,’aca.i.cfg’]; totaltime = [‘20min’,’2h’]’: Will observe the 12-m array in whatever Cycle 2 configuration yields a zenith synthesized beam as close as possible to 0.5 arcsec (at the center frequency of your skymodel) for 20 minutes and the ACA 7-m array for 2 hours.

  • antennalist = [‘alma.cycle1.2.cfg’,’aca.cycle2.i.cfg’]; totaltime = ‘20min’: Will observe the 12-m array in Cycle 1 configuration 2 for 20 minutes and the ACA 7-m array for the default of 2×(12-m time) = 1h20min. This parameter setting will also generate a warning that the user is combining configurations from different ALMA Cycles (but the simulation will run despite that).

Total power can either be included along with interferometric configurations e.g., antennalist = [‘alma.cycle1.2.cfg’,’aca.cycle2.i.cfg’,’aca.tp.cfg’], or by using the tpnant and tptime parameters. The latter is preferred since it allows greater control (in particular the number of total power antennas to use – if more than one is used, multiple total power observations will be generated and combined in imaging).

Field Setup

There are two ways to setup pointings, i.e., Rectangle Setup and Multi-Pointing.In the Rectangle Setups, pointings are automatically calculated from the pointing center (direction) and the map size. A rectangular map region is covered by a hexagonal grid (maptype = ‘alma’) with Nyquist sampling, i.e., 0.48 primary beam (PB) spacing (where PB ≡ 1.2 λ / D), in both ALMA 12-m and ACA 7-m array simulations. A slightly larger area is mapped in ACA total power simulations for later combination with interferometer visibilities. The map area is extended by 1 PB in each direction and covered by a lattice grid with 0.225 PB spacing.

In Multi-Pointing, a list of pointings is defined in the direction parameter or read from a file (when setpointings = False; note that simobserve can read ALMA OT pointing files in the old and new format but the latter only when they are saved as sexagesimal absolute positions). The ALMA 12-m and ACA 7-m arrays observe the specified directions. The ACA total power simulations map either (1) square regions of 2 PB extent centered at each of the pointings, or (2) a rectangle region that covers all the pointings. Either (1) or (2), whichever can be done with the smaller number of points, is selected. The pointing spacing in total power simulations is, again, 0.225 PB in lattice grids.

It is advisable that for Total Power Simulations, the field is chosen sufficiently large, maybe padding at least 1-2 primary beams on each side.

Integration time

The total observation time of each component or configuration is defined by the totaltime parameter as noted above. A scalar will trigger use of the Cycle 2 default time multipliers, 1:0.5:2:4 for the first 12-m configuration, any additional 12-m configurations, any 7-m configuration, and any total power observation.

In general, the integration time (dump interval) of simulations is defined by the integration parameter with an exception. Since the ACA total power array always observes larger areas compared to the ALMA 12-m and ACA 7-m arrays, it is possible that the ACA total power array cannot cover all pointings in the given observation time. In such a case, the integration time in the total power simulation is scaled so that the all pointings are observed at least once in its observation time, i.e., integration_TP = tptime / (the number of total power pointings).

Imaging and combination of ALMA with ACA

The CLEAN algorithm is used in simalma to generate images from visibilities. The visibilities are weighted to UV-plane using Briggs weighting.When ACA observations are simulated, visibilities of ACA 7-m are weighted by the relative sensitivities to ALMA 12-m visibilities, and both data sets are concatenated before imaging. The relative weight of ACA 7-m visibilities is defined in proportion to the ratio of beam areas squared, i.e., \((7/12)^{4} = 0.11\). This is because simalma uses a bandwidth and an integration time common to both ALMA 12-m and ACA 7-m simulations.

The interferometer and total power images are combined using feather task when total power observations are included. The total power image is scaled by the interferometer primary beam coverage before combination. The final image product is the combined image corrected for the interferometer primary beam coverage. The output image of the feather task is divided by the interferometer primary beam coverage in the final step.


simutil

simutil contains numerous utility methods which can assist users in generic ephemeris and geodesy calculations to aid in performing simulations and other activities in CASA, as well as some methods used internally by simobserve and simanalyze. Several of these tasks directly call the simulator tool in an attempt to lessen the amount of scripting required and to make it easier for the user. It is used by import and instantiation, similarly to testhelper and recipes:

from simutil import simutil

u=simutil()
help u.readantenna

Antenna configuration files are important for several tasks in simutil and other simulator tools. Below is an example of a properly formatted configuration file.

#observatory=ALMA
#COFA=-67.75,-23.02
#coordsys=LOC (local tangent plane)
# uid___A002_Xdb6217_X55ec_target.ms
# x             y               z             diam  station  ant
-5.850273514   -125.9985379    -1.590364043   12.   A058     DA41
-19.90369337    52.82680653    -1.892119601   12.   A023     DA42
13.45860758    -5.790196849    -2.087805181   12.   A035     DA43
5.606192499     7.646657746    -2.087775605   12.   A001     DA44
24.10057423    -25.95933768    -2.08466565    12.   A036     DA45

The observatory, COFA (center of array), coordsys (coordinate system), x, y, z, diam (diameter) and name will be interpreted as header keys or value pairs if they contain “=” and begin with #, and as comments otherwise. Other possible header keys are: zone, datum, or hemisphere. If no sixth column is provided, antenna names will default to station names. If no fifth column is provided, station names will default to A0x where x is the zero-indexed row number. To find the observatory name, one can check the known observatories list by using the measures tool command me.obslist. If an unknown observatory is specified, then one either must use absolute positions (coordsys, XYZ [Cartesian coordinates], UTM [Universal Transverse Mercator]), or specify COFA (longitude and latitude). coordsys can be XYZ (Earth-centered), UTM (easting, northing, and altitude), or LOC (xoffset, yoffset, and height). Files for many observatories can be found in the directory returned by the following command:

casa.values()[0]['data']+"/alma/simmos"

Tsys and Noise

simutil.noisetemp

Noise temperature and efficiencies can be calculated for several telescopes: ALMA, ACA, EVLA, VLA, and SMA. The inputs for simutil.noisetemp method include: telescope, e.g., “ALMA”, freq (observing frequency) as a quantity string, e.g., “300GHz”, diam (optional - knows diameters for arrays above), e.g., “12m”, epsilon = RMS surface accuracy in microns (also optional - this method contains the engineering specification values for each telescope). The outputs produced \(\eta_p\) phase efficiency (from Ruze formula), \(\eta_s\) spill (main beam) efficiency, \(\eta_b\) geometrical blockage efficiency, \(\eta_t\) taper efficiency, \(\eta_q\) correlator efficiency including quantization, \(t_{rx}\) receiver temperature. Where the total antenna efficiency can be calculated from these outputs as such: \(\epsilon = \eta_p * \eta_s * \eta_b * \eta_t\).

NOTE: VLA correlator efficiency includes waveguide loss. EVLA correlator efficiency is probably optimistic at 0.88.

simutil.sensitivity

This method is used to calculate the noise in an observation by adding noise to visibilities in exactly the same way as sm.corrupt (if doimnoise=True) and also creates a simulated image from which to measure noise. The inputs to calculate sensitivity are: freq, bandwidth = channel width, e.g., “1GHz”, etime = exposure time / length of track, e.g., “500sec”, integration = scan time, e.g., “10sec”, elevation, e.g., “80deg”; either an antennalist (a simobserve-format antenna configuration filename) must be given or the parameters telescope, diam, and nant (number of antennas) must be set. Other optional inputs include: pwv in mm, doimnoise uses the simulator task sm.corrupt to create an MS and image it to measure the noise, integration or integration time (units required) e.g., “10s”, debug, method which is equivalent to the mode parameter in the simulator task sm.setnoise (options: “tsys-atm” (default) or “tsys-manual”), tau0 or the zenith atmospheric opacity (must be set if method=”tsys-manual”), and t_sky (default=200 (K) when method=”tsys-manual”).

Geodesy and Antenna Positions

NOTE: For more information on geodesy and pointing and other helper functions that are useful and available at https://www.ngs.noaa.gov/TOOLS/program_descriptions.html.

The ITRF frame mentioned in several of the following tasks is not the official ITRF (International Terrestrial Reference Frame), just a right-handed Cartesian system with X going through 0 latitude and 0 longitude, and Z going through the north pole.

simutil.readantenna

simutil.readantenna is a helper function to read antenna configuration files, using the antab parameter as an input. Outputs will be: earth-centered x,y,z, diameter, name, observatory_name, observatory_measure_dictionary.

NOTE: The observatory_measure_dictionary output was added between CASA 4.7 and 5.0.

simutil.baselineLengths

When given an antenna configfile, this method will return the zenith baseline lengths.

simutil.approxBeam

When given an antenna configfile and freq (in GHz), this method will return the approximate beam size at zenith from the 90th percentile baseline length.

simutil.long2xyz

This method returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at geodetic latitude (parameter lat) and longitude (parameter lon) [radians] and elevation [m].

simutil.xyz2long

When given ITRF Earth-centered (X, Y, Z, using the parameters x, y, and z) coordinates [m] for a point, this method returns geodetic latitude and longitude [radians] and elevation [m]. Elevation is measured relative to the closest point to the (latitude, longitude) on the WGS84 (World Geodetic System 1984) reference ellipsoid.

simutil.locxyz2itrf

This method returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at “local” (x, y, z, using the parameters locx, locy, and locz) [m] measured at geodetic latitude (lat) and longitude (longitude) [degrees] and altitude (alt) of the reference point. The “local” (x, y, z) are measured relative to the closest point on the WGS84 reference ellipsoid, with z normal to the ellipsoid and y pointing north.

simtuil.itrf2loc

Given Earth-centered ITRF (X, Y, Z, using the parameters x, y, and z) coordinates [m] and the Earth-centered coords of the center of array (using the parameters cx, cy, and cz), this method returns local (x, y, z) [m] relative to the center of the array, oriented with x and y tangent to the closest point at the COFA (latitude, longitude) on the WGS84 reference ellipsoid, with z normal to the ellipsoid and y pointing north.

simutil.itrf2locname

Given Earth-centered ITRF (X, Y, Z) coordinates [m] and the name of an known array using the obsname parameter (see me.obslist), the method simutil.itrf2locname returns local (x, y, z) [m] relative to the center of the array, oriented with x and y tangent to the closest point at the COFA (latitude, longitude) on the WGS84 reference ellipsoid, with z normal to the ellipsoid and y pointing north.

simutil.utm2xyz

This method returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at UTM easting, northing, elevation [m], and zone of a given datum (e.g., ‘WGS84’) and north/south flag nors (“N” or “S”, denotes northern or southern hemisphere). The ITRF frame used is not the official ITRF, just a right-handed Cartesian system with X going through 0 latitude and 0 longitude, and Z going through the north pole.

simutil.utm2long

The method simutil.utm2long converts UTM coordinates to GPS longitude and latitude (in radians). This task has the following parameters: east, north, zone, datum, and nors.

Pointing and Directions

simutil.calc_pointings2

This method is used to calculate mosaic pointings to cover a region. This returns a hexagonally packed list of pointings determined by the size (either [size[0],size[1]] or [size,size] if a single value is given) parameter separated by parameter spacing and fitting inside an area specified by direction and maptype. If multiple pointings can not be fit to the given parameters, a single pointing will be returned. ]If direction is a list, the task simply returns the direction and the number of pointings in it. The 3 options for maptype are: “HEX”agonal (default), “SQU”are, and “ALM”A (triangular tiling). The hexagonal packing starts with a horizontal row centered on direction, and the other rows alternate being horizontally offset by a half spacing. For hexagonal or square maptypes, the relmargin (default=0.5) parameter affects the number of pointings returned in the mosaic pattern. For triangular maptypes, the beam parameter is used to determine the number of pointings returned in the mosaic pattern, although this parameter is optional.

simutil.read_pointings

This method will read a pointing list from a file using the parameter filename. The input file (ASCII) should contain at least 3 fields separated by a space which specify positions with epoch, RA and Dec (in degrees / minutes / seconds or hours / minutes / seconds). The optional field and time columns should be a list of decimal numbers which specifies integration time at each position (in units of seconds). The lines which start with ‘#’ are ignored and can be used as comment lines. Example of a file:

#Epoch     RA          DEC      TIME(optional)
 J2000 23h59m28.10 -019d52m12.35 10.0
 J2000 23h59m32.35 -019d52m12.35 10.0
 J2000 23h59m36.61 -019d52m12.35 60.0

simutil.write_pointings

This method will write a list of pointings out to a file (example above), given by the parameter filename. The optional parameter time can be an array of integration times.

simutil.average_direction

This method will return the average of directions (default=None) as a string, and relative offsets.

simutil.median_direction

This method will return the median of directions (default=None) as a string, and relative offsets.

simutil.ephemeris

This method calculates the elevation of a source on a given date, in a given direction, seen from a given telescope. The date should be given in the format YEAR / MO / DY / TI:ME. The time given is referenced with the International Atomic Time, or TAI (from the French name name temps atomique international). Other optional parameters include: usehourangle (boolean parameter which sets or unsets the reference time at transit, essentially centering the plot), ms (uses the information from the OBSERVATION table in the given MeasurementSet and plots the entire range of the observation), and cofa (allows the user to change the center of the array position). The cofa parameter must be set if using an unknown observatory. A list of known observatories can be found by using the measures tool command me.obslist.

Utility

simutil.statim

This method will plot an image and calculate its statistics. Optional parameters: plot (default True), incell, disprange (low and high values for pl.imshow), bar (show colorbar, default=True), showstats (show stats on the image, default=True).

simutil.plotants

An alternate antenna configuration plotting routine that takes arrays of x,y=local offset from the array center, z=altitude, d=diameter, and name. This method routine either plots points or, if the array is compact enough to see the diameters, plots to the actual scaled size of the dishes.

simutil.modifymodel

simutil.modifymodel is a method that converts a model image into a 4D-coordinate image that can be used in CASA, with axes in space, stokes, spectral order, which the Toolkit requires (e.g., sm.predict in the simulator tool). The input parameters inimage and outimage allow the user to specify the names of the input and output. Values that are absent in the input, or that the user wishes to override, can be input as quantity strings with the in* parameters (inbright, indirection, incell, incenter, inwidth, innchan). e.g., inbright=”4Jy/pixel” will scale outimage to have 4Jy/pixel peak, incell=”0.2arcsec” will set the cell size in outimage to 0.2arcsec. The flatimage parameter allows one to also generate a flat (2D, integrated intensity) image from inimage, which can be useful for display purposes.

simutil.convimage

Given a (2D) model (modelflat) image, this method will regrid it to the scale of the outflat image, and convolve it to the beam of the outflat image. This is useful to compare a skymodel with a simulated output image. The optional parameter complist allows the user to import a componentlist to add unresolved components to the outflat image. Information on creating a component list can be found in the CASA guides here.

simutil.imtclean

This wrapper function is the method by which the standard CASA imaging task tclean is called for simulated image reconstruction inside the task simanalyze[. It replaces the deprecated method simutil.imtclean. If dryrun=True, this method only creates a template ‘[imagename.config].tclean.last’ file for users to reference in their custom calls to tclean. The cell parameter expects a list of qa.quantity objects. ]Selecting individual fields for imaging is not supported.