Merge branch 'develop' into 'master'

Updates to preview images and logs

Closes #41, #39, #38, and #37

See merge request sdc/grisinv!28

README.md

@@ -128,6 +128,7 @@ mpiexec -n 1 vfisv \
   --preview='test_preview.png'  \
   --errors='test_errors.fits' \
   --diagnose='test_diagnose.fits' \
+  --log='test_log.txt'
 ```
 
 ### Python interface

grisinv/data_files/line_names.csv

 Fe I,15648.514,15648.510,3.0
+Fe I,15662.017,15662.014,1.5
 Si I,10827.091,10827.089,1.5
 Ca I,10838.970,10838.970,1.5
\ No newline at end of file

grisinv/data_files/vfisv_info.txt

+Very Fast Inversion of the Stokes Vector version 5.0
+Node for Spectrograph data: VFISV_SPEC
+J.M.Borrero | Leibniz Institut fuer Sonnenphysik.
+This software is distributed under GPL 2.0 license.
+You can use/modifiy/distribute as you wish as long
+as it is under the GPL 2.0 license. Please cite:
+Borrero et al. (2011) Sol.Phys. 273, 267 in your work.

grisinv/invert.py

@@ -926,7 +926,7 @@ class GrisFigure(Figure):
         secax_y.set_ylabel("y (arcsec)")
 
 
-def create_map(data,header,filename,unit_vmin=0,unit_vmax=1,cmap='Gray',grid='white',cfmt=None,title=None,dpi=100,cticks=None,cbar_title=None):
+def create_map(data,header,filename,unit_vmin=0,unit_vmax=1,cmap='Gray',grid='white',cfmt=None,title=None,dpi=100,cticks=None,cbar_title=None,annotate_starttime=False,annotate_endtime=False):
     """Create individual maps"""
     wcs = WCS(header)
     aspect_ratio = data.shape[0] / data.shape[1]
@@ -940,6 +940,8 @@ def create_map(data,header,filename,unit_vmin=0,unit_vmax=1,cmap='Gray',grid='wh
         data.value, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
     )
 
+
+
     # ax.set_title(title)
     fig.suptitle(title)
 
@@ -949,10 +951,23 @@ def create_map(data,header,filename,unit_vmin=0,unit_vmax=1,cmap='Gray',grid='wh
     # additonal axis
     set_additional_axis(ax, header)
 
+    # For now annotate time for maps
+    if annotate_starttime:
+        start_time = datetime.fromisoformat(header["DATE-BEG"]).strftime("%H:%M:%S")
+        ax.text(0.0, 1., f'{start_time} |', ha='right', va='center', fontsize=8,
+                color='gray', transform=ax.transAxes, bbox=dict(facecolor='white', alpha=1., edgecolor='None', pad=1.))
+
+    if annotate_endtime:
+        end_time = datetime.fromisoformat(header["DATE-END"]).strftime("%H:%M:%S")
+        ax.text(1.0, 1., f'| {end_time}', ha='left', va='center', fontsize=8,
+                color='gray', transform=ax.transAxes, bbox=dict(facecolor='white', alpha=1.,edgecolor='None',pad=0.5))
+
+
     cbar = plt.colorbar(cb, cax=cax, format=cfmt)
     cbar.set_ticks(cticks)
     cbar.ax.set_ylabel(cbar_title + f" ({data.unit})")
     plt.savefig(filename, bbox_inches='tight')
+    plt.close('all')
 
 
 def plot_image(
@@ -986,11 +1001,9 @@ def plot_image(
         + header["POINT_ID"]
     )
 
-
-
     if header["RMAPS"] == 1:
 
-        unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,0)
+        unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax)
 
         # a quick correction for plots of field
         if data[0, ...].unit == Unit.Gauss:
@@ -999,17 +1012,31 @@ def plot_image(
                 if (unit_vmax != None) and (unit_vmax > np.max(unit_data.value))
                 else unit_vmax
             )
+            if unit_vmax > 1.:
+                unit_vmax = np.around(unit_vmax)
+                cticks = None
+            else:
+                unit_vmax = (np.mean(unit_data.value) + np.max(unit_data.value)) / 2.
 
-        create_map(unit_data, header, preview, unit_vmin=unit_vmin, unit_vmax=unit_vmax, cmap=cmap, grid=grid, cfmt=cfmt, title=title,cticks=cticks,cbar_title=cbar_title)
+                unit_vmax = np.around(unit_vmax, decimals=1) * 1000
+                unit_data = unit_data.to(Unit.Gauss)
+                cticks = None
+                cfmt = '%3d'
+
+        create_map(unit_data[0,...], header, preview, unit_vmin=unit_vmin, unit_vmax=unit_vmax, cmap=cmap, grid=grid, cfmt=cfmt, title=title,cticks=cticks,cbar_title=cbar_title,annotate_starttime=True,annotate_endtime=True)
 
     else:
         # create a temp directory
         td = TemporaryDirectory()
-
-        ims = []
-        for imap in range(header["RMAPS"]):
+        rmaps = header["RMAPS"]
+        for imap in tqdm(range(rmaps),
+                        bar_format=termcolors.text
+                        + "Creating animation   {n_fmt}/{total_fmt}: {desc}|{bar:30}|"
+                        + termcolors.end
+                        + "{percentage:3.0f}%  ({elapsed_s:3.1f}s)"
+                        ):
             # TODO: Clean this
-            unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,imap)
+            unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax)
 
             # a quick correction for plots of field
             if data[imap, ...].unit == Unit.Gauss:
@@ -1018,10 +1045,27 @@ def plot_image(
                     if (unit_vmax != None) and (unit_vmax > np.max(unit_data.value))
                     else unit_vmax
                 )
+                if unit_vmax>1.:
+                    unit_vmax = np.around(unit_vmax)
+                    cticks = None
+                else:
+                    unit_vmax = (np.mean(unit_data.value)+np.max(unit_data.value))/2.
+
+                    unit_vmax = np.around(unit_vmax, decimals=1)*1000
+                    unit_data = unit_data.to(Unit.Gauss)
+                    cticks = None
+                    cfmt   = '%3d'
 
             outpath = os.path.join(td.name, f'{imap:03d}_'+preview)
-            create_map(unit_data, header, outpath, unit_vmin=unit_vmin, unit_vmax=unit_vmax, cmap=cmap, grid=grid, cfmt=cfmt,
-                       title=title+f'_{imap:03d}',cticks=cticks,cbar_title=cbar_title)
+            annotate_starttime=False
+            annotate_endtime=False
+            if imap==0:
+                annotate_starttime = True
+            if imap==header["RMAPS"]-1:
+                annotate_endtime = True
+
+            create_map(unit_data[imap,...], header, outpath, unit_vmin=unit_vmin, unit_vmax=unit_vmax, cmap=cmap, grid=grid, cfmt=cfmt,
+                       title=title+f'_{imap:03d}',cticks=cticks,cbar_title=cbar_title,annotate_starttime=annotate_starttime,annotate_endtime=annotate_endtime)
 
         filelist = glob(os.path.join(td.name, "*"))
 
@@ -1029,6 +1073,22 @@ def plot_image(
         td.cleanup()
 
     return None
+def convert_units(data,vmin,vmax):
+    '''Convert the units based on vmin/vmax'''
+
+    if vmax == None:
+    # if vmax is not provided with a unit
+        unit_data = data
+        unit_vmin = vmin
+        unit_vmax = vmax
+    else:
+    # if vmax is provided with a unit.
+        unit_data = data.to(vmax.unit)
+        unit_vmin = vmin.value
+        unit_vmax = vmax.value
+
+    return unit_data, unit_vmin, unit_vmax
+
 
 def merge_to_gif(filelist, outfile):
     """Convert a list of images to a single gif"""
@@ -1037,107 +1097,78 @@ def merge_to_gif(filelist, outfile):
     command = [
         "convert",
         "-delay",
+        f"{delay*2}",
+        filelist[0],
+        "-delay",
         f"{delay}",
+        *filelist,
         "-loop",
         "0",
-        *filelist,
         outfile,
     ]
     subprocess.call(command)
 
 
-def convert_units(data,vmin,vmax,imap):
-    '''Convert the units based on vmin/vmax'''
-
-    if vmax == None:
-    # if vmax is not provided
-        unit_data = data[imap, ...]
-        unit_vmin = vmin
-        unit_vmax = vmax
-    else:
-    # if vmax is not provided with a unit.
-        unit_data = data[imap, ...].to(vmax.unit)
-        unit_vmin = vmin.value
-        unit_vmax = vmax.value
-
-    return unit_data, unit_vmin, unit_vmax
-
-def change_grid_spacing(axis,direction,canvas):
-    """Check for labels and and change the spacing"""
-
-    spacing=4.
-    divider=2.
-    grid=True
 
-    while direction not in axis.ticklabels.text:
-        plt.rcParams.update({'figure.max_open_warning': 0})
-        axis.set_ticks(spacing=spacing * Unit.arcsec, color='dimgray')
-        spacing /= divider
+def switch_coordinates(lon,lat,canvas):
+    """in case the the labels are nor proper, switch the axis"""
 
-        if spacing < 1:
-            grid=False
+    # Check if the left or the bottom labels are drwan.
+    # If not, then this is possibly due to slit aligned horizontally.
+    if ("l" not in lat.ticklabels.text) or ("b" not in lon.ticklabels.text):
+        lon.set_ticks_position('l')
+        lon.set_ticklabel_position('l')
+        lon.set_axislabel_position('l')
+        lat.set_ticks_position('b')
+        lat.set_ticklabel_position('b')
+        lat.set_axislabel_position('b')
+        canvas.draw()
+    else:
+        lat.set_ticks_position('l')
+        lat.set_ticklabel_position('l')
+        lat.set_axislabel_position('l')
+        lon.set_ticks_position('b')
+        lon.set_ticklabel_position('b')
+        lon.set_axislabel_position('b')
         canvas.draw()
-        if spacing < 0.125:
-            grid = False
-            return grid
-
-    return grid
-
-
 
 def set_axis_lon_lat(ax,grid=None):
-    """Format the Longitude and Latitide axes"""
+    """Format the Longitude and Latitude axes"""
+    # get the longitude and latitude
     lon = ax.coords['hpln']
     lat = ax.coords['hplt']
 
+    # label the longitude and latitude
     lon.set_axislabel("HPLN-TAN")
     lat.set_axislabel("HPLT-TAN")
+
+    # set the ticks for the longitude and latitude
     lon.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
     lat.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
+
+    # draw the figure
     canvas = ax.figure.canvas
     canvas.draw()
-    latgrid = change_grid_spacing(lat, "l", canvas)
-    longrid = change_grid_spacing(lon, "b", canvas)
-    # if 'l' not in lat.ticklabels.text:
-    #     print('left not present')
-    #     lat.set_ticks(spacing=2 * Unit.arcsec, color='dimgray')
-    #     canvas.draw()
-    # if 'b' not in lon.ticklabels.text:
-    #     print('bottom not present')
-    #     lon.set_ticks(spacing=2 * Unit.arcsec, color='dimgray')
-    #     canvas.draw()
-
-    if longrid:
-        lon.grid(color=grid, alpha=0.5, linestyle='dashed')
-    if latgrid:
-        lat.grid(color=grid, alpha=0.5, linestyle='dashed')
-
-    lon.tick_params(
-        axis='both',
-        which='both',
-        bottom=True,
-        top=False,
-        right=False,
-        left=False,
-        labelbottom=True)
-
-    lat.tick_params(
-        axis='both',
-        which='both',
-        bottom=False,
-        top=False,
-        right=False,
-        left=True,
-        labelleft=True)
 
+    # Switch axis if the labels are not created.
+    # in cases where the slits are horionzatally-oriented.
+    switch_coordinates(lon, lat, canvas)
 
+    # Add gris lines for longitude and latitude
+    lon.grid(color=grid, alpha=0.5, linestyle='dashed')
+    lat.grid(color=grid, alpha=0.5, linestyle='dashed')
 
 def set_additional_axis(ax,header):
     """Format the x/y axes"""
+    # Get the delta_x and delta_y
     xdeltarcsec = attach_unit(abs(header['CDELT1']),Unit.degree,convert=Unit.arcsecond,value=True)
     ydeltarcsec = attach_unit(abs(header['CDELT2']), Unit.degree, convert=Unit.arcsecond, value=True)
+
+    # Append the x/y axes
     secax_x = ax.secondary_xaxis('top', functions=(lambda x: x * xdeltarcsec, lambda x: x / xdeltarcsec))
     secax_y = ax.secondary_yaxis('right', functions=(lambda x: x * ydeltarcsec, lambda x: x / ydeltarcsec))
+
+    # set labels
     secax_x.set_xlabel("x (arcsec)")
     secax_y.set_ylabel("y (arcsec)")
 
@@ -1159,6 +1190,8 @@ def attach_unit(data,units,convert=None,value=False):
 
 def create_preview(data, header, preview):
     """Create preview image"""
+
+    # Loop through the default quantities to plot
     for iquant in vfisv_previews:
         if iquant=='b_para':
             quant_data = np.abs(data['bmag']*np.cos(data['incl']))
@@ -1173,6 +1206,7 @@ def create_preview(data, header, preview):
             mag_mask = data["bmag"].value < 1000.
             quant_data -= np.mean(quant_data[mag_mask])
 
+        # plot the quantity
         plot_image(
             quant_data,
             header,
@@ -1188,6 +1222,35 @@ def create_preview(data, header, preview):
     return None
 
 
+def write_log_file(log,header,data,stokes,total_time,numproc):
+
+    vfisv_info = (
+        pkgutil.get_data("grisinv", os.path.join("data_files", "vfisv_info.txt"))
+        .decode()
+    )
+    f = open(log, "w")
+    f.write("SDC:GRIS Inversion\n")
+    f.write("------------------\n\n")
+    f.write("The inversion was carried out using:\n\n")
+    f.write(vfisv_info)
+    f.write("\nDetails of the run:\n\n")
+
+    for run_args in header["PRPARA1"].split(','):
+        argu,value = run_args.split('=')
+
+        f.write(f"{argu.capitalize()}\t: {value}\n")
+    f.write("\n")
+    f.write(f"Total time for inversion : {total_time:6.2f} s\n")
+    f.write(f"Number of processors     : {numproc}\n")
+    f.write(f"Total profiles inverted  : {data['vlos'].size}\n\n")
+
+    f.write(f"Mean  chi_sqr  : {data['chsq'].value.mean():5.2f}\n")
+    f.write(f"Best  chi_sqr  : {data['chsq'].value.min():5.2f}\n")
+    f.write(f"Worst chi_sqr  : {data['chsq'].value.max():5.2f}\n")
+
+    f.close()
+
+
 def copy_keywords(from_header, to_header, keyword):
     """Copy keyword and comments from one header to another"""
     try:
@@ -1285,6 +1348,7 @@ class VFISVpackage:
         self.map_props = None
         self.imap = None
         self.istep = None
+        self.time = None
         self.date_obs = None
         self.date_beg = None
         self.date_end = None
@@ -1458,7 +1522,7 @@ class VFISVpackage:
 
 
             # Set the pix continnum to pix_end
-            self.pix_cont = self.pix_end
+            self.pix_cont = self.pix_ini
 
             # Set the start of the wavelenght array for
             self.wave_start = wave_array[self.pix_ini]
@@ -1496,8 +1560,9 @@ class VFISVpackage:
         self.SV = np.zeros_like(self.SI)
         self.ICONT = np.zeros([len(self.filenames), nsteps], dtype="float32")
 
-        self.imap = np.zeros(len(self.filenames))
-        self.istep = np.zeros(len(self.filenames))
+        self.imap = np.zeros(len(self.filenames), dtype="int")
+        self.istep = np.zeros(len(self.filenames), dtype="int")
+        self.time = np.zeros(len(self.filenames), dtype='datetime64[s]')
 
         for iff in tqdm(
             range(len(self.filenames)),
@@ -1514,6 +1579,7 @@ class VFISVpackage:
                 self.ICONT[iff, :],
                 self.imap[iff],
                 self.istep[iff],
+                self.time[iff]
             ) = self.read_fits(self.filenames[iff])
 
         self.SI = np.asfortranarray(self.SI)
@@ -1559,6 +1625,7 @@ class VFISVpackage:
             sq = data[1, self.pix_ini : self.pix_end + 1, :, 0].copy()
             su = data[2, self.pix_ini : self.pix_end + 1, :, 0].copy()
             sv = data[3, self.pix_ini : self.pix_end + 1, :, 0].copy()
+
             ic = data[0, self.pix_cont, :, 0].copy()
 
         return (
@@ -1569,6 +1636,7 @@ class VFISVpackage:
             ic.T,
             hdul[0].header["IMAP"],
             hdul[0].header["ISTEP"],
+            np.datetime64(hdul[0].header["DATE-OBS"]),
         )
 
     def check_fits(self, header):
@@ -1599,7 +1667,8 @@ class VFISVpackage:
 @click.option("--preview", type=str, help="Filename to save the plot")
 @click.option("--errors", type=str, help="Filename to save the uncertainities")
 @click.option("--diagnose", type=str, help="Filename to write the inversion fits")
-def main(path, id, out, numproc, line, width, preview, errors, diagnose):
+@click.option("--log", type=str, help="Filename to write the log file")
+def main(path, id, out, numproc, line, width, preview, errors, diagnose,log):
     """SDC:GRIS Inversion Pipeline code"""
 
     # Get the MPI Communicator
@@ -1649,6 +1718,13 @@ def main(path, id, out, numproc, line, width, preview, errors, diagnose):
         if diagnose:
             write_stokes_fits(stokes, header, diagnose)
 
+        # write log file
+        if log:
+            write_log_file(log,header,data,stokes,total_time,numproc)
+
+
+
+
 
 if __name__ == "__main__":
     main()