Merge branch 'develop' into 'master'

Preview cleaned after first archive run

Closes #5, #24, #34, #35, and #36

See merge request sdc/grisinv!22

grisinv/data_files/header_list.csv

-Observation Information,POINT_ID,OBSERVER,OBS_TRGT,IOBS,NMAPS,NSTEPS
+Observation Information,POINT_ID,OBSERVER,OBS_TRGT,IOBS,NMAPS,NSTEPS,RMAPS
 Temporal Information,DATE,DATE-OBS,DATE-BEG,DATE-END,DATE-AVG
 Instrument and Processing Info,ORIGIN,OBSRVTRY,TELESCOP,INSTRUME,GRATING,CAMERA,OBSGEO-X,OBSGEO-Y,OBSGEO-Z,AWAVLNTH,AWAVMAX,AWAVMIN,WAVEUNIT
 Detector Readout and Data Scaling,OBS_MODE,IMGSYS,XPOSURE,TEXPOSUR,NSUMEXP

grisinv/invert.py

@@ -17,9 +17,14 @@ from datetime import date, datetime
 import matplotlib.pyplot as plt
 import matplotlib.gridspec as gridspec
 import matplotlib.animation as animation
+from matplotlib.figure import Figure
 
 from scipy import interpolate
 
+from tempfile import TemporaryDirectory
+from glob import glob
+import subprocess
+
 #supress fits obsgeo warnings
 warnings.filterwarnings('ignore', category=FITSFixedWarning)
 
@@ -847,6 +852,9 @@ def add_l1_keys(hdr, vfisv_params):
             ]:
                 copy_keywords(vfisv_params.header, hdr, ikeys)
 
+    # Add RMAPS
+    hdr.comments["RMAPS"] = "Number of completed maps"
+
 
 def make_header(vfisv_params):
     """Create the header according to L1 keywords."""
@@ -882,8 +890,73 @@ def make_header(vfisv_params):
     return hdr
 
 
+class GrisFigure(Figure):
+    """A custom figure class for previews"""
+
+    def __init__(self, *args, **kwargs):
+        self.params = {
+            "data": kwargs.pop("data", []),
+            "header": kwargs.pop("header", [])
+        }
+        self.data = self.params["data"]
+        self.header = self.params["header"]
+        self.wcs = WCS(header)
+        super().__init__(*args, **kwargs)
+
+        aspect_ratio = self.data.shape[1] / self.data.shape[2]
+        shape_inches = (max([8 / aspect_ratio + 0.5, 3]) + 1.0, 8)
+
+        self.spec = gridspec.GridSpec(ncols=3, nrows=2, figure=self, width_ratios=[0.1, shape_inches[1], 0.1])
+
+
+    def show_plot(self):
+        self.pax = self.add_subplot(self.spec[:, 1], projection=self.wcs, slices=["x", "y", 0])
+        self.plot = self.pax.imshow(self.data.squeeze(), origin="lower")
+
+    def add_colorbar(self):
+        cax = self.add_subplot(self.spec[1, 2])
+        cbar = self.colorbar(self.plot, cax=cax)
+
+    def add_second_axis(self):
+        xdeltarcsec = 1
+        ydeltarcsec = 1
+        secax_x = self.pax.secondary_xaxis('top', functions=(lambda x: x * xdeltarcsec, lambda x: x / xdeltarcsec))
+        secax_y = self.pax.secondary_yaxis('right', functions=(lambda x: x * ydeltarcsec, lambda x: x / ydeltarcsec))
+        secax_x.set_xlabel("x (arcsec)")
+        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):
+    """Create individual maps"""
+    wcs = WCS(header)
+    aspect_ratio = data.shape[0] / data.shape[1]
+    shape_inches = (max([8 / aspect_ratio + 0.5, 3]) + 1.0, 8)
+    fig = plt.figure(figsize=shape_inches, dpi=dpi, constrained_layout=True)
+    spec = gridspec.GridSpec(nrows=2, ncols=2, figure=fig, width_ratios=[data.shape[1], min([5.,0.1*data.shape[1]/2.])], wspace=0.05, hspace=0.)
+    ax = fig.add_subplot(spec[:, 0], projection=wcs, slices=["x", "y", 0])
+    cax = fig.add_subplot(spec[1, 1])
+
+    cb = ax.imshow(
+        data.value, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
+    )
+
+    # ax.set_title(title)
+    fig.suptitle(title)
+
+    # set lon_lat axis
+    set_axis_lon_lat(ax, grid=grid)
+
+    # additonal axis
+    set_additional_axis(ax, header)
+
+    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')
+
+
 def plot_image(
-    data, header, preview, cbar_title=None, dpi=100, vmin=None, vmax=None, cmap=None, grid=None, cfmt=None,cticks=None
+    data, header, preview, cbar_title=None, vmin=None, vmax=None, cmap=None, grid=None, cfmt=None,cticks=None
 ):
     """
     Plot preview images of the inversion data based on the header params
@@ -906,7 +979,6 @@ def plot_image(
     """
 
     # TODO: Clean this. make it less cluttered.
-    wcs = WCS(header)
 
     title = (
         f"GRIS Inversions - {header['WAVE_STR']} {header['WAVELNTH']}"
@@ -914,122 +986,130 @@ def plot_image(
         + header["POINT_ID"]
     )
 
-    aspect_ratio = data.shape[1] / data.shape[2]
-    shape_inches = (max([8 / aspect_ratio + 0.5, 3]) + 1.0, 8)
-    fig = plt.figure(figsize=shape_inches, dpi=dpi, constrained_layout=True)
-    spec = gridspec.GridSpec(nrows=2, ncols=2, figure=fig, width_ratios=[data.shape[2],5.],wspace=0.1,hspace=0.)
-    ax = fig.add_subplot(spec[:, 0], projection=wcs, slices=["x", "y", 0])
-    cax = fig.add_subplot(spec[1, 1])
 
-    if header["NMAPS"] == 1:
 
-        unit_data, data_units, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,0)
+    if header["RMAPS"] == 1:
+
+        unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,0)
 
         # a quick correction for plots of field
         if data[0, ...].unit == Unit.Gauss:
             unit_vmax = (
-                np.max(unit_data)
-                if (unit_vmax != None) and (unit_vmax > np.max(unit_data))
+                np.max(unit_data.value)
+                if (unit_vmax != None) and (unit_vmax > np.max(unit_data.value))
                 else unit_vmax
             )
 
