EOVSA Data Analysis Tutorial: Difference between revisions

From EOVSA Wiki
Jump to navigation Jump to search
No edit summary
 
(518 intermediate revisions by 4 users not shown)
Line 1: Line 1:
=Software=
=Browsing and Obtaining EOVSA data=
[[file:fig-browser.png|thumb|right|500px|Browsing EOVSA data in RHESSI Browser]]
 
Currently the most convenient way for browsing EOVSA data is through the [http://sprg.ssl.berkeley.edu/~tohban/browser/ RHESSI Browser]. First, check the "EOVSA Radio Data" box on the data selection area (top-left corner). Then select year/month/date to view the overall EOVSA dynamic spectrum. Note if a time is selected at early UTC hours (e.g., 0-3 UT), it will show the EOVSA dynamic spectrum from the previous day. Also note that EOVSA data were not commissioned for spectroscopic imaging prior to April 2017. An example of the overview EOVSA dynamic spectrum for 2017 Aug 21 is shown in the figure on the right.
 
The overview EOVSA dynamic spectra are from the median of the <b>uncalibrated</b> cross-power visibilities at a few short baselines, which are not (but a good proxy of) the total-power dynamic spectra. The effects of spatial information blended in the cross-power visibilities can be clearly seen as the "U"-shaped features throughout the day, which are due to the movement of the Sun across the sky that effectively changes the length and orientation of the baselines. Flare times can be easily seen in the EOVSA dynamic spectra, which usually appear as vertical bright features across many frequency bands. More information can be found on [http://ovsa.njit.edu/data-browsing.html this page].
 
We are currently working on a pipeline to create quicklook full-disk images for implementation at the RHESSI Browser in a way similar as the RHESSI quicklook images. Please stay tuned.
 
Once you identify the flare time, you can find the full-resolution (1-s cadence) uncalibrated visibility files (in Miriad format) at [http://www.ovsa.njit.edu/fits/IDB/ this link]. Each data file is usually 10 minutes in duration. Name convention is YYYYMMDD/IDBYYYYMMDDHHMMSS, where the time in the file name indicates the start time of the visibility data.
 
=Calibrating EOVSA data=
 
*<b>Calibration</b>: Calibrating EOVSA data is an involved process. If you are interested the [[Practical Calibration Tutorial]] describes the procedure.  We are currently working on a semi-automatic pipeline for calibrating the visibility data. At this moment, however, please contact the EOVSA team if you wish to have calibrated visibility data for a specific event. In the near future, we envision the pipeline-calibrated visibility data to be available for download on a regular basis.
 
*<b>Self-calibration</b>: For improving the quality of the imaging, self-calibration is usually needed. This is still a work-in-progress, but we have had some successful practice in self-calibrating some flare events. Please refer to [[Self-Calibrating Flare Data]] [[Tohban_EOVSA_Imaging_Tutorial_A-Z]] for examples and detailed discussions on our current practice for self-calibrating EOVSA flare data.
 
=Analyzing EOVSA data=
==Software==
We have developed two packages for EOVSA data processing and analysis:
We have developed two packages for EOVSA data processing and analysis:
* [https://github.com/binchensun/suncasa SunCASA] A wrapper around CASA for imaging and visualizing spectral imaging data of the Sun. More information about CASA can be found on [https://casa.nrao.edu/ NRAO's CASA website ]. Note, (Sun)CASA is ONLY AVAILABLE on UNIX-BASED PLATFORMS (sorry Windows users).
* [https://github.com/suncasa/suncasa SunCASA] A wrapper around [https://casa.nrao.edu/ CASA (the Common Astronomy Software Applications package)] for synthesis imaging and visualizing solar spectral imaging data. CASA is one of the leading software tool for "supporting the data post-processing needs of the next generation of radio astronomical telescopes such as ALMA and VLA", an international effort led by the [https://public.nrao.edu/ National Radio Astronomy Observatory]. The current version of CASA uses Python (2.7) interface. More information about CASA can be found on [https://casa.nrao.edu/ NRAO's CASA website ]. Note, CASA is available ONLY on UNIX-BASED PLATFORMS (and therefore, so is SunCASA).  
* [https://github.com/Gelu-Nita/GSFIT GSFIT] A IDL-widget(GUI)-based spectral fitting package called gsfit, which provides a user-friendly display of EOVSA image cubes and an interface to fast fitting codes (via platform-dependent shared-object libraries).  
* [https://github.com/Gelu-Nita/GSFIT GSFIT] A IDL-widget(GUI)-based spectral fitting package called gsfit, which provides a user-friendly display of EOVSA image cubes and an interface to fast fitting codes (via platform-dependent shared-object libraries).  


There are two approaches in accessing our software packages. One is through our Amazon AWS server (recommended for participants of the EOVSA tutorial at the [http://rhessi18.umn.edu/ RHESSI 18 Workshop]). Another is to install them on your own machine. We discuss the first approach in Section 1.1, and the second in Section 1.2.1 (SunCASA) and 1.2.2 (GSFIT).
For this tutorial, there are two approaches in accessing our software packages:
* 1. Through our '''Amazon Cloud server'''. See [[#Using Software on AWS server|Section 3.1.1]] for detailed instructions.
* 2. Install on '''your own machine'''. See [[#SunCASA Installation|Section 3.1.2]] (SunCASA) and [[#GSFIT Installation|Section 3.1.3]] (GSFIT) for instructions.
 
===Using Software on Amazon Cloud Server===
We use an [https://aws.amazon.com/lightsail/ Amazon AWS Lightsail] cloud server for lightweight data processing and testing purposes. The server has 2 CPUs, 8 GB RAM, and 160 GB SSD storage. It runs CentOS 7 (1901-01) Linux.


==Connecting to AWS server==
''NOTE: the guest account is <b>ONLY INTENDED FOR THE EOVSA TUTORIAL, WHICH WILL BE DISABLED SHORTLY AFTER THE RHESSI WORKSHOP </b> (you will be notified prior to that). '' If you wish to continue using the server for testing purposes after the tutorial, please contact [mailto:bin.chen@njit.edu Bin Chen] to set up a personal account.  
We use an Amazon AWS Lightsail server for testing purposes. The server has 2 CPUs, 8 GB RAM, and 160 GB SSD storage. It runs CentOS 7 (1901-01) Linux. Please limit your usage to '''LIGHTWEIGHT DATA PROCESSING ONLY'''.


'''NOTE: THE ACCOUNT IS ONLY INTENDED FOR THE EOVSA TUTORIAL, BUT *NOT* FOR CARRYING OUT ACTUAL DATA REDUCTION'''
''DISCLAIMER: We are NOT RESPONSIBLE for any loss, damage, or leak of your personal data on the server.''
* Obtain SSH Key from [mailto:bin.chen@njit.edu Bin Chen]
* Obtain SSH Key to the guest account from [mailto:bin.chen@njit.edu Bin Chen]
* Put it under a secure location on your own machine.
* Put it under a secure location on your own machine.
* Follow remaining directions depending on your client machine.
* Follow remaining directions depending on your client machine.


===Linux / Mac===
====Connecting via command line====
Recommend to use "~/.ssh" (create if it does not exist by "mkdir ~/.ssh").  
* Edit the permission of the SSH key and the directory you choose to place it. Here I use ~/.ssh as an example (create it first by ''mkdir ~/.ssh'' if it does not exist.)
* Edit the permission of ~/.ssh and the key (here I use ~/.ssh as the directory to place your key)
<pre>
<pre>
chmod 700 ~/.ssh
chmod 700 ~/.ssh
chmod 400 ~/.ssh/guest-virgo.pem
chmod 400 ~/.ssh/guest-virgo.pem
</pre>
</pre>
* Log on to test AWS server (password-less)
* Log on to AWS server (password-less)
<pre>
<pre>
ssh -X -i ~/.ssh/guest-virgo.pem guest@virgo.arcs.az.njit.edu
ssh -X -i ~/.ssh/guest-virgo.pem guest@virgo.arcs.az.njit.edu
</pre>
</pre>
Now you are connected to '''virgo'''.


===Windows (MobaXterm)===
====Connecting via Windows (MobaXterm)====
Recommend to use "Documents\MobaXterm\home\.ssh", which should exist if you have already installed the free MobaXterm[https://mobaxterm.mobatek.net/].
Recommend to use "Documents\MobaXterm\home\.ssh", which should exist if you have already installed the free MobaXterm[https://mobaxterm.mobatek.net/].
* Create new session, click SSH, enter virgo.arcs.az.njit.edu for Remote host, guest for username.
* Create new session, click SSH, enter virgo.arcs.az.njit.edu for Remote host, guest for username.
* On advanced SSH settings tab, click Use private key, navigate to and select file guest_key.pem
* On advanced SSH settings tab, click Use private key, navigate to and select file guest_virgo.pem
* Close setup window and click the new sessions icon, which will log you in.
* Close setup window and click the new session's icon, which will log you in.


Once the connection is setup, you will have access to SunCASA and GSFIT (included in the sswIDL installation):
====Startup SunCASA and sswIDL====
Enter SunCASA
Once you are connected to '''virgo''', you will have access to SunCASA and GSFIT (included in the sswIDL installation). To not interfere with others (who share the same "guest" account), please <b>create your own directory and work under it</b>. For easier identification, please kindly use <b>the initial of your first name and your full last name as the name of your directory</b> (here "bchen" is used as an example).
<pre>
[guest@ip-172-26-5-203 ~]$ mkdir bchen
[guest@ip-172-26-5-203 ~]$ cd bchen
</pre>
To enter SunCASA
<pre>
<pre>
[guest@ip-172-26-5-203 ~]$ suncasa
[guest@ip-172-26-5-203 ~/bchen]$ suncasa
</pre>
</pre>
Enter sswIDL
To enter sswIDL
<pre>
<pre>
[guest@ip-172-26-5-203 ~]$ sswidl
[guest@ip-172-26-5-203 ~/bchen]$ sswidl
</pre>
</pre>


==Installation on Your Own Machine==
===[[SunCASA Installation]]===
===SunCASA===
Please [http://www.ovsa.njit.edu/wiki/index.php/SunCASA_Installation follow this link] for details regarding installation of SunCASA on your own machine (only available on Unix-bases OS). This will take you to another page.
====MacOS====
We have packed a standard disk image (dmg) for installation on MacOS. It has been tested to work under Mojave (macOS v10.14). YMMV for earlier versions of macOS. Installation steps below:
* Download [https://drive.google.com/open?id=1EFOQRFyM5Ih3SJn1ugn1fKOQp1uaDIyT SunCASA disk image (SunCASA-0.7.4_Pre-release.OSX.dmg)].
* If you do not have Java SE Development Kit (JDK) installed on your Mac, please download from the [https://www.oracle.com/technetwork/java/javase/downloads/jdk12-downloads-5295953.html official site] and install it before SunCASA installation. The latest version (JDK 12) was tested to work properly. Earlier versions may also work but YMMV.
* Open the disk image file (if your browser does not do so automatically).
* Drag the SunCASA application to the Applications folder of your hard disk.
* Eject the SunCASA disk image.
* Double-click the SunCASA application to run it for the first time. If the OS does not allow you to install apps from non-Apple sources, please Change the settings in "System Preferences-> Security & Privacy -> General" and "Allow applications downloaded from: Mac App store and identified developers". If the OS still reports that the app is damaged, please turn off the OS Gatekeeper for now by running command ''sudo spctl --master-disable'' in Terminal. Then double-click the SunCASA application again.
* Initialize SunCASA. To do so, run ''!install_suncasa'' from a SunCASA prompt. This step also allow you to create symbolic links to the SunCASA version and its executables (Administrator privileges are required), which will allow you to run ''suncasa'', ''casaviewer'', ''casaplotms'', etc. from any terminal command line.
* '''Important''': Make sure to turn the OS Gatekeeper back on after the installation by running ''sudo spctl --master-enable'' in Terminal.
* Optional: To update the data repository, run !update-data from the SunCASA prompt.
* Restart SunCASA by exiting the current SunCASA prompt, and run ''suncasa'' in Terminal or double-click the SunCASA application icon.
* Now SunCASA has been up and running on your Mac.


====Linux====
===[[GSFIT Installation]]===
* To install SunCASA for Linux, we have packaged up a binary distribution of SunCASA which is available as a downloadable tar file. We have tested the package under Scientific Linux 6 (derived from RHEL 6) and CentOS 7 (equivalent to RHEL 7). They may work under some other Linux distributions (e.g. Ubuntu), but YMMV.  
Please [http://www.ovsa.njit.edu/wiki/index.php/GSFIT_Installation follow this link] for details regarding installation of GSFIT on your own machine (which supports multiple platforms). This will take you to another page.
* The following prerequisite packages that need to be installed to be able to install SunCASA. Use command "''yum install '''yourpackagename'''''" in Terminal to install.
** xauth
** libXft
** libXi
** libXrandr
** libXfixes
** libXcursor
** libXinerama
** libGL
** libXpm


* Download the packages
==Example Data and Scripts==
** [https://drive.google.com/file/d/1pQVZ4rTUsFOKps7uuLIVBf-Ir1zAgJw_/view?usp=sharing Package for RedHat6 or equivalent]
In this tutorial, we will demonstrate steps for spectral imaging (self-)calibrated visibility data for a C-class flare on 2017 Aug 21 at ~20:20 UT ([[#Spectral Imaging with SunCASA|Section 3.3]]), followed by a tutorial on visualizing and analyzing the resulting spectral imaging cube ([[#Spectral Fitting with GSFIT|Section 3.4]]).
** [https://drive.google.com/file/d/1UgKhdF9rgm5V5VSUhpAMfwRbMXm9XGk6/view?usp=sharing Package for RedHat7 or equivalent]
 
* Untar the package
The visibility data can be accessed by
<pre>tar -xzvf SunCASA-release-##version##.tar.gz</pre>
* Downloading from [https://drive.google.com/open?id=1ZMCY9Y3FTv-m0QFyBBwrvqaSuS5pAc3l this Google Drive link].
You do not have to have root or sudo permission to install or run SunCASA. The package is self-contained in the (untarred) directory "SunCASA-release-##version##". You can move/delete the SunCASA directory as you wish. But to use the executables, do the following:
* If you are on our AWS cloud server "virgo" (means for accessing the server are discussed in [[#Using Software on AWS server|Section 3.1.1]]), it is located at /common/data/eovsa_tutorial/IDB20170821201800-202300.4s.slfcaled.ms.
* All executables, including ''suncasa'' are in the ''SunCASA-release-##version##/bin'' directory. Include these executables to your path for once (examples below are in bash)
For completeness, the pre-slfcalibration visibility data can be downloaded from [https://drive.google.com/open?id=1ZMCY9Y3FTv-m0QFyBBwrvqaSuS5pAc3l this Google Drive link] (not used by this tutorial).
<pre>
cd SunCASA-release-##version##/bin
''To reduce the amount of data for this tutorial, the time resolution has been reduced by a factor of 4, to 4 s per sample.  The actual time resolution available is 1 s.'' If you are on virgo, please copy it over to your own working directory that you created earlier under the guest account (e.g., ~/bchen/).
PATH=`pwd`:$PATH
 
<pre style="background-color:#FCEBD9;">
cd ''your_working_directory''
cp -r /common/data/eovsa_tutorial/IDB20170821201800-202300.4s.slfcaled.ms ./
</pre>
</pre>
*Then you can initialize SunCASA with the command in bash
 
<pre>
Example scripts for doing the analysis covered in this tutorial are available as [https://github.com/binchensun/eovsa-tutorial/tree/master/rhessi18 this Github repository]. On Virgo, the scripts can be found under /common/data/eovsa_tutorial/.
install_suncasa
 
==Spectral Imaging with SunCASA==
EOVSA data is handled in [https://casa.nrao.edu CASA] tables system, known as a Measurement Set (MS). The actual visibility data are stored in a MAIN table that contains a number of rows, each of which is effectively a single timestamp for a single spectral window and a single baseline. Within SunCASA, in addition to everything already available in CASA, you have access to additional tools that allow you to explore and utilize the new radio dynamic spectroscopic imaging data from EOVSA. We also had success in using SunCASA for analyzing data from the [https://science.nrao.edu/facilities/vla Karl G. Jansky Very Large Array]. The steps are very similar to those described here.
To start SunCASA, use:
<pre style="background-color:#FCEBD9;">
suncasa
</pre>
</pre>
*Optional: To update the data repository, run !update-data from the SunCASA prompt.
Within SunCASA, you are under the IPython environment. Everything you know about (I)Python should be applicable here. The installation comes with frequently used packages including Matplotlib, Numpy, SciPy, AstroPy, SunPy. However, it is not very intuitive to add (compatible) Python packages within (Sun)CASA. If you choose to do so, you may want to check [http://docs.astropy.org/en/stable/install.html#installing-astropy-into-casa this page] on how to install AstroPy (which is not originally shipped with CASA) within CASA. This method is generally applicable to adding other packages (but not thoroughly tested). If you need some specific packages for your analysis, and it does not require direct interaction with (Sun)CASA, we recommend you to use the standard Python environment. On Virgo, we have installed Anaconda 3, which can be accessed by, e.g., typing "ipython" in a terminal window.
*Now SunCASA has been successfully installed on your machine. Open a new Terminal and run ''suncasa''.


===GSFIT===
===Cross-Power Dynamic Spectrum===
We have developed a IDL-widget(GUI)-based spectral fitting package called [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_GUI_Application ''gsfit''], which provides a user-friendly display of EOVSA image cubes and an interface to fast fitting codes (via platform-dependent shared-object libraries). Fits to individual spectra can be done quickly for manual investigation, while parallel/multi-core batch processing of selected blocks of data can also be performed using a command prompt application called [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFITCP_Batch_Mode_Application ''gsfitcp'']. A helper routine called [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIVIEW_GUI_Application ''gsfitview''] allows further display and investigation of the fitting results.
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


====GSFIT SSW_UPGRADE Instalation====
## define the polarizations to show (here I use XX and YY)
The release version of the IDL GSFIT package is intended to be distributed through the SSW IDL repository [https://www.lmsal.com/solarsoft/].
pol='XXYY'
Although installed as a stand-alone package, the GSFIT code relies on a series of IDL support routines that are part of the gx_simulator package. Therefore, in addition to installing the gsfit package , the gx_simulator package must be also installed. This installation can be performed by issuing an upgrade command, i.e.
## 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.


<pre>
===Quick-Look Imaging===
IDL> ssw_upgrade,/gsfit,/gx_simulator,/gen,/spawn,/loud,/passive_ftp
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>
</pre>


====GSFIT Manual SWW Installation and Setup====
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).  
However, if as of today the GSFIT package has not been yet released through SSW, you may perform a manual install by copying the directory structure located at  https://github.com/Gelu-Nita/GSFIT to your local machine SSW  directory, $ssw/packages/gsfit/'
If such manual installation is performed, the $ssw/gen/setup/setup.ssw_env script must be edited by adding or altering the following lines:
<pre>
setenv SSW_GSFIT $SSW_PACKAGES/gsfit$


setenv SSW_PACKAGES_INSTR "packages/gsfit [...]"
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>


setenv SSW_PACKAGES_ALL "$SSW_GSFIT [...]"
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>
</pre>


where [...] denote the definitions already existent in the original setup.ssw_env script.


====GSFIT SSW Instrument Setup====
===Quick-Look Imaging series===
Please note that, regardless the method chosen to install any of these two packages, to ensure proper functionality, one should make sure that sswidl.bat (on Windows platforms) or cshrc (on Unix or Mac platforms) scripts are properly updated such as to include gx_simulator and gsfit in the SSW_INSTR path declaration, i.e.
<pre style="background-color: #FCEBD9;">
<pre>
# In SunCASA
set SSW_INSTR=gx_simulator, gsfit [...]
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>
</pre>
where [..] denotes any other SSW packages already installed.


====GSFIT Manual Stand-Alone Installation====
===Batch-Mode Imaging===
Although '''not reccomended''', one may in principle choose to install the gsfit and gx_simulator packages in standalone directories that are not part of the SSW repositories. However if this (not recommended) option is used, the gsfit and gx_simulator paths should be explicitly added to the IDL !path global variable, i.e.
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>


and explicitly add this path to your IDL path structure.
[[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>
<pre>
IDL> !path='..//gsfit/idl:'+!path
################### USER INPUT GOES IN THIS BLOK ################
IDL> !path='..//gx_simulator/idl:'+!path
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>
</pre>


where '''../''' should be replaced with the explicit paths to the respective local machine repositories.
Then, run the script in SunCASA via
<pre style="background-color:#FCEBD9;">
CASA <##>: execfile('imaging_example.py')
</pre>


====GSFIT Linux / Mac OS Requirements====
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.
On Mac and Unix platforms, the gfortran compiler must be also installed to allow automatic compilation of the source code located in the ..//gsfit/unix directory, in order to generate a series of shared libraries, when gsfit is launched for the first time. WARNING: For this action to be successfully completed, the user must have writing rights to the ..//gsfit/unix directory.
<pre>                       
NAXIS  =                    4 / number of array dimensions                   
NAXIS1  =                  128                                                 
NAXIS2  =                  128                                                 
NAXIS3  =                  30                                                 
NAXIS4  =                    2 
</pre>


====GSFIT Windows OS Requirements====
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]].
The GSFIT package is distributed along with a set of dynamic link libraries located in the ..//gsfit/win subfolder, which have been compiled assuming a WIN64 architecture. For Win32, please contact us to inquire about alternative options.


=Spectral Imaging with SunCASA=
====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. 


=Spectral Fitting with GSFIT=
<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>
==GSFIT GUI Application==
The GSFIT interactive GUI application may be launched as follows


<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>
<pre>
IDL> gsfit [,nthreads]
################### 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>
</pre>


where '''nthreads''' is an optional argument that indicates the number of parallel asynchronous threads to be used when performing the fit tasks. By default, GSFT launches with only one thread, but the user may interactively add or delete threads as needed at the run-time up to the number of CPUs available on the system.
Then, run the script in SunCASA by
<pre style="background-color:#FCEBD9;">
CASA <##>: execfile('imaging_timeseries_example.py')
</pre>


[[File:GSFIT.png|border|1200px|GSFIT GUI]]
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].


===GSFIT GUI Organization and Functionality===
====Plotting Images in SSWIDL====
The GSFIT GUI is organized in three main sections:
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)


* [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Top_Main_Menu_Toolbar GSFIT Top Main Menu Toolbar], which may be used to perform most of the actions.
First, start SSWIDL from command line:
* [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Left_Panel GSFIT Left Panel], which displays a microwave map corresponding to the observing frequency and time instance selected via two horizontal scroll bars.
<pre  style="background-color:#FCEBD9;">
* [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Middle_Panel GSFIT Middle Panel], which displays the microwave spectrum corresponding to the selected pixel in the right panel map.
sswidl
* [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Right_Panel GSFIT Right Panel], which is contain three main elements:
</pre>
**[http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_List_Scrollbar GSFIT List Scrollbar], which allows the user to locate spatially and temporarily the fit already performed and inspect the fit results
**[http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Settings GSFIT Settings page], which allows the user to setup a series of the spectral fit options and perform some actions via a local Toolbar menu.
**[http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Task_Queue GSFIT Task Queue page], which displays a list of completed, active, and pending fitting tasks.
====GSFIT Top Main Menu Toolbar====
This toolbar contain a series of action buttons described below.


=====[[File:open.png|border]] Upload EOVSA Map Cube=====
Find and convert the fits files:
EOVSA map cubes are IDL sav files containg a [ntimes x nfrequencies] brightness temperature (Tb) map structure.  
<pre style="background-color:#FCEBD9;">
Follow this [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#Spectral_Imaging_with_SunCASA link] to learn how such EOVSA map cubes are produced. When uploaded, GSFIT converts these Tb maps to Solar Flux (total intensity or circular polarization) maps.
; cd to your path that contains casa_readfits.pro and casa_fits2map.pro. If on Virgo, just add the path
=====[[File:color.png|border]] Select Color Table=====
IDL> add_path, '/common/data/eovsa_tutorial/'
Use this button to select a color table for displaying the microwave maps and spectra. GSFIT uses as default the IDL color table #39 (Rainbow+White)
IDL> fitsdir = '/common/data/eovsa_tutorial/image_series/'
=====[[File:mousemode.png|border]]Select Mouse Mode=====
IDL> fitsfiles = file_search(fitsdir+'*_ALLBD.fits')
Use this set of mutual exclusive radio buttons to switch between four modes of mouse tool operation.
; The resulting fitfiles contains all multi-frequency FITS files under "fitsdir"
*'''FOV/Pixel Selection Mode''': Use this mouse mode to zoom in and out the field of view (FOV) of the microwave maps displayed on the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Left_Panel GSFIT Left Panel], or to select the pixel for which the microwave spectrum is displayed in the[http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Middle_Panel GSFIT Middle Panel].
; The following command converts the fits files to an array of map structures.  
**To zoom in, press, hold, move, and release the left mouse button to select the desired rectangular FOV (rubber-band selection process). To revert to the full FOV, use a single left-button click.
; "casa_fits2map.pro" also work well on a single FITS file.  
**Use the right button to select a map pixel on the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Left_Panel GSFIT Left Panel]. The currently selected pixel is indicated by a set og horizontal and vertical dotted lines. Alternatively, the pixel may be selected using the ''Cursor X'' and ''Cursor Y'' numerical controls located below the map plot area.
; (Optional) keyword "calcrms" is to calculate rms and dynamic range (SNR) of the maps in a user-defined "empty" region
*'''Rectangular ROI Selection Mode''': Use this mouse mode to define a rectangular region of interest (ROI) through a left button rubber-band selection process. The selected ROI is preserved and displayed as long it is not replaced by another ROI selection, or deleted by a single left button click.
;        of the map specified by "rmsxran" and "rmsyran"
*'''Polygonal ROI Selection Mode''': Use this mouse mode to define the vertices of polygonal ROI through a series of sequential left button clicks. To finalize the process, use a right button click to define the last vertex, which will be automatically connected with the first one. The selected ROI is preserved and displayed as long it is not replaced by another ROI selection, or deleted by a single left button click.
;;; Here we load 10 time frames as a demonstration
*'''Free-hand ROI Selection Mode''': Use this mouse mode to define an arbitrarily shaped ROI contour through a continues movement of the mouse cursor while the left button is pressed. To finalize the process, release the left button to define the last vertex of the countour, which will be automatically connected with the first one. The selected ROI is preserved and displayed as long it is not replaced by another ROI selection, or deleted by a single left button click.
IDL> casa_fits2map, fitsfiles[45:54], eomaps [, /calcrms, rmsxran=[260., 320.], rmsyran=[-70., 0.]]
</pre>


=====[[File:SaveROI.png|border]]Save ROI=====
The resulting maps have a shape of [# of frequencies, # of polarizations, # of times]
Use this button to save to an IDL file the current ROI selection as a set of two arrays, xroi and yroi, containing its vertex heliocentric coordinates.
<pre  style="background-color:#FCEBD9;">
=====[[File:RestoreROI.png|border]]Restore ROI=====
IDL> help,eomaps
Use this button to import from an IDL file a set of xroi and yroi ROI vertex coordinates previously saved by GSFIT, or generated by alternative means.
EOMAPS          STRUCT    = -> <Anonymous> Array[30, 1, 10]
=====[[File:FreqRange.png|border]]Select Fit Frequency Range=====
</pre>
Use this set off drop lists controls to select the frequency range to be used for spectral fitting. By default, the full frequency range is selected when a map data cube is uploaded, but the user may choose to restrict this default frequency range for various reasons, such as data quality, or apparent source nonuniformity.
They have compatible keywords recognized by, e.g., plot_map.pro. Example of the output map keywords:
=====[[File:FitPixel.png|border]]Fit Selected Pixel=====
<pre  style="background-color:#FCEBD9;">
Use this button to fit the microwave spectrum displayed on the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Middle_Panel GSFIT Middle Panel] using the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#Select_Fit_Frequency_Range selected frequency range] and the settings displayed on the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Settings GSFIT Settings] panel.
IDL> help,eomaps,/str
=====[[File:FitRange.png|border]]Fit Range=====
** Structure <36009f8>, 24 tags, length=131336, data length=131332, refs=1:
Use this button to perform one or more of the following actions:
  DATA            DOUBLE    Array[128, 128]
*Define a rectangular ROI starting at the current cursor position, if no ROI has been previously defined
  XC              DOUBLE          378.99175
*Select a time range and time resolution for producing a sequence of fits for each ROI pixel
  YC              DOUBLE          49.001472
*Ignore the current ROI and Time Range selections and choose instead to fit the entire map cube
  DX              DOUBLE          2.0000000
*Use or change the current FIT Frequency Range selection
  DY              DOUBLE          2.0000000
*Enqueue and start processing the fitting of all ROI pixels over the user-defined time interval (default)
  TIME            STRING    '21-Aug-2017 20:20:59.500'
*Enqueue and Wait for future processing of the fit task of all ROI pixels over the user-defined time interval
  ID              STRING    'EOVSA XX 3.413GHz'
*Generate and save the fit task for future batch processing by the GSFIT Command Prompt ([http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFITCP_Batch_Mode_Application GSFITCP]) application.
  DUR            DOUBLE          3.9999997
Depending on whether or not a ROI is already defined, the user is exposed to a different dialog box to define a fitting task.
  XUNITS          STRING    'arcsec'
<gallery widths=500px heights=400px>
  YUNITS          STRING    'arcsec'
File:FitCube.png|Fit Task Dialog Box if no ROI has been already defined
  ROLL_ANGLE      DOUBLE          0.0000000
File:FitROI.png|Fit Task Dialog Box if a ROI has been already defined
  ROLL_CENTER    DOUBLE    Array[2]
</gallery>
  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>


To enqueue the fit tasks for immediate or future processing by GSFIT, the following rules are followed:
[[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>


*All ROI spectra from a given time frame are processed simultaneously by being as evenly as possible assigned to all active asynchronous parallel threads
==Spectral Fitting with GSFIT==
In this section we discuss using GSFIT (now a package of SSWIDL) to perform spectral fit based on the resulting FITS images produced from the previous steps. We will use the FITS images produced in [[#Producing a time series of 30-band maps (a 4-d cube)|Section 3.3.3.2]] as the input. On Virgo, we prepared an IDL save file located at /common/data/eovsa_tutorial/21August2017.sav. In addition, this folder contains three more demo files that may be used to demonstrate the GSFIT functionality.
For those intending to use their own machines for this part of the tutorial, these IDL fits files may be downloaded from [https://www.dropbox.com/sh/tjez92wfph1jp7x/AAAJBO6AJ0GC5H3GCZypc2BYa?dl=0 this location]


*The individual time frames are enqueued and processed sequentially
===GSFIT GUI Application===
The GSFIT GUI application may be launched as follows


An example about how the tasks are scheduled is provided in the  [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Task_Queue GSFIT Task Queue] section.
<pre>
IDL> gsfit [,nthreads]
</pre>


'''NOTE:''' If a fit task is saved for future batch processing, the scheduling is performed  by the GSFIT Command Prompt ([http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFITCP_Batch_Mode_Application GSFITCP]) application using the number of active parallel threads available to it.
where '''nthreads''' is an optional argument that indicates the number of parallel asynchronous threads to be used when performing the fit tasks. By default, GSFIT launches with only one thread, but the user may interactively add or delete threads as needed at the run-time up to the number of CPUs available on the system (<b>Note, if you are working on the AWS server, please use only one thread</b>). After some delay while the interface loads, the GUI below should appear.


=====[[File:ProcessPendingTasks.png|border]]Process Pending Tasks=====
<gallery  widths=600px heights=300px>
Use this button to start processing all pending tasks in the [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Task_Queue GSFIT Queue]
File:GSFITwin.png|GSFIT GUI appearance on Windows OS
File:GSFITunix.png|GSFIT GUI appearance on Linux/Mac OS
</gallery>


=====[[File:AbortDeleteTasks.png|border]]Abort/Delete Tasks=====
A detailed description of the GSFIT functionality is provided in the page linked below.
Use this button to select one of the options from the Delete Options dialog box.
====[[GSFIT Help]]====


[[File:DeleteOptions.png|400px|Abort/Delete Tasks Dialog]]
===GSFITCP Batch Mode Application===
GSFITCP is the command prompt counterpart of the GSFIT GUI microwave spectral fitting application. Once started, GSFITCP is designed to run in unattended mode until all fitting tasks assigned to it are completed and log on the disk in an user-defined *.log file, the content of which may be visualized using the GSFITVIEW GUI application. When run remotely on an Linux/Mac platform, the GSFITCP may be launched on a detached [https://linuxize.com/post/how-to-use-linux-screen/ screen], which allows the remote user to logout without stopping the process in which GSFITCP runs.


====GSFIT Left Panel====
GSFITCP is launched using the following call:
<gallery widths=500px heights=500px>
<pre>
File:LeftPanel.png|A: Selected EOVSA map displayed with no Integrated Flux Threshold applied
IDL> gsfitcp, taskfile, nthreads, /start
File:LeftPanelTh.png|B: Selected EOVSA map displayed with 1 sfu Integrated Flux Threshold applied
</pre>
</gallery>
where '''''taskfile''''' is a path to a file in which a GSFIT task has been previously saved, as explained in the [[GSFIT Help]] page, '''''nthreads''''' is an optional argument indicating the number of parallel asynchronous threads to be used, and the optional keyword '''''/start''''', if set, requests immediate start of the batch processing.
 
A detailed description of the GSFITCP functionality is provided in the page linked below.
====[[GSFITCP Help]]====


====GSFIT Middle Panel====
===GSFITVIEW GUI Application===
<gallery mode="packed-hover">
The GSFITVIEW GUI application may be launched as follows
Image:MiddlePanel.png|A:
<pre>
Image:MiddlePanelG.png|B:
IDL> gsfitview [,gsfitmaps]
Image:MiddlePanelFC.png|C:
</pre>
</gallery>
where the optional '''gsfitmaps''' argument is either the filename of an IDL ''*.sav'' file containg a GSFIT Parameter Map Cube structure produced by the GSFIT or GSFITCP applications, or an already restored such structure.
<gallery widths=600px heights=400px>
File:GSFITVIEWwin.png
File:GSFITVIEWunix.png
]</gallery>


====GSFIT Right Panel====
===== GSFIT List Scrollbar=====
=====GSFIT Settings=====
=====GSFIT Task Queue=====
This page displays the current status of the already processed, active, or pending fit tasks.
The [http://www.ovsa.njit.edu/wiki/index.php/EOVSA_Data_Analysis_Tutorial#GSFIT_Top_Main_Menu_Toolbar GSFIT top menu toolbar] provides several options for controlling the GSFIT Task Queue.


'''Here we provide an example of how the tasks are enqueued and processed in two diffrent circumstances:'''
A detailed description of the GSFITVIEW functionality is provided in the page linked below.  
 
====[[GSFITVIEW Help]]====
If a 4x4 ROI; 4 time frames is enqueued while only 1 parallel thread is active, only 16 pixels from the same time frame are processed at any given time.
If, instead, 4 active threads are used, batches of 4 pixels are processed by 4 asynchronous threads at any given time. Note that, in the case of 4 active threads, the time needed to fully process a given time frame is dictated by the longest time needed to process one of the 4 pixel batched of that frame, which is about half the time needed to process all 16 pixels in a single thread.  
<gallery widths=450px heights=450px>
File:Queue1thread.png|A: GSFIT Queue active processing: 4x4 ROI; 4 time frames; 1 parallel thread
File:Queue4threads.png|B: GSFIT Fit Queue active processing: 4x4 ROI; 4 time frames;  4 parallel threads
File:Queue4threadsAll.png|C: GSFIT Fit Queue finalized task: 4x4 ROI; 4 time frames;  4 parallel threads
</gallery>


==GSFIVIEW GUI Application==
===[[GSFIT Data Format]]===
==GSFITCP Batch Mode Application==

Latest revision as of 18:02, 20 May 2022

Browsing and Obtaining EOVSA data

Browsing EOVSA data in RHESSI Browser

Currently the most convenient way for browsing EOVSA data is through the RHESSI Browser. First, check the "EOVSA Radio Data" box on the data selection area (top-left corner). Then select year/month/date to view the overall EOVSA dynamic spectrum. Note if a time is selected at early UTC hours (e.g., 0-3 UT), it will show the EOVSA dynamic spectrum from the previous day. Also note that EOVSA data were not commissioned for spectroscopic imaging prior to April 2017. An example of the overview EOVSA dynamic spectrum for 2017 Aug 21 is shown in the figure on the right.

The overview EOVSA dynamic spectra are from the median of the uncalibrated cross-power visibilities at a few short baselines, which are not (but a good proxy of) the total-power dynamic spectra. The effects of spatial information blended in the cross-power visibilities can be clearly seen as the "U"-shaped features throughout the day, which are due to the movement of the Sun across the sky that effectively changes the length and orientation of the baselines. Flare times can be easily seen in the EOVSA dynamic spectra, which usually appear as vertical bright features across many frequency bands. More information can be found on this page.

We are currently working on a pipeline to create quicklook full-disk images for implementation at the RHESSI Browser in a way similar as the RHESSI quicklook images. Please stay tuned.

Once you identify the flare time, you can find the full-resolution (1-s cadence) uncalibrated visibility files (in Miriad format) at this link. Each data file is usually 10 minutes in duration. Name convention is YYYYMMDD/IDBYYYYMMDDHHMMSS, where the time in the file name indicates the start time of the visibility data.

Calibrating EOVSA data

  • Calibration: Calibrating EOVSA data is an involved process. If you are interested the Practical Calibration Tutorial describes the procedure. We are currently working on a semi-automatic pipeline for calibrating the visibility data. At this moment, however, please contact the EOVSA team if you wish to have calibrated visibility data for a specific event. In the near future, we envision the pipeline-calibrated visibility data to be available for download on a regular basis.
  • Self-calibration: For improving the quality of the imaging, self-calibration is usually needed. This is still a work-in-progress, but we have had some successful practice in self-calibrating some flare events. Please refer to Self-Calibrating Flare Data Tohban_EOVSA_Imaging_Tutorial_A-Z for examples and detailed discussions on our current practice for self-calibrating EOVSA flare data.

Analyzing EOVSA data

Software

We have developed two packages for EOVSA data processing and analysis:

  • SunCASA A wrapper around CASA (the Common Astronomy Software Applications package) for synthesis imaging and visualizing solar spectral imaging data. CASA is one of the leading software tool for "supporting the data post-processing needs of the next generation of radio astronomical telescopes such as ALMA and VLA", an international effort led by the National Radio Astronomy Observatory. The current version of CASA uses Python (2.7) interface. More information about CASA can be found on NRAO's CASA website . Note, CASA is available ONLY on UNIX-BASED PLATFORMS (and therefore, so is SunCASA).
  • GSFIT A IDL-widget(GUI)-based spectral fitting package called gsfit, which provides a user-friendly display of EOVSA image cubes and an interface to fast fitting codes (via platform-dependent shared-object libraries).

For this tutorial, there are two approaches in accessing our software packages:

Using Software on Amazon Cloud Server

We use an Amazon AWS Lightsail cloud server for lightweight data processing and testing purposes. The server has 2 CPUs, 8 GB RAM, and 160 GB SSD storage. It runs CentOS 7 (1901-01) Linux.

NOTE: the guest account is ONLY INTENDED FOR THE EOVSA TUTORIAL, WHICH WILL BE DISABLED SHORTLY AFTER THE RHESSI WORKSHOP (you will be notified prior to that). If you wish to continue using the server for testing purposes after the tutorial, please contact Bin Chen to set up a personal account.

DISCLAIMER: We are NOT RESPONSIBLE for any loss, damage, or leak of your personal data on the server.

  • Obtain SSH Key to the guest account from Bin Chen
  • Put it under a secure location on your own machine.
  • Follow remaining directions depending on your client machine.

Connecting via command line

  • Edit the permission of the SSH key and the directory you choose to place it. Here I use ~/.ssh as an example (create it first by mkdir ~/.ssh if it does not exist.)
chmod 700 ~/.ssh
chmod 400 ~/.ssh/guest-virgo.pem
  • Log on to AWS server (password-less)
ssh -X -i ~/.ssh/guest-virgo.pem guest@virgo.arcs.az.njit.edu

Now you are connected to virgo.

Connecting via Windows (MobaXterm)

Recommend to use "Documents\MobaXterm\home\.ssh", which should exist if you have already installed the free MobaXterm[1].

  • Create new session, click SSH, enter virgo.arcs.az.njit.edu for Remote host, guest for username.
  • On advanced SSH settings tab, click Use private key, navigate to and select file guest_virgo.pem
  • Close setup window and click the new session's icon, which will log you in.

Startup SunCASA and sswIDL

Once you are connected to virgo, you will have access to SunCASA and GSFIT (included in the sswIDL installation). To not interfere with others (who share the same "guest" account), please create your own directory and work under it. For easier identification, please kindly use the initial of your first name and your full last name as the name of your directory (here "bchen" is used as an example).

[guest@ip-172-26-5-203 ~]$ mkdir bchen
[guest@ip-172-26-5-203 ~]$ cd bchen

To enter SunCASA

[guest@ip-172-26-5-203 ~/bchen]$ suncasa

To enter sswIDL

[guest@ip-172-26-5-203 ~/bchen]$ sswidl

SunCASA Installation

Please follow this link for details regarding installation of SunCASA on your own machine (only available on Unix-bases OS). This will take you to another page.

GSFIT Installation

Please follow this link for details regarding installation of GSFIT on your own machine (which supports multiple platforms). This will take you to another page.

Example Data and Scripts

In this tutorial, we will demonstrate steps for spectral imaging (self-)calibrated visibility data for a C-class flare on 2017 Aug 21 at ~20:20 UT (Section 3.3), followed by a tutorial on visualizing and analyzing the resulting spectral imaging cube (Section 3.4).

The visibility data can be accessed by

  • Downloading from this Google Drive link.
  • If you are on our AWS cloud server "virgo" (means for accessing the server are discussed in Section 3.1.1), it is located at /common/data/eovsa_tutorial/IDB20170821201800-202300.4s.slfcaled.ms.

For completeness, the pre-slfcalibration visibility data can be downloaded from this Google Drive link (not used by this tutorial).

To reduce the amount of data for this tutorial, the time resolution has been reduced by a factor of 4, to 4 s per sample. The actual time resolution available is 1 s. If you are on virgo, please copy it over to your own working directory that you created earlier under the guest account (e.g., ~/bchen/).

cd ''your_working_directory''
cp -r /common/data/eovsa_tutorial/IDB20170821201800-202300.4s.slfcaled.ms ./

Example scripts for doing the analysis covered in this tutorial are available as this Github repository. On Virgo, the scripts can be found under /common/data/eovsa_tutorial/.

Spectral Imaging with SunCASA

EOVSA data is handled in CASA tables system, known as a Measurement Set (MS). The actual visibility data are stored in a MAIN table that contains a number of rows, each of which is effectively a single timestamp for a single spectral window and a single baseline. Within SunCASA, in addition to everything already available in CASA, you have access to additional tools that allow you to explore and utilize the new radio dynamic spectroscopic imaging data from EOVSA. We also had success in using SunCASA for analyzing data from the Karl G. Jansky Very Large Array. The steps are very similar to those described here. To start SunCASA, use:

suncasa

Within SunCASA, you are under the IPython environment. Everything you know about (I)Python should be applicable here. The installation comes with frequently used packages including Matplotlib, Numpy, SciPy, AstroPy, SunPy. However, it is not very intuitive to add (compatible) Python packages within (Sun)CASA. If you choose to do so, you may want to check this page on how to install AstroPy (which is not originally shipped with CASA) within CASA. This method is generally applicable to adding other packages (but not thoroughly tested). If you need some specific packages for your analysis, and it does not require direct interaction with (Sun)CASA, we recommend you to use the standard Python environment. On Virgo, we have installed Anaconda 3, which can be accessed by, e.g., typing "ipython" in a terminal window.

Cross-Power Dynamic Spectrum

The first module we introduce is 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 time range, spectral windows/channels, antenna baseline, or uvrange. The selection syntax follows the CASA convention. More information of CASA selection syntax may be found in the above links or the Measurement Selection Syntax.

Figure 1: EOVSA cross power dynamic spectrum at stokes XX and YY
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

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 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).

EOVSA full-Sun single-band quicklook image
## 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)

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 qlookplot, it is easy to engage solar data from SDO/AIA in the summary plot. Setting plotaia=True in the qlookplot command will download SDO/AIA data at the given time to current directory and add it to the summary plot.

EOVSA full-Sun single-band quicklook image with SDO/AIA as background
## in SunCASA
ql.qlookplot(vis=msfile, specfile=specfile, timerange=timerange, spw=spw, \
    stokes=stokes, plotaia=True)

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.

                        
NAXIS   =                    4                                                  
NAXIS1  =                  512/ Nx
NAXIS2  =                  512/ Ny                                             
NAXIS3  =                    1/  number of frequency                                           
NAXIS4  =                    2/  number of polarization

By default, 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.

EOVSA multi-frequency synthesis quicklook image
## 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)

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).

EOVSA multi-band quicklook image
## 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)

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.

                         
NAXIS   =                    4 / number of array dimensions                     
NAXIS1  =                  128                                                  
NAXIS2  =                  128                                                  
NAXIS3  =                   30                                                  
NAXIS4  =                    2  


Quick-Look Imaging series

# 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)

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 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.

# in SunCASA
CASA <##>: cd your_working_directory
CASA <##>: !cp /common/data/eovsa_tutorial/imaging_example.py ./
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.

################### 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

Then, run the script in SunCASA via

CASA <##>: execfile('imaging_example.py')

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.

                         
NAXIS   =                    4 / number of array dimensions                     
NAXIS1  =                  128                                                  
NAXIS2  =                  128                                                  
NAXIS3  =                   30                                                  
NAXIS4  =                    2  

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 Section 3.3.3.3.

Producing a time series of 30-band maps (a 4-d cube)

An example script can be found at 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.

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.

# in SunCASA
CASA <##>: cd your_working_directory
CASA <##>: !cp /common/data/eovsa_tutorial/imaging_timeseries_example.py ./

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.

Example multi-frequency images at one time frame in the output javascript movie
################### 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

Then, run the script in SunCASA by

CASA <##>: execfile('imaging_timeseries_example.py')

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 the previous section. The summary yields the time, frequency, and path to every image.

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

The last block of the example script uses SunPy.map to generate a series of plots and wraps them as a 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 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 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:

  • casa_readfits.pro: read an array of multi-frequency FITS files into FITS header (index) and data.
  • 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:

sswidl

Find and convert the fits files:

; 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.]]

The resulting maps have a shape of [# of frequencies, # of polarizations, # of times]

IDL> help,eomaps
EOMAPS          STRUCT    = -> <Anonymous> Array[30, 1, 10]

They have compatible keywords recognized by, e.g., plot_map.pro. Example of the output map keywords:

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'
Example multi-frequency EOVSA images at one time plotted in SSWIDL

An example for plotting all images at a selected time:

; 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.]

Spectral Fitting with GSFIT

In this section we discuss using GSFIT (now a package of SSWIDL) to perform spectral fit based on the resulting FITS images produced from the previous steps. We will use the FITS images produced in Section 3.3.3.2 as the input. On Virgo, we prepared an IDL save file located at /common/data/eovsa_tutorial/21August2017.sav. In addition, this folder contains three more demo files that may be used to demonstrate the GSFIT functionality. For those intending to use their own machines for this part of the tutorial, these IDL fits files may be downloaded from this location

GSFIT GUI Application

The GSFIT GUI application may be launched as follows

IDL> gsfit [,nthreads]

where nthreads is an optional argument that indicates the number of parallel asynchronous threads to be used when performing the fit tasks. By default, GSFIT launches with only one thread, but the user may interactively add or delete threads as needed at the run-time up to the number of CPUs available on the system (Note, if you are working on the AWS server, please use only one thread). After some delay while the interface loads, the GUI below should appear.

A detailed description of the GSFIT functionality is provided in the page linked below.

GSFIT Help

GSFITCP Batch Mode Application

GSFITCP is the command prompt counterpart of the GSFIT GUI microwave spectral fitting application. Once started, GSFITCP is designed to run in unattended mode until all fitting tasks assigned to it are completed and log on the disk in an user-defined *.log file, the content of which may be visualized using the GSFITVIEW GUI application. When run remotely on an Linux/Mac platform, the GSFITCP may be launched on a detached screen, which allows the remote user to logout without stopping the process in which GSFITCP runs.

GSFITCP is launched using the following call:

IDL> gsfitcp, taskfile, nthreads, /start

where taskfile is a path to a file in which a GSFIT task has been previously saved, as explained in the GSFIT Help page, nthreads is an optional argument indicating the number of parallel asynchronous threads to be used, and the optional keyword /start, if set, requests immediate start of the batch processing.

A detailed description of the GSFITCP functionality is provided in the page linked below.

GSFITCP Help

GSFITVIEW GUI Application

The GSFITVIEW GUI application may be launched as follows

IDL> gsfitview [,gsfitmaps]

where the optional gsfitmaps argument is either the filename of an IDL *.sav file containg a GSFIT Parameter Map Cube structure produced by the GSFIT or GSFITCP applications, or an already restored such structure.


A detailed description of the GSFITVIEW functionality is provided in the page linked below.

GSFITVIEW Help

GSFIT Data Format