|
|
Line 38: |
Line 38: |
| =Plot light-curve and spectrum= | | =Plot light-curve and spectrum= |
| =Quick-Look Imaging= | | =Quick-Look Imaging= |
| =Cross-Power Dynamic Spectrum=
| |
| The first module we introduce is [https://github.com/suncasa/suncasa/blob/master/utils/dspec.py ''dspec'']. This module allows you to generate a cross-power dynamic spectrum from an MS file, and visualize it. You can select a subset of data by specifying a [https://casa.nrao.edu/Release3.3.0/docs/UserMan/UserMansu112.html time range], [https://casaguides.nrao.edu/index.php/Selecting_Spectral_Windows_and_Channels spectral windows/channels], [https://casaguides.nrao.edu/index.php/Antenna/Baseline_Selection_Syntax_with_or_without_Autocorrelations antenna baseline], or [https://casa.nrao.edu/Release3.3.0/docs/UserMan/UserMansu113.html uvrange]. The selection syntax follows the ''CASA'' convention. More information of CASA selection syntax may be found in the above links or the [https://casa.nrao.edu/casadocs/casa-5.4.0/data-selection/data-selection-in-a-measurementset Measurement Selection Syntax].
| |
| [[file:fig-dspec.png|thumb|right|300px|Figure 1: EOVSA cross power dynamic spectrum at stokes XX and YY]]
| |
| <pre style="background-color: #FCEBD9;">
| |
| from suncasa.utils import dspec as ds
| |
| import matplotlib.pyplot as plt
| |
| ## define the visbility data file
| |
| msfile = 'IDB20170821201800-202300.4s.slfcaled.ms'
| |
| ## define the output filename of the dynamic spectrum
| |
| specfile = msfile + '.dspec.npz'
| |
| ## select relatively short baselines within a length (here I use 0.15~0.5km),
| |
| ## and take a median cross all of them (with the domedian parameter)
| |
| ## alternatively, you can use the "bl" parameter to select individual baseline(s)
| |
| uvrange = '0.15~0.5km'
| |
| ## this step generates a dynamic spectrum and saves it to specfile
| |
| dspec=ds.get_dspec(vis=msfile, specfile=specfile, uvrange=uvrange, domedian=True)
| |
| ## dspec is a Python dictionary that contains the resulting dynamic spectrum.
| |
| ## A copy is saved in "specfile" as a numpy npz file.
| |
| ## Other optional parameters are available for more selection criteria
| |
| ## such as frequency range ("spw"), and time range ("timeran")
| |
| ## Use "ds.get_dspec?" to see more options
| |
|
| |
| ## define the polarizations to show (here I use XX and YY)
| |
| pol='XXYY'
| |
| ## The following command displays the resulting cross-power dynamic spectrum
| |
| ds.plt_dspec(specfile, pol=pol) # alternatively you can use the output from the previous step "dspec" as the input
| |
| </pre>
| |
| Now you should have a popup window showing the dynamic spectrum. Hover your mouse over the dynamic spectrum, you can read the time and frequency information at the bottom of the window.
| |
|
| |
| ===Quick-Look Imaging===
| |
| Imaging EOVSA data involves image cleaning, as well as solar coordinate transformation and image registration. We bundled a number of these steps ino a module named [https://github.com/suncasa/suncasa/blob/master/utils/qlookplot.py ''qlookplot''], allowing users to generate an observing summary plot showing the cross power dynamic spectrum, GOES light curves and EOVSA quick-look images. Now let us start with making a summary plot of EOVSA image at the spectral window 5 (5.4 GHz).
| |
| [[file:fig-qlookplot0.png|thumb|right|300px|EOVSA full-Sun single-band quicklook image]]
| |
| <pre style="background-color: #FCEBD9">
| |
| ## in SunCASA
| |
| from suncasa.utils import qlookplot as ql
| |
| msfile = 'IDB20170821201800-202300.4s.slfcaled.ms'
| |
| ## (Optional) Supply the npz file of the dynamic spectrum from previous step.
| |
| ## If not provided, the program will generate a new one from the visibility.
| |
| specfile = msfile + '.dspec.npz'
| |
| ## set the time interval
| |
| timerange = '20:21:10~20:21:18'
| |
| ## select frequency range from 5 GHz to 6 GHz
| |
| spw = '5~6GHz'
| |
| ## select stokes XX
| |
| stokes = 'XX'
| |
| ## turn off AIA image plotting, default is True
| |
| plotaia = False
| |
|
| |
| ql.qlookplot(vis=msfile, specfile=specfile, timerange=timerange, spw=spw, \
| |
| stokes=stokes, plotaia=plotaia)
| |
| </pre>
| |
|
| |
| The empty panels are there as we only selected one polarization for imaging/spectroscopy. Feel free to explore by adjusting the above parameters, e.g. use a different time range ("timerange" parameter) and/or frequency range ("spw" parameter).
| |
|
| |
| With [https://github.com/suncasa/suncasa/blob/master/utils/qlookplot.py ''qlookplot''], it is easy to engage solar data from SDO/AIA in the summary plot. Setting ''plotaia=True'' in the [https://github.com/suncasa/suncasa/blob/master/utils/qlookplot.py ''qlookplot''] command will download SDO/AIA data at the given time to current directory and add it to the summary plot.
| |
| [[file:fig-qlookplot1.png|thumb|right|300px|EOVSA full-Sun single-band quicklook image with SDO/AIA as background]]
| |
| <pre style="background-color: #FCEBD9">
| |
| ## in SunCASA
| |
| ql.qlookplot(vis=msfile, specfile=specfile, timerange=timerange, spw=spw, \
| |
| stokes=stokes, plotaia=True)
| |
| </pre>
| |
| The resulting radio image is a 4-D datacube (in solar X-pos, Y-pos, frequency, and polarization), which is, by default, saved as a fits file ''msfile + '.outim.image.fits''' under your working directory. The name of the output fits file can be specified using the "outfits" parameter. In this example, we combine all selected frequencies (specified in keyword "spw") into one image (i.e., multi-frequency synthesis). Therefore the third axis only has one plane. The fourth axis contains polarization. In this example, however, only the first one is populated ("XX"). Here is the relevant information in the header of the resulting FITS file.
| |
| <pre>
| |
| NAXIS = 4
| |
| NAXIS1 = 512/ Nx
| |
| NAXIS2 = 512/ Ny
| |
| NAXIS3 = 1/ number of frequency
| |
| NAXIS4 = 2/ number of polarization
| |
| </pre>
| |
| By default, [https://github.com/suncasa/suncasa/blob/master/utils/qlookplot.py ''qlookplot''] produces a full sun radio image (512x512 with a pixel size of 5"). If you know where the radio source is (e.g., from the previous full-Sun imaging), you can make a partial solar image around the source by specifying the image center ("xycen"), pixel scale ("cell"), and image field of view ("fov"). Here we show an example that images a 8-s interval around 20:21:14 UT using multi-frequency synthesis in 12-14 GHz and a smaller restoring beam. The microwave source is show to bifurcate into two components, which correspond pretty well with the double flare ribbons in SDO/AIA.
| |
| [[file:fig-qlookplot3.png|thumb|right|300px|EOVSA multi-frequency synthesis quicklook image]]
| |
| <pre style="background-color: #FCEBD9">
| |
| ## in SunCASA
| |
| xycen = [375, 45] ## image center for clean in solar X-Y in arcsec
| |
| cell=['2.0arcsec'] ## pixel scale
| |
| imsize=[128] ## number of pixels in X and Y. If only one value is provided, NX = NY
| |
| fov = [100,100] ## field of view of the zoomed-in panels in unit of arcsec
| |
| ql.qlookplot(vis=msfile, specfile=specfile, timerange='20:21:10~20:21:18', \
| |
| spw='12GHz~14GHz', stokes='XX', restoringbeam=['6arcsec'],\
| |
| imsize=imsize,cell=cell,xycen=xycen,fov=fov,calpha=1.0)
| |
| </pre>
| |
|
| |
| Next, we will make images for every single spectral window in this data set (from spw 1 to 30, spw 0 is in the visibility data, but is not calibrated for this event).
| |
| [[file:fig-qlookplot_mbds.png|thumb|right|300px|EOVSA multi-band quicklook image]]
| |
| <pre style="background-color: #FCEBD9">
| |
| ## in SunCASA
| |
| xycen = [375, 45] ## image center for clean in solar X-Y in arcsec
| |
| cell=['2.0arcsec'] ## pixel size
| |
| imsize=[128] ## x and y image size in pixels.
| |
| fov = [100,100] ## field of view of the zoomed-in panels in unit of arcsec
| |
| spw = ['{}'.format(s) for s in range(1,31)]
| |
| clevels = [0.5, 1.0] ## contour levels to fill in between.
| |
| calpha=0.35 ## now tune down the alpha
| |
| restoringbeam=['6arcsec']
| |
| ql.qlookplot(vis=msfile, specfile=specfile, timerange=timerange, spw=spw, stokes=stokes, \
| |
| restoringbeam=restoringbeam,imsize=imsize,cell=cell, \
| |
| xycen=xycen,fov=fov,clevels=clevels,calpha=calpha)
| |
| </pre>
| |
| The output FITS file as a 4-d cube is saved to ''msfile + '.outim.image.fits''' under your working directory. The 30 spectral windows used for spectral imaging are combined in the output FITS file. So NAXIS3 (frequency axis) is 30.
| |
| <pre>
| |
| NAXIS = 4 / number of array dimensions
| |
| NAXIS1 = 128
| |
| NAXIS2 = 128
| |
| NAXIS3 = 30
| |
| NAXIS4 = 2
| |
| </pre>
| |
|
| |
|
| |
| ===Quick-Look Imaging series===
| |
| <pre style="background-color: #FCEBD9;">
| |
| # In SunCASA
| |
| msfile = 'IDB20170821201800-202300.4s.slfcaled.ms'
| |
| specfile = msfile + '.dspec.npz'
| |
| ## set the time interval
| |
| timerange = '2017/08/21/20:21:10~2017/08/21/20:21:18'
| |
| ## Bin width for time averaging
| |
| twidth = 1
| |
| ## frequency range
| |
| spw = ['0~3']
| |
| ## image center for clean in solar X-Y in arcsec
| |
| xycen = [375, 45]
| |
| ## number of pixels in X and Y. If only one value is provided, NX = NY
| |
| imsize = [128]
| |
| ## field of view of the zoomed-in panels in unit of arcsec
| |
| fov = [100., 100.]
| |
| ## pixel scale
| |
| cell = ['2.0arcsec']
| |
| ## select stokes XX
| |
| stokes = 'XX'
| |
| ## for EOVSA data, set usemsphacenter to False
| |
| usemsphacenter = False
| |
| ## set True if make a movie
| |
| mkmovie = True
| |
| ## set True if generate compressed fits
| |
| docompress = True
| |
| #### ---- Control knobs for AIA plotting ---- ####
| |
| ## set True if plot AIA images as the background
| |
| plotaia = True
| |
| ## provide the path to the directory where the AIA fits files are located. Otherwise, set it to be ''
| |
| aiadir = 'Directory-where-AIA-fits-files-are-located'
| |
| ## AIA passband. The options are [171,131,304,335,211,193,94,1600,1700]
| |
| aiawave = 171
| |
| ## numbers of CPU threads for computing
| |
| ncpu = 2
| |
|
| |
|
| |
| ql.qlookplot(vis=msfile, specfile=specfile, timerange=timerange, spw=spw,
| |
| xycen=xycen, imsize=imsize, fov=fov, cell=cell, usemsphacenter=usemsphacenter,
| |
| plotaia=plotaia, aiadir=aiadir, aiawave=aiawave,
| |
| mkmovie=mkmovie, twidth=twidth, ncpu=ncpu, docompress=docompress, stokes=stokes)
| |
|
| |
| </pre>
| |
|
| |
| ===Batch-Mode Imaging===
| |
| This section is for interested users who wish to generate FITS files with full control on all parameters being used for synthesis imaging. We provide one example SunCASA script for generating 30-band spectral imaging maps, and another for iterating over time to produce a time series of these maps.
| |
| ====Producing a 30-band map cube for a given time====
| |
| An example script can be found at [https://github.com/binchensun/eovsa-tutorial/blob/master/rhessi18/imaging_example.py this Github link]. If you are on the AWS server Virgo, it is under /common/data/eovsa_tutorial/imaging_example.py. First, download or copy the script to your own working directory and cd to your directory.
| |
| <pre style="background-color:#FCEBD9;">
| |
| # in SunCASA
| |
| CASA <##>: cd your_working_directory
| |
| CASA <##>: !cp /common/data/eovsa_tutorial/imaging_example.py ./
| |
| </pre>
| |
|
| |
| [[file:fig-specimg.png|thumb|right|300px|Example multi-frequency images at a single time integration]]
| |
|
| |
| Second (optional), change inputs in the following block in your copy of the "imaging_example.py" script and save the changes. This block has definitions for time range, image center and FOV, antennas used, cell size, number of pixels, etc.
| |
| <pre>
| |
| ################### USER INPUT GOES IN THIS BLOK ################
| |
| vis = 'IDB20170821201800-202300.4s.slfcaled.ms' # input visibility data
| |
| trange = '2017/08/21/20:21:00~2017/08/21/20:21:30' # time range for imaging
| |
| xycen = [380., 50.] # center of the output map (in solar X and Y, Unit: arcsec)
| |
| xran = [340., 420.] # plot range in solar X. Unit: arcsec
| |
| yran = [10., 90.] # plot range in solar Y. Unit: arcsec
| |
| antennas = '' # Use all 13 EOVSA 2-m antennas by default
| |
| npix = 256 # number of pixels in the image
| |
| cell = '1arcsec' # pixel scale in arcsec
| |
| pol = 'XX' # polarization to image, use XX for now
| |
| pbcor = True # correct for primary beam response?
| |
| outdir = './images' # Specify where to save the output fits files
| |
| outimgpre = 'EO' # Something to add to the output image name
| |
| </pre>
| |
|
| |
| Then, run the script in SunCASA via
| |
| <pre style="background-color:#FCEBD9;">
| |
| CASA <##>: execfile('imaging_example.py')
| |
| </pre>
| |
|
| |
| The output is a combined 30-band FITS file saved under "outdir". The naming conversion is outimgpre + YYYYMMDDTHHMMSS_ALLBD.fits. In this example, it is "./images/EO20170821T202115.000_ALLBD.fits". Here is the detailed information of the axes of the output FITS file.
| |
| <pre>
| |
| NAXIS = 4 / number of array dimensions
| |
| NAXIS1 = 128
| |
| NAXIS2 = 128
| |
| NAXIS3 = 30
| |
| NAXIS4 = 2
| |
| </pre>
| |
|
| |
| From this point, you can use your favorite language to read the fits files and plot them. The last block of the example script uses SunPy.map to generate the plot shown on the right. For SSWIDL users, please refer to [[#Plotting Images in SSWIDL|Section 3.3.3.3]].
| |
|
| |
| ====Producing a time series of 30-band maps (a 4-d cube)====
| |
| An example script can be found at [https://github.com/binchensun/eovsa-tutorial/blob/master/rhessi18/imaging_timeseries_example.py this Github link]. If you are on the AWS server Virgo, it is under /common/data/eovsa_tutorial/imaging_timeseries_example.py. The script enables parallel-processing (based on a customized task "ptclean3" -- do not ask me why we have "3" in the task name). To invoke more CPU processes for parallel-processing, change "nthreads" in the preambles from 1 to the number of threads you wish to use.
| |
|
| |
| <span style="color:#ff0000">Caution: This script is very computationally intensive. '''Please do not try to run this script on the AWS server during the RHESSI 18 live tutorial''', as it has limited resources that everyone must share. </span>
| |
|
| |
| <pre style="background-color:#FCEBD9;">
| |
| # in SunCASA
| |
| CASA <##>: cd your_working_directory
| |
| CASA <##>: !cp /common/data/eovsa_tutorial/imaging_timeseries_example.py ./
| |
| </pre>
| |
| Similar as the previous script, change inputs in the following block of your copy of the "imaging_timeseries_example.py" script and save the changes.
| |
|
| |
| [[file:fig-specmovie.png|thumb|right|300px|Example multi-frequency images at one time frame in the [https://web.njit.edu/~sjyu/download/eovsa-tutorial/movie.html output javascript movie]]]
| |
| <pre>
| |
| ################### USER INPUT GOES IN THIS BLOK ####################
| |
| vis = 'IDB20170821201800-202300.4s.slfcaled.ms' # input visibility data
| |
| specfile = vis + '.dspec.npz' ## input dynamic spectrum
| |
| nthreads = 1 # Number of processing threads to use
| |
| overwrite = True # whether to overwrite the existed fits files.
| |
| trange = '' # time range for imaging, default to all times in the data
| |
| twidth = 1 # make one image out of every 2 time integrations
| |
| xycen = [380., 50.] # center of the output map in solar X and Y. Unit: arcsec
| |
| xran = [340., 420.] # plot range in solar X. Unit: arcsec
| |
| yran = [10., 90.] # plot range in solar Y. Unit: arcsec
| |
| antennas = '' # default is to use all 13 EOVSA 2-m antennas.
| |
| npix = 128 # number of pixels in the image
| |
| cell = '2arcsec' # pixel scale in arcsec
| |
| pol = 'XX' # polarization to image, use XX for now
| |
| pbcor = True # correct for primary beam response?
| |
| grid_spacing = 5. * u.deg # grid spacing in degrees
| |
| outdir = './image_series/' # Specify where to save the output fits files
| |
| imresfile = 'imres.npz' # File to write summary of the imaging results
| |
| outimgpre = 'EO' # Something to add to the image name
| |
| </pre>
| |
|
| |
| Then, run the script in SunCASA by
| |
| <pre style="background-color:#FCEBD9;">
| |
| CASA <##>: execfile('imaging_timeseries_example.py')
| |
| </pre>
| |
|
| |
| The output fits images and a summary file imresfile are saved under "outdir" (in this example "./image_timeseries"). The naming convention of output fits images is the same as [[#Producing a 30-band map at a given time|the previous section]]. The summary yields the time, frequency, and path to every image.
| |
| <pre style="background-color:#FCEBD9;">
| |
| CASA <##>: imres = np.load(outdir + imresfile)['imres'].item()
| |
| CASA <##>: imres.keys()
| |
| ['FreqGHz', 'ImageName', 'Succeeded', 'BeginTime', 'EndTime']
| |
| CASA <##>: ls image_series/*.fits | head -5
| |
| image_series/EO20170821T201800.500_ALLBD.fits
| |
| image_series/EO20170821T201804.500_ALLBD.fits
| |
| image_series/EO20170821T201808.500_ALLBD.fits
| |
| image_series/EO20170821T201812.500_ALLBD.fits
| |
| image_series/EO20170821T201816.500_ALLBD.fits
| |
| </pre>
| |
| The last block of the example script uses SunPy.map to generate a series of plots and wraps them as a [https://web.njit.edu/~sjyu/download/eovsa-tutorial/movie.html javascript movie].
| |
|
| |
| ====Plotting Images in SSWIDL====
| |
| All the output files are in standard FITS format in the Helioprojective Cartesian coordinate system (that most spacecraft solar image data adopt; FITS header CTYPE is HPLN-TAN and HPLT-TAN). They are fully compatible with [https://hesperia.gsfc.nasa.gov/rhessidatacenter/complementary_data/maps/ the SSWIDL map suite] which deals with FITS files. We have prepared SSWIDL routines to convert the single time or time-series FITS files to an array of SSWIDL map structure. The scripts are available in the [https://github.com/binchensun/eovsa-tutorial/tree/master/rhessi18 Github repository of our tutorial]. For those working on Virgo, local copies are placed under /common/data/eovsa_tutorial/. Two scripts are relevant to this tutorial:
| |
| * [https://github.com/binchensun/eovsa-tutorial/blob/master/rhessi18/casa_readfits.pro casa_readfits.pro]: read an array of multi-frequency FITS files into FITS header (index) and data.
| |
| * [https://github.com/binchensun/eovsa-tutorial/blob/master/rhessi18/casa_fits2map.pro casa_fits2map.pro]: convert an array of multi-frequency FITS files into an array of SSWIDL map structure (adapted from P. T. Gallagher's hsi_fits2map.pro)
| |
|
| |
| First, start SSWIDL from command line:
| |
| <pre style="background-color:#FCEBD9;">
| |
| sswidl
| |
| </pre>
| |
|
| |
| Find and convert the fits files:
| |
| <pre style="background-color:#FCEBD9;">
| |
| ; cd to your path that contains casa_readfits.pro and casa_fits2map.pro. If on Virgo, just add the path
| |
| IDL> add_path, '/common/data/eovsa_tutorial/'
| |
| IDL> fitsdir = '/common/data/eovsa_tutorial/image_series/'
| |
| IDL> fitsfiles = file_search(fitsdir+'*_ALLBD.fits')
| |
| ; The resulting fitfiles contains all multi-frequency FITS files under "fitsdir"
| |
| ; The following command converts the fits files to an array of map structures.
| |
| ; "casa_fits2map.pro" also work well on a single FITS file.
| |
| ; (Optional) keyword "calcrms" is to calculate rms and dynamic range (SNR) of the maps in a user-defined "empty" region
| |
| ; of the map specified by "rmsxran" and "rmsyran"
| |
| ;;; Here we load 10 time frames as a demonstration
| |
| IDL> casa_fits2map, fitsfiles[45:54], eomaps [, /calcrms, rmsxran=[260., 320.], rmsyran=[-70., 0.]]
| |
| </pre>
| |
|
| |
| The resulting maps have a shape of [# of frequencies, # of polarizations, # of times]
| |
| <pre style="background-color:#FCEBD9;">
| |
| IDL> help,eomaps
| |
| EOMAPS STRUCT = -> <Anonymous> Array[30, 1, 10]
| |
| </pre>
| |
| They have compatible keywords recognized by, e.g., plot_map.pro. Example of the output map keywords:
| |
| <pre style="background-color:#FCEBD9;">
| |
| IDL> help,eomaps,/str
| |
| ** Structure <36009f8>, 24 tags, length=131336, data length=131332, refs=1:
| |
| DATA DOUBLE Array[128, 128]
| |
| XC DOUBLE 378.99175
| |
| YC DOUBLE 49.001472
| |
| DX DOUBLE 2.0000000
| |
| DY DOUBLE 2.0000000
| |
| TIME STRING '21-Aug-2017 20:20:59.500'
| |
| ID STRING 'EOVSA XX 3.413GHz'
| |
| DUR DOUBLE 3.9999997
| |
| XUNITS STRING 'arcsec'
| |
| YUNITS STRING 'arcsec'
| |
| ROLL_ANGLE DOUBLE 0.0000000
| |
| ROLL_CENTER DOUBLE Array[2]
| |
| FREQ DOUBLE 3.4125694
| |
| FREQUNIT STRING 'GHz'
| |
| STOKES STRING 'XX'
| |
| DATAUNIT STRING 'K'
| |
| DATATYPE STRING 'Brightness Temperature'
| |
| BMAJ DOUBLE 0.0097377778
| |
| BMIN DOUBLE 0.0097377778
| |
| BPA DOUBLE 0.0000000
| |
| RSUN DOUBLE 948.03989
| |
| L0 FLOAT 0.00000
| |
| B0 DOUBLE 6.9299225
| |
| COMMENT STRING 'Converted by CASA_FITS2MAP.PRO'
| |
| </pre>
| |
|
| |
| [[file:fig-specimg_idl.png|thumb|right|300px|Example multi-frequency EOVSA images at one time plotted in SSWIDL]]
| |
| An example for plotting all images at a selected time:
| |
| <pre style="background-color:#FCEBD9;">
| |
| ; choose a time pixel
| |
| IDL> plttime = anytim('2017-08-21T20:21:15')
| |
| IDL> dt = min(abs(anytim(eomaps[0,0,*].time)-plttime),tidx)
| |
| ; use maps at this time, first polarization, and all bands
| |
| IDL> eomap = reform(eomaps[*,0,tidx])
| |
| IDL> window,0,xs=1000,ys=800
| |
| IDL> !p.multi=[0,6,5]
| |
| IDL> loadct, 3
| |
| IDL> for i=0, 29 do plot_map,eomap[i],grid=5.,$
| |
| tit=string(eomap[i].freq,format='(f5.2)')+' GHz', $
| |
| chars=2.0,xran=[340.,420.],yran=[10.,90.]
| |
| </pre>
| |
|
| |
|
| =Self-calibration (optional)= | | =Self-calibration (optional)= |
| =Downloading fits files to your local system= | | =Downloading fits files to your local system= |