-        cb = ax.imshow(
-            unit_data, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
-        )
-
-        #ax.set_title(title)
-        fig.suptitle(title)
+        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)
 
-        # set lon_lat axis
-        set_axis_lon_lat(ax,grid=grid)
-
-        # additonal axis
-        set_additional_axis(ax, header)
-
-        cbar = plt.colorbar(cb, cax=cax,format=cfmt)
-        cbar.set_ticks(cticks)
-        cbar.ax.set_ylabel(cbar_title+ f" ({data_units})")
-        plt.savefig(preview,bbox_inches='tight')
     else:
+        # create a temp directory
+        td = TemporaryDirectory()
+
         ims = []
-        for imap in range(header["NMAPS"]):
+        for imap in range(header["RMAPS"]):
             # TODO: Clean this
-            unit_data, data_units, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,imap)
+            unit_data, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,imap)
 
             # a quick correction for plots of field
             if data[imap, ...].unit == Unit.Gauss:
                 unit_vmax = (
-                    np.max(unit_data)
-                    if (unit_vmax != None) and (unit_vmax > np.max(unit_data))
+                    np.max(unit_data.value)
+                    if (unit_vmax != None) and (unit_vmax > np.max(unit_data.value))
                     else unit_vmax
                 )
 
-            cb = ax.imshow(
-                unit_data,
-                origin="lower",
-                vmin=unit_vmin,
-                vmax=unit_vmax,
-                cmap=cmap,
-                animated=True,
-            )
-            if imap == 0:
-
-                # set lon_lat axis
-                set_axis_lon_lat(ax,grid=grid)
-
-                # additonal axis
-                set_additional_axis(ax,header)
+            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)
 
-                cbar = plt.colorbar(cb, cax=cax,format=cfmt)
-                cbar.set_ticks(cticks)
-                cbar.ax.set_ylabel(cbar_title+ f" ({data_units})")
-                fig.canvas.draw()
+        filelist = glob(os.path.join(td.name, "*"))
 
-            fig.suptitle(title)
+        merge_to_gif(filelist,preview.replace('.png','.gif'))
+        td.cleanup()
 
-            ims.append([cb])
+    return None
 
-        ani = animation.ArtistAnimation(fig, ims, interval=500.0, blit=True)
-        if ".png" in preview:
-            ani.save(preview.replace(".png", ".gif"))
-        elif ".jpeg" in preview:
-            ani.save(preview.replace(".jpeg", ".gif"))
-        elif ".jpg" in preview:
-            ani.save(preview.replace(".jpeg", ".gif"))
-        else:
-            ani.save(preview+".gif")
+def merge_to_gif(filelist, outfile):
+    """Convert a list of images to a single gif"""
+    nfiles = len(filelist)
+    delay = max(20, 100 / nfiles)
+    command = [
+        "convert",
+        "-delay",
+        f"{delay}",
+        "-loop",
+        "0",
+        *filelist,
+        outfile,
+    ]
+    subprocess.call(command)
 
-    return None
 
 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, ...].value
+        unit_data = data[imap, ...]
         unit_vmin = vmin
         unit_vmax = vmax
-        data_units = data.unit
     else:
     # if vmax is not provided with a unit.
-        unit_data = data[imap, ...].to(vmax.unit).value
-        data_units = data[imap, ...].to(vmax.unit).unit
+        unit_data = data[imap, ...].to(vmax.unit)
         unit_vmin = vmin.value
         unit_vmax = vmax.value
 
-    return unit_data, data_units, unit_vmin, unit_vmax
+    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
+
+        if spacing < 1:
+            grid=False
+        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"""
     lon = ax.coords['hpln']
     lat = ax.coords['hplt']
+
     lon.set_axislabel("HPLN-TAN")
     lat.set_axislabel("HPLT-TAN")
     lon.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
     lat.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
-    if grid:
+    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(
@@ -1051,6 +1131,7 @@ def set_axis_lon_lat(ax,grid=None):
         labelleft=True)
 
 
+
 def set_additional_axis(ax,header):
     """Format the x/y axes"""
     xdeltarcsec = attach_unit(abs(header['CDELT1']),Unit.degree,convert=Unit.arcsecond,value=True)
@@ -1442,9 +1523,19 @@ class VFISVpackage:
 
         self.ICONT = np.asfortranarray(self.ICONT)
 
-        # Change nmaps to reflect the fullfilled maps
+        # Truncate unfullfilled maps
+        self.truncate_maps()
+
+    def truncate_maps(self):
+        """Truncate unfullfileld maps"""
 
+        # Change nmaps to only the number of maps that were fullfilled
         self.header['NMAPS'] = max(self.imap).astype("int")
+        self.header['RMAPS'] = self.header['NMAPS']
+        # In case the steps were Change nmaps to reflect the fullfilled steps for the last map
+        if self.header['NMAPS'] > 1:
+            if self.header['NSTEPS'] > max(self.istep[self.imap == self.header['NMAPS']]).astype("int"):
+                self.header['RMAPS'] = self.header['NMAPS'] -1
 
     def read_fits(self, filename):
         """