Sunspot Band6 Feathering for CASA 6.5.4: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Line 15: Line 15:


==Prepare the Single-Dish Data for Combination==
==Prepare the Single-Dish Data for Combination==
We start the process on CASA by defining some names of the data.
We by defining names for our images.


<source lang='python'>
<source lang='python'>
Line 27: Line 27:
[[File:sunspot_sdspot.png|thumb|right|'''Fig. 2.''' The trimmed map from the full-sun map.]]
[[File:sunspot_sdspot.png|thumb|right|'''Fig. 2.''' The trimmed map from the full-sun map.]]


We shift the coordinate of the full-sun map created from the single-dish data, based on the RA/Dec coordinate obtained in the previous section, as follows.
We shift the coordinate of the full-sun map created from the single-dish data, based on the RA/Dec coordinate obtained in the previous section, with {{imhead_6.5.4}} as follows.


<source lang='python'>
<source lang='python'>
Line 38: Line 38:
</source>
</source>


The unit of the full-sun map is Kelvin. On the other hand, the unit of the synthesized image is Jy/beam. So, the conversion of the unit has to be done. In the tutorial, we convert the unit of the full-sun image to Jy/beam, using the following commands.
The unit of the full-sun map is Kelvin. On the other hand, the unit of the synthesized image is Jy/beam. We convert the unit of the full-sun image to Jy/beam, using the following commands with {{immath_6.5.4}} and {{imhead_6.5.4}}.


<source lang='python'>
<source lang='python'>
# In CASA
# In CASA
# derive a scaling factor based on frequency and beam size
crval4S= imhead(sd_img_base+'.rescl', mode='get',hdkey='crval4')
crval4S= imhead(sd_img_base+'.rescl', mode='get',hdkey='crval4')
freqSD = crval4S['value']
freqSD = crval4S['value']
sbema = imhead(sd_img_base+'.rescl', mode='get',hdkey='beammajor')
sbema = imhead(sd_img_base+'.rescl', mode='get', hdkey='beammajor')
sbemi = imhead(sd_img_base+'.rescl', mode='get',hdkey='beamminor')
sbemi = imhead(sd_img_base+'.rescl', mode='get', hdkey='beamminor')
sbmajor = sbema['value']
sbmajor = sbema['value']
sbminor = sbemi['value']
sbminor = sbemi['value']
Line 53: Line 55:


imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='crval3', hdvalue=['I'])
imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='crval3', hdvalue=['I'])
imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='bunit', hdvalue='Jy/beam')
imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='bunit', hdvalue='Jy/beam')
</source>
</source>

Revision as of 16:14, 12 April 2024

Last checked on CASA Version 6.5.4

Overview

This guide features CARTA, the “Cube Analysis and Rendering Tool for Astronomy,” which is the new NRAO visualization tool for images and cubes. The CASA viewer (imview) has not been maintained for a few years and will be removed from future versions of CASA. We strongly recommend using CARTA, as it provides a much more efficient, stable, and feature rich user experience. A comparison of the CASA viewer and CARTA, as well as instructions on how to use CARTA at NRAO, is provided in the CARTA section of the CASA docs.

This portion of the Sunspot Band6 guide covers the combination process between the synthesized solar image and the full-sun image obtained with the TP array. It is the final step after: 1) calibrating the IF visibility data; 2) imaging the IF data; 3) calibrating and imaging the single-dish data. If you wish, you may skip parts 1 and 2, download the calibrated data package, Sunspot_Band6_CalibratedData.tgz, and complete part 3.

If you completed parts 1 and 2, we assume that you are working in the directory Sunspot_Band6_UncalibratedData. If you skip parts 1 and 2, and use the provided calibrated data, we assume that the working directory is Sunspot_Band6_CalibratedData.

Obtain the RA/Dec Coordinate of the Sun Center at the Reference Time

Before starting the combination process in CASA, we need to know the RA/Dec coordinate of the sun center at the reference time, which is used during the interferometric observation, because the observing time of the TP-array is not exactly the same. The Sun was moving on the RA/Dec coordinate frame during the gap between these observations. We need to shift the coordinate of the full-sun image based on the coordinate of the sun center at the reference time. The reference time of the observation is:
19:49:00UT, December 18 2015.

We obtain the RA/Dec coordinate of the sun center estimated from JPL's HORIZONS system:
RA: 266.01898 degree, Dec: -23.38648 degrees.

Prepare the Single-Dish Data for Combination

We by defining names for our images.

# In CASA
int_base = 'AR12470_B6AllSpw_I'                         # basename of the interferometric image
sd_img_base = 'uid___A002_Xae00c5_X2e6b.PM03.StkI.Spw3' # basename of the single dish image
out_base = 'AR12470_B6AllSpw_I_Feather'                 # basename of the combined products
Fig. 1. The full-sun map with 'Jy/beam' unit.
Fig. 2. The trimmed map from the full-sun map.

We shift the coordinate of the full-sun map created from the single-dish data, based on the RA/Dec coordinate obtained in the previous section, with imhead as follows.

