Merge pull request #310 from sageyu123/master

Major Update: Integrate WSClean into EOVSA Synoptic Imaging Pipeline
and fix an issue with ephemeris retrieval in task_ptclean6.py

Author: sageyu123

suncasa/dspec/dspec.py

@@ -1315,7 +1315,7 @@ def format_coord(x, y):
                 elif pol == 'IV':
                     spec_plt_1 = I_plot
                     spec_plt_2 = V_plot
-                    cmap2 = 'gray'
+                    cmap2 = 'RdBu_r'
                     if (vmax2 is None) and (vmin2 is None):
                         vmax2 = np.nanmax(np.abs(spec_plt_2))
                         vmin2 = -vmax2
@@ -1328,7 +1328,7 @@ def format_coord(x, y):
                     # this is for Stokes I + polarization degree
                     spec_plt_1 = I_plot
                     spec_plt_2 = V_plot / I_plot
-                    cmap2 = 'gray'
+                    cmap2 = 'RdBu_r'
                     if (vmax2 is None) and (vmin2 is None):
                         vmax2 = 1.
                         vmin2 = -1.

suncasa/eovsa/eovsa_diskmodel.py

@@ -1102,11 +1102,24 @@ def pipeline_run(vis, outputvis='', workdir=None, slfcaltbdir=None, imgoutdir=No
         return None
 
     # Copy original ms to local directory
-    if os.path.exists(os.path.basename(vis)):
-        shutil.rmtree(os.path.basename(vis))
-    print('Copy {} to working directory {}/'.format(vis, os.getcwd()))
-    shutil.copytree(vis, os.path.basename(vis))
-    vis = os.path.basename(vis)
+    if vis.lower().endswith('.ms'):
+        if os.path.exists(os.path.basename(vis)):
+            shutil.rmtree(os.path.basename(vis))
+        print('Copy {} to working directory {}/'.format(vis, os.getcwd()))
+        shutil.copytree(vis, os.path.basename(vis))
+        vis = os.path.basename(vis)
+    elif vis.lower().endswith('.tar.gz'):
+        if os.path.exists(os.path.basename(vis)):
+            shutil.rmtree(os.path.basename(vis))
+        if os.path.exists(os.path.basename(vis.rstrip('.tar.gz'))):
+            shutil.rmtree(os.path.basename(vis.rstrip('.tar.gz')))
+        print(f'Extracting {vis} to working directory {os.getcwd()}/')
+        os.system(f'tar -xzf {vis} -C {os.getcwd()}')
+        vis = os.path.basename(vis.rstrip('.tar.gz'))
+    else:
+        print('Input vis file must be either a .ms or .tar.gz file.')
+        return None
+
     # Generate calibrated visibility by self calibrating on the solar disk
     ##ms_slfcaled, diskxmlfile = disk_slfcal(vis, slfcaltbdir=slfcaltbdir)
     flagmanager(vis, mode='save', versionname='pipeline_init')

suncasa/eovsa/eovsa_fitsutils.py

@@ -94,138 +94,90 @@ def rewriteImageFits(datestr, verbose=False, writejp2=False, overwritejp2=False,
                 ndfits.write_j2000_image(fj2name, data[::-1, :], hdu.header)
     return
 
+def main(dateobj=None, ndays=1, overwritejp2=False, overwritefits=False):
+    """
+    Main pipeline for creating compressed FITS and JP2 files of EOVSA daily full-disk images.
+
+    :param dateobj: The starting datetime for processing. If None, defaults to two days before now.
+    :type dateobj: datetime, optional
+    :param ndays: Number of days to process (spanning from dateobj - ndays + 1 to dateobj), defaults to 1.
+    :type ndays: int, optional
+    :param overwritejp2: If True, overwrite existing EOVSA JP2 files, defaults to False.
+    :type overwritejp2: bool, optional
+    :param overwritefits: If True, overwrite existing EOVSA FITS files, defaults to False.
+    :type overwritefits: bool, optional
+    :raises Exception: If an error occurs during processing.
+    :return: None
+    :rtype: None
+    """
+    from datetime import timedelta
+    import numpy as np
+    from astropy.time import Time
 
-def main(year=None, month=None, day=None, ndays=1, overwritejp2=False, overwritefits=False):
-    # tst = datetime.strptime("2017-04-01", "%Y-%m-%d")
-    # ted = datetime.strptime("2019-12-31", "%Y-%m-%d")
-    if year:
-        ted = datetime(year, month, day)
-    else:
+    # Use dateobj if provided; otherwise, default to two days before now.
+    if dateobj is None:
         ted = datetime.now() - timedelta(days=2)
+    else:
+        ted = dateobj
+
+    # Compute the start date (tst) for processing based on ndays.
     tst = Time(np.fix(Time(ted).mjd) - ndays + 1, format='mjd').datetime
-    print("Running pipeline_fitsutils for date from {} to {}".format(tst.strftime("%Y-%m-%d"),
-                                                                     ted.strftime("%Y-%m-%d")))
+    print("Running pipeline_fitsutils for date from {} to {}.".format(
+        tst.strftime("%Y-%m-%d"), ted.strftime("%Y-%m-%d")))
     dateobs = tst
     while dateobs <= ted:
         datestr = dateobs.strftime("%Y-%m-%d")
-        rewriteImageFits(datestr, verbose=True, writejp2=True, overwritejp2=overwritejp2, overwritefits=overwritefits)
+        rewriteImageFits(datestr, verbose=True, writejp2=True,
+                          overwritejp2=overwritejp2, overwritefits=overwritefits)
         dateobs = dateobs + timedelta(days=1)
 
 
 if __name__ == '__main__':
-    '''
-    Name: 
-    eovsa_fitsutils --- pipeline for created the compressed fits and jp2 files of EOVSA daily full-disk images.
-
-    Synopsis:
-    eovsa_fitsutils.py [options]... [DATE_IN_YY_MM_DD]
-
-    Description:
-    Plot EOVSA daily full-disk images at multi frequencies of the date specified
-    by DATE_IN_YY_MM_DD (or from ndays before the DATE_IN_YY_MM_DD if option --ndays/-n is provided).
-    If DATE_IN_YY_MM_DD is omitted, it will be set to 2 days before now by default. 
-    The are no mandatory arguments in this command.
-
-    -c, --clearcache
-            Remove temporary files
-
-    -n, --ndays
-            Processing the date spanning from DATE_IN_YY_MM_DD-ndays to DATE_IN_YY_MM_DD. Default is 30
-
-    -o, --overwritejp2
-            If True, overwrite eovsa jp2 files.
-            Syntax: True, False, T, F, 1, 0
-
-    -O, --overwritefits
-            If True, overwrite eovsa fits files.
-            Syntax: True, False, T, F, 1, 0                          
-
-
-    Example: 
-    eovsa_fitsutils.py -c True -n 2 -o True -O True 2020 06 10
-    '''
-    import sys
-    import numpy as np
-    import getopt
+    import argparse
     from datetime import datetime, timedelta
+    from astropy.time import Time
+
+    parser = argparse.ArgumentParser(
+        description='Pipeline for creating compressed FITS and JP2 files of EOVSA daily full-disk images.'
+    )
+    # Default date is set to two days before the current date at 20:00 UT (YYYY-MM-DDT20:00).
+    default_date = (datetime.now() - timedelta(days=2)).strftime('%Y-%m-%dT20:00')
+    parser.add_argument(
+        '--date', type=str, default=default_date,
+        help='Date to process in YYYY-MM-DDT20:00 format, defaults to 20:00 UT two days before the current date.'
+    )
+    parser.add_argument(
+        '--ndays', type=int, default=1,
+        help='Process data spanning from DATE minus ndays to DATE (default: 1 day).'
+    )
+    parser.add_argument(
+        '--overwritejp2', action='store_true',
+        help='Overwrite existing EOVSA JP2 files.'
+    )
+    parser.add_argument(
+        '--overwritefits', action='store_true',
+        help='Overwrite existing EOVSA FITS files.'
+    )
+    # Optional positional date arguments: year month day (overrides --date if provided)
+    parser.add_argument(
+        'date_args', type=int, nargs='*',
+        help='Optional date arguments: year month day. If provided, overrides --date.'
+    )
+
+    args = parser.parse_args()
+
+    # Determine the processing date.
+    if len(args.date_args) == 3:
+        year, month, day = args.date_args
+        dateobj = datetime(year, month, day, 20)  # Use 20:00 UT for the specified date.
+    else:
+        dateobj = Time(args.date).datetime
 
-    # import subprocess
-    # shell = subprocess.check_output('echo $0', shell=True).decode().replace('\n', '').split('/')[-1]
-    # print("shell " + shell + " is using")
-
-    print(sys.argv)
-    year = None
-    month = None
-    day = None
-    ndays = 1
-    clearcache = True
-    opts = []
-    overwritejp2 = False
-    overwritefits = False
-    try:
-        argv = sys.argv[1:]
-        opts, args = getopt.getopt(argv, "c:n:o:O:", ['clearcache=', 'ndays=', 'overwritejp2=', 'overwritefits='])
-        print(opts, args)
-        for opt, arg in opts:
-            if opt in ['-c', '--clearcache']:
-                if arg in ['True', 'T', '1']:
-                    clearcache = True
-                elif arg in ['False', 'F', '0']:
-                    clearcache = False
-                else:
-                    clearcache = np.bool_(arg)
-            elif opt in ('-n', '--ndays'):
-                ndays = int(arg)
-            elif opt in ('-o', '--overwritejp2'):
-                if arg in ['True', 'T', '1']:
-                    overwritejp2 = True
-                elif arg in ['False', 'F', '0']:
-                    overwritejp2 = False
-                else:
-                    overwritejp2 = np.bool_(arg)
-            elif opt in ('-O', '--overwritefits'):
-                if arg in ['True', 'T', '1']:
-                    overwritefits = True
-                elif arg in ['False', 'F', '0']:
-                    overwritefits = False
-                else:
-                    overwritefits = np.bool_(arg)
-        nargs = len(args)
-        if nargs == 3:
-            year = int(args[0])
-            month = int(args[1])
-            day = int(args[2])
-        else:
-            year = None
-            month = None
-            day = None
-    except getopt.GetoptError as err:
-        print(err)
-        print('Error interpreting command line argument')
-        year = None
-        month = None
-        day = None
-        ndays = 1
-        clearcache = True
-        opts = []
-        overwritejp2 = False
-        overwritefits = False
-
-    # ##debug
-    # year = 2023
-    # month = 1
-    # day = 5
-    # ndays = 1
-    # clearcache = False
-    # overwritejp2 = True
-    # overwritefits = True
-
-    print("Running eovsa_fitsutils for date {}-{}-{}.".format(year, month, day))
-    kargs = {'ndays': ndays,
-             'clearcache': clearcache,
-             'overwritejp2': overwritejp2,
-             'overwritefits': overwritefits}
-    for k, v in kargs.items():
-        print(k, v)
-
-    main(year, month, day, ndays, overwritejp2=overwritejp2, overwritefits=overwritefits)
+    print(f"Running eovsa_fitsutils for date {dateobj.strftime('%Y-%m-%d')}.")
+    print("Arguments:")
+    print(f"  ndays: {args.ndays}")
+    print(f"  overwritejp2: {args.overwritejp2}")
+    print(f"  overwritefits: {args.overwritefits}")
+
+    # Call the main function with the parsed datetime object.
+    main(dateobj, args.ndays, overwritejp2=args.overwritejp2, overwritefits=args.overwritefits)

suncasa/eovsa/eovsa_flare_pipeline.py

@@ -68,6 +68,55 @@ def get_user_confirmation(prompt):
             print("Invalid input. Please enter 'y' or 'n'.")
 
 class FlareSelfCalib():
+    """
+    FlareSelfCalib provides a pipeline for self-calibration and imaging of EOVSA solar flare data.
+
+    This class handles the detection of flare times and locations, self-calibration of visibility data,
+    and final imaging steps for flare events, as well as renaming and moving the output files to designated
+    web directories.
+
+    :param vis: Full path for the input visibility data or a pre-processed flare-calibrated dataset, defaults to None
+    :type vis: str or object, optional
+    :param workpath: Working directory path for pipeline outputs, defaults to './'
+    :type workpath: str, optional
+    :param logfile: Path to the logfile. If not provided, a default log file is generated based on the current time, defaults to None
+    :type logfile: str, optional
+
+    :raises ValueError: If the input visibility (when provided as a string) does not exist.
+
+    Example:
+        from suncasa.eovsa import eovsa_flare_pipeline
+        from eovsapy.util import Time
+
+        trange_str = ['2023-12-14T16:54:00', '2023-12-14T17:10:00']
+        flare_id = '20231214170000'
+
+        # Initialize the pipeline with a given time range
+        trange = Time(trange_str)
+        fp = eovsa_flare_pipeline.FlareSelfCalib(vis=trange)
+
+        # Run the self-calibration and imaging pipeline
+        fp.slfcal_pipeline(doselfcal=True, doimaging=True)
+
+        # Rename and move the output files to the designated web directories
+        fp.rename_move_files(
+            flare_id=flare_id,
+            fitsdir_web_tp='/data1/eovsa/fits/flares/',
+            movdir_web_tp='/common/webplots/SynopticImg/eovsamedia/eovsa-browser/',
+            msdir_web_tp='/data1/eovsa/fits/flares/',
+            dorename_fits=True,
+            domove_fits=True,
+            dorename_mov=True,
+            domove_mov=True,
+            domove_ms=True,
+            dormworkdir=True,
+            docopy=True
+        )
+
+    :return: None
+    :rtype: None
+    """
+
     def __init__(self, vis=None, workpath='./', logfile=None):
         ##========================= initial setups =================================
         self.workpath = workpath