Moment maps - eakadams/ihow-hi GitHub Wiki

After creating a clean cube, the next step is to create moments based on the cube, which provide 2D views of the data and are useful in the data analysis and interpretation stages.

Moment zero map

The first map we will make is the moment zero map. This is a total intensity image of the HI. Here we will make the moment zero map two different ways so that you can see the impact of masking on the moment map.

First, create a moment zero map from the cube, over the channels with signal (previously identified as channels 7~66) without any additional masking:

immoments(imagename='sexb_cube.image',moments=[0],chans='7~66',outfile='mom0_nomask.im')

Note that there will be warnings as this runs; the beam size varies slightly from channel to channel so the task convolves to a common beam and warns about that because the differences are small, requiring a small kernel for convolving. You could alternativey use the "common resolution" cube here that we created at the end of the previous step. This is your choice to explore and see.

Also create a moment zero map masking with the clean mask. In this case, because we are using a mask, we don't have to worry about the channel range selection:

immoments(imagename='sexb_cube.image',moments=[0],mask='sexb_cube.mask',outfile='mom0_cleanmask.im')

Examine these images in the viewer. What do you notice that is similar or different between them? Below is a side-by-side display of the two moment zero maps: Note that you can fix the intensity display in CASA to be global, which is a nice way to compare.

Depending on your science goals and your source, you may want to apply different types of masking to create your moment zero map. For example, when looking for low surface brightness emission, not masking is often the best way to pull that emission out. But when studying bright emission, masking can help increase your S/N. And by restricting to cleaned regions, you can also be sure that all your emission is in the same beam units.

In addition, it could be that you want to create a mask that is more specific than the clean mask, e.g., places more stringent limits. You generally don't want to use a mask that is more generous than your clean mask because then you have uncleaned emission included.

Converting to column density

The units of a moment zero total intensity map are Jy/beam km/s. These are not the most intuitive units, and it's often best to convert to column density, which measures neutral hydrogen atoms per cm^2. You can find equations for how to do this at this helpful webpage.

First, you will need to know the beamsize, so find that for the moment zero map. We will go ahead and use the masked moment zero map for now, but you could also use the unmasked map.

imhead(imagename='mom0_cleanmask.im')

The restoring beam is 45.90"x43.05".

Now we can do the conversion to column density:

immath(imagename='mom0_cleanmask.im',mode='evalexpr',expr='IM0*606*1.823e21/(45.90*43.05)',outfile='nhi_cleanmask.im')

The image below shows the NHI map with contours overlaid at [1, 2, 4, 8, 12]x10^20 atoms cm^-2:

Primary beam correction

Note that our telescope does not have the same sensitivity across the full field of view; the primary beam pattern of the individual dishes impacts our measured response. Thus, we need to corret for this before deriving information such as fluxes from out data when our sources are extended. This can be done in the tclean task itself, but it also changes the noise profile of the cube, as noise is elevated further from the center of the image. Thus, I often prefer to manually run the primary beam correction at the stage I am ready for it, which is for moment-zero maps. CASA provides a task for this, but note that it can also be done yourself in simple cases - the primary beam correction is simply dividing your image by the primary beam response.

Our primary beam image is a cube, and since we want to correct an image, we need to extract a single channel of this data cube. We will pick a central channel, near the center of our source emission. The primary beam response does change with frequency, but is neglible for the narrow frequency range we are working with.

immath(imagename='sexb_cube.pb',chans='50',expr='IM0',outfile='pb_onechan.im')

Now we can do the primary beam correction. Since this does elevate the noise, we will only do it to already masked images.

impbcor(imagename='mom0_cleanmask.im',pbimage='pb_onechan.im',outfile='mom0_cleanmask_pbcorr.im')

impbcor(imagename='nhi_cleanmask.im',pbimage='pb_onechan.im',outfile='nhi_cleanmask_pbcorr.im')

You can plot these in viewer and compare contours between the primary beam corrected and non-corrected to see what the differences are. To further understand the (lack of) differences you can also plot contours of the primary beam image to see what level of primary beam response you are at. What level of the primary beam response is the galaxy contained within?

Moment 1 and 2 maps