# In CASA
import numpy as np
sun_ra = str(266.01898 * np.pi/180.)
sun_dec = str(-23.38648 * np.pi/180.)
imhead(imagename=sd_img_base+'.rescl', mode='put', hdkey='crval1', hdvalue=sun_ra)
imhead(imagename=sd_img_base+'.rescl', mode='put', hdkey='crval2', hdvalue=sun_dec)

The unit of the full-sun map is Kelvin. On the other hand, the unit of the synthesized image is Jy/beam. We convert the unit of the full-sun image to Jy/beam, using the following commands with immath and imhead.

# In CASA

# derive a scaling factor based on frequency and beam size
crval4S= imhead(sd_img_base+'.rescl', mode='get',hdkey='crval4')
freqSD = crval4S['value']
sbema = imhead(sd_img_base+'.rescl', mode='get', hdkey='beammajor')
sbemi = imhead(sd_img_base+'.rescl', mode='get', hdkey='beamminor')
sbmajor = sbema['value']
sbminor = sbemi['value']
convKJb = str(13.6 * (300./(freqSD/1.e9))**2. * (1./sbmajor)*(1./sbminor))

immath(imagename=sd_img_base+'.rescl', expr='IM0/'+convKJb, outfile=sd_img_base+'.jyb.im')

imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='crval3', hdvalue=['I'])
imhead(imagename=sd_img_base+'.jyb.im', mode='put', hdkey='bunit', hdvalue='Jy/beam')

The result of the process is shown in Figure 1. A look of the image is not changed, but you can find the changing of the values and unit in the 'Sunspot_Band6_SingleDish_for_CASA_4.7' tutorial.

Then we do the re-gridding and trimming to fit the map area to that of the synthesized image. The area of the trimming is set to the same as the CLEAN box, in the tutorial.

# In CASA
#re-griding of SD data
imregrid(imagename=sd_img_base+'.jyb.im', template=int_base+'.image', axes=[0, 1], output=sd_img_base+'.jyb.regrid')

# Trimming to exclude noisy edge regions in the TP image
imsubimage(imagename=sd_img_base+'.jyb.regrid', outfile=sd_img_base+'.jyb.regrid.subim', box='100,100,400,400')

As the final process of the preparation of the single-dish data, we multiply the trimmed image by the primary beam response, as follows.

# In CASA
imsubimage(imagename=int_base+'.pb',outfile=int_base+'.flux.subim',box='100,100,400,400')

immath(imagename=[sd_img_base+'.jyb.regrid.subim',int_base+'.flux.subim'], expr='IM0*IM1', outfile=sd_img_base+'.jyb.regrid.subim.depb')

Prepare the Synthesized Image for Combination

Fig. 3 The trimmed map from the synthesized image.

To fit the map area to that of the trimmed single-dish image, we do the trimming for the synthesized image.

# In CASA
imsubimage(imagename=int_base+'.image', outfile=int_base+'.image.subim', box='100,100,400,400')

Combination and Primary Beam Correction

The combining of the maps can be done using 'feathering' task, as the same as the non-solar data.

# In CASA
feather(imagename=out_base+'.image', highres=int_base+'.image.subim',lowres=sd_img_base+'.jyb.regrid.subim.depb')

CAUTION: In this tutorial and the packages of the solar SV data released on 2017/01/18, we use the default parameters of the task, as shown above. As a result, the averaged brightness temperature of the combined image is always larger a few % than the temperature brightness at the same position in the single-dish map, although the values have to be the same basically. It means that you need to tune the “sdfactor" parameter of the task for obtaining consistent images, before using combined images for your science. The detail of the parameter is described in Section 5.6 of “CASA User Reference & Cookbook Release 4.7.0”[1].

After the feathering process, we do the primary beam correction, as follows.

# In CASA
impbcor(imagename=out_base+'.image', pbimage=int_base+'.flux.subim', outfile=out_base+'.pbcor', mode='divide')

To convert the unit from Jy/beam to Kelvin, we use the following commands.

Fig. 4. The trimmed map from the synthesized image.
# In CASA
crval4I= imhead(int_base+'.image', mode='get',hdkey='crval4')
freq = crval4I['value']
fbema = imhead(out_base+'.image', mode='get',hdkey='beammajor')
fbemi = imhead(out_base+'.image', mode='get',hdkey='beamminor')
fbmajor = fbema['value']
fbminor = fbemi['value']

convKJb2 = str(13.6 * (300./(freq/1.e9))**2. * (1./fbmajor)*(1./fbminor))

immath(imagename=out_base+'.pbcor',expr='IM0*'+convKJb2, outfile= out_base+'_K.pbcor')

imhead(imagename=out_base+'_K.pbcor', mode='put', hdkey='bunit', hdvalue='K')

Finally, we create the FITS file of the combined and primary corrected map.

# In CASA
exportfits(imagename=out_base+'_K.pbcor', fitsimage=out_base+'_K.fits')

The result of the combining is shown in Figure 4.


Caution: The single-dish images and feathering images in the solar SV data packages (ex. Sunspot_Band6_ReferenceImages.tgz) were not applied the re-scaling to the single-dish data (see Sunspot_Band6_SingleDish_for_CASA_6.5.4#Calibration_.28re-scaling.29_after_Imaging). Please do not use the reference image for your science.