After the total intensity map, the moment one and moment two maps are then the most useful. The moment one map is the intensity weighted by the velocity, and thus represents the velocity field of the galaxy (e.g., shows the recessional velocity for each pixel). The second moment map is weighted by the velocity squared and represents the dispersion of the gas, if there is a single spectral component. (If there are two or more spectral components, than the second moment is impacted by this.) These moments are made analagously to the first moment, simply specifying a different order. However, these moments are much more sensitive to noise and thus it is advisable to clip them at a higher level. A good practice is to use the moment zero map as a mask - pick a threshold/significance level of that moment zero map and use it to the mask the higher order moments. In this case, we will use 1.8x10^20, which is ~4-sigma (based on noise in the unmasked moment zero map - an advantage of that map).

First, get the moment one and two maps, applying the clean mask. While we could do this in one call, specifying more than one moment order, than CASA appiles some automatic naming. Instead, we will run separate calls.

immoments(imagename='sexb_cube.image',moments=[1],mask='sexb_cube.mask',outfile='mom1_cleanmask.im')
immoments(imagename='sexb_cube.image',moments=[2],mask='sexb_cube.mask',outfile='mom2_cleanmask.im')

Then mask the moment one and two based on the N_HI map:

immath(imagename='mom1_cleanmask.im',mode='evalexpr',expr='IM0',mask='nhi_cleanmask.im>=1.8e20',outfile='mom1_nhimask.im')
immath(imagename='mom2_cleanmask.im',mode='evalexpr',expr='IM0',mask='nhi_cleanmask.im>=1.8e20',outfile='mom2_nhimask.im')

Explore the moment maps in the viewer. Can you see the impact of adding the moment zero based clipping?

The moment one and two maps:

Write moment maps to fits

By writing your maps to fits, you can then use them in other programs. For example, you can load into astropy and plot HI contours on an optical image. Or you can make more tailored publication quality images by using python.

exportfits(imagename='nhi_cleanmask.im',fitsimage='sexb_nhi.fits',velocity=True,optical=True,dropdeg=True)
exportfits(imagename='mom1_nhimask.im',fitsimage='sexb_mom1.fits',velocity=True,optical=True,dropdeg=True)
exportfits(imagename='mom2_nhimask.im',fitsimage='sexb_mom2.fits',velocity=True,optical=True,dropdeg=True)

Note that specifying dropdeg=True is very helpful as it makes sure you have a 2D image rather than 3- or 4-D, even though the Stokes and frequency axes are empty, which can cause problems when using tools designed around 2D optical images (e.g., astropy). In this case, we also specify optical=True since we don't have a velocity axis anymore and want to generally try and ensure the best comparison to multi-wavelength, which is generally optical/wavelength-based.

Get the HI line flux and calculate the HI mass

You can use imstat on the masked, primary-beam corrected moment zero map to calculate the total line flux (within the clean mask)

imstat(imagename='mom0_cleanmask_pbcorr.im') 

The output you should see is something like:

2023-09-07 11:01:36 INFO imstat	         -- flux [flux]:                            79.6524 Jy.km/s

You can then check the VLA ANGST paper to see how your line flux compares to theirs. Is it the same? How different is it? Can you think of reasons why they might be different?

Their value is 91.0 Jy km/s, while the reported single-dish is 72.9 Jy km/s. Thus, we are ~12% lower in line flux. Remember that ~10% is commonly taken as the level of flux accuracy with the VLA, so we are formally consistent. In addition, masking adds an additional level of uncertainty as how you select the mask can impact the final line flux. You could experiment with using a different mask, for example, to see how that changes the line flux you recover.

With an HI line flux and a distance, you can always directly compute an HI mass (assuming the emission is optically thin, which is generally true). This is much more accurate than stellar masses, for example, where the assumptions and models of the stellar population need to be used. The formula for calculating the HI mass is:

M_HI (Msun) = 2.365e5 * D(Mpc)^2 * S (Jy km/s)

Taking the distance of 1.39 Mpc for Sextans B from the VLA ANGST paper, what do you find for the HI mass?

You should find an answer around ~3.5x10^7 Msun.

Explore overlays

If you have time, you can go further by exploring how to overlay HI contours on an optical image. There is a small notebook for doing this in the notebook directory of this repository.