Adding more metadata to fits files and EISMap

eispac/core/eiscube.py

@@ -4,6 +4,7 @@
 import copy
 import numpy as np
 import astropy.units as u
+from astropy.nddata import StdDevUncertainty
 from astropy.convolution import convolve, CustomKernel
 from astropy.coordinates import SkyCoord
 from ndcube import __version__ as ndcube_ver
@@ -228,7 +229,7 @@ def apply_radcal(self, input_radcal=None):
                 new_radcal = new_radcal/self.radcal
 
         new_data = self.data.copy()*new_radcal
-        new_errs = self.uncertainty.array.copy()*new_radcal
+        new_errs = StdDevUncertainty(self.uncertainty.array.copy()*new_radcal)
         new_meta = copy.deepcopy(self.meta)
         new_meta['mod_index']['bunit'] = 'erg / (cm2 s sr)'
         new_meta['notes'].append('Applied radcal to convert photon counts to intensity')
@@ -260,7 +261,7 @@ def remove_radcal(self):
             return self
 
         new_data = self.data.copy()/self.radcal
-        new_errs = self.uncertainty.array.copy()/self.radcal
+        new_errs = StdDevUncertainty(self.uncertainty.array.copy()/self.radcal)
         new_meta = copy.deepcopy(self.meta)
         new_meta['mod_index']['bunit'] = 'photon'
         new_meta['notes'].append('Removed radcal to convert intensity to photon counts')
@@ -471,6 +472,7 @@ def smooth_cube(self, width=3, **kwargs):
         if self.uncertainty is not None:
             sm_errs = np.sqrt(convolve(self.uncertainty.array**2,
                                        sm_kernel, **kwargs))
+            sm_errs = StdDevUncertainty(sm_errs)
         else:
             sm_errs = None
         sm_data_mask = np.logical_or(np.isnan(sm_data), sm_data < 0)

eispac/core/eisfitresult.py

@@ -6,6 +6,7 @@
 from scipy.ndimage import shift as shift_img
 from scipy.interpolate import interp1d
 import matplotlib.pyplot as plt
+from astropy.nddata import StdDevUncertainty
 
 from eispac import __version__ as eispac_version
 import eispac.core.fitting_functions as fit_fns
@@ -471,16 +472,18 @@ def get_map(self, component=0, measurement='intensity', **kwargs):
         if self.meta is None or 'mod_index' not in self.meta.keys():
             print("Error: Missing mod_index containing pointing information.")
             return None
+
+        if 'date_obs' not in self.meta.keys():
+            self.meta['date_obs'] = None
+
+        if 'duration' not in self.meta.keys():
+            self.meta['duration'] = None
 
         # Fetch index from the meta structure, cut out spectral data and update
         hdr_dict = copy.deepcopy(self.meta['mod_index'])
-        void = hdr_dict.pop('cname3', None)
-        void = hdr_dict.pop('crval3', None)
-        void = hdr_dict.pop('crpix3', None)
-        void = hdr_dict.pop('cdelt3', None)
-        void = hdr_dict.pop('ctype3', None)
-        void = hdr_dict.pop('cunit3', None)
-        void = hdr_dict.pop('naxis3', None)
+        for KEY in ['cname3', 'crval3', 'crpix3', 
+                    'cdelt3', 'ctype3', 'cunit3', 'naxis3']:
+            void = hdr_dict.pop(KEY, None)
         hdr_dict['naxis'] = 2
         hdr_dict['line_id'] = self.fit['line_ids'][gauss_ind]
         code_ver = self.eispac_version
@@ -507,7 +510,13 @@ def get_map(self, component=0, measurement='intensity', **kwargs):
                  +' with a value of "intensity", "velocity", or "width"')
             return None
 
-        output_map = EISMap(data_array, hdr_dict, uncertainty=err_array, **kwargs)
+        date_obs_arr = self.meta['date_obs']
+        exptime_arr = self.meta['duration']
+        output_map = EISMap(data_array, hdr_dict, 
+                            uncertainty=StdDevUncertainty(err_array), 
+                            step_date_obs=date_obs_arr, 
+                            step_exptime=exptime_arr, 
+                            **kwargs)
 
         return output_map
 

eispac/core/eismap.py

@@ -1,11 +1,13 @@
 """
 `~sunpy.map.Map` subclass for the EUV Imaging Spectrometer (EIS) on Hinode
 """
+import sys
 import pathlib
 
 import numpy as np
 import astropy.units as u
 from astropy.io import fits
+from astropy.nddata import StdDevUncertainty
 from astropy.visualization import ImageNormalize, AsinhStretch, LinearStretch
 import sunpy.map
 from sunpy.map.mapbase import SpatialPair
@@ -44,6 +46,9 @@ class EISMap(sunpy.map.GenericMap):
     """
 
     def __init__(self, data, header=None, **kwargs):
+        step_date_obs = None
+        step_exptime = None
+
         # Check for a fits file containing a binary table with the errs
         # NB: errs are in a table with y-axis number of rows, each with x-axis
         #     number of values
@@ -53,24 +58,60 @@ def __init__(self, data, header=None, **kwargs):
                     data = hdul[0].data
                     header = hdul[0].header
                     if len(hdul) >= 2:
-                        kwargs['uncertainty'] = hdul[1].data['errors']
-
+                        data_errors = StdDevUncertainty(hdul[1].data['errors'])
+                        kwargs['uncertainty'] = data_errors
+                    if len(hdul) >= 3:
+                        step_date_obs = parse_time(hdul[2].data['step_date_obs'])
+                        step_exptime = hdul[2].data['step_exptime']
+
+        # Get user-input step_date_obs and step_exptime
+        # NB: this will overwrite any values from an input fits file
+        user_step_date_obs = kwargs.pop('step_date_obs', None)
+        user_step_exptime = kwargs.pop('step_exptime', None)
+        if user_step_date_obs is not None and header is not None:
+            if len(user_step_date_obs) == header['naxis1']:
+                step_date_obs = parse_time(user_step_date_obs)
+            else:
+                print(f'WARNING incorrect number of "step_date_obs" values!'
+                     +f' This EIS observation has {header["naxis1"]} steps.',
+                      file=sys.stderr)
+        if user_step_exptime is not None and header is not None:
+            if len(user_step_exptime) == header['naxis1']:
+                step_exptime = user_step_exptime
+            else:
+                print(f'WARNING: incorrect number of "step_exptime" values!'
+                     +f' This EIS observation has {header["naxis1"]} steps.',
+                      file=sys.stderr)
+
+        # Initalize the map
         super().__init__(data, header, **kwargs)
+        self._step_date_obs = step_date_obs
+        self._step_exptime = step_exptime
 
-        # Setup plot settings
+        # Setup plot settings and get default data masks
+        # This includes adjusting colormap and normalization depending on whether 
+        # the map contains intensity, velocity, or line width data
+        default_mask = None
         self.plot_settings['aspect'] = self.meta['cdelt2'] / self.meta['cdelt1']
-        # Adjust colormap and normalization depending on whether the map contains
-        # intensity, velocity, or line width data
         if self.meta['measrmnt'].lower().startswith('int'):
             self.plot_settings['cmap'] = 'Blues_r'
             self.plot_settings['norm'] = ImageNormalize(stretch=AsinhStretch())
+            default_mask = self.data == 0
         elif self.meta['measrmnt'].lower().startswith('vel'):
             self.plot_settings['cmap'] = 'RdBu_r'
             # Autoscale color range to 3*std (rounded to nearest multiple of 5)
             vlim = 5*round(3*self.data.std()/5)
             self.plot_settings['norm'] = ImageNormalize(vmin=-vlim, vmax=vlim)
+            if self.uncertainty is not None:
+                # Note: velocities of 0 may be valid UNLESS the errors are NaN
+                default_mask = (self.data == 0) & np.isnan(self.uncertainty.array)
         elif self.meta['measrmnt'].lower().startswith('wid'):
             self.plot_settings['cmap'] = 'viridis'
+            default_mask = self.data == 0
+
+        # Set the default mask (ignored if the user input their own mask)
+        if self.mask is None:
+            self.mask = default_mask
 
     @property
     def spatial_units(self):
@@ -79,7 +120,7 @@ def spatial_units(self):
 
     @property
     def processing_level(self):
-        return self.meta.get('lvl_num', 3)
+        return self.meta.get('lvl_num', 2)
 
     @property
     def waveunit(self):
@@ -141,6 +182,54 @@ def reference_date(self):
         """
         return self.date_start
 
+    @property
+    def duration(self):
+        """Total duration of the observation in units of [min]
+        """
+        total_time = (np.datetime64(self.meta['date_end']) 
+                    - np.datetime64(self.meta['date_obs']))
+        total_time = total_time / np.timedelta64(60, 's') # convert to [min] 
+        return total_time * u.Unit('min')
+
+    @property
+    def step_date_obs(self):
+        """date_obs timestamp of each step along the x-axis
+        """
+        if self._step_date_obs is not None:
+            return self._step_date_obs
+        else:
+            print(f'WARNING: exact "step_date_obs" values are unknown!' 
+                 +f' Estimating based on the observation start and end times.',
+                  file=sys.stderr)
+            total_time = self.duration.to('s').value
+            est_cad = total_time / self.meta['naxis1']
+            est_date_obs = (np.datetime64(self.meta['date_obs']) 
+                           + np.arange(self.meta['naxis1'])
+                           * np.timedelta64(int(est_cad*1000), 'ms'))
+
+            if self.meta['nraster'] == 1: 
+                # Sit-and-stare timestamps inc left to right
+                return parse_time(est_date_obs)
+            else:
+                # Normal raster timestamps inc from right to left (scan dir)
+                return parse_time(np.flip(est_date_obs))
+
+    @property
+    def step_exptime(self):
+        """Exposure time of each step along the x-axis
+        """
+        if self._step_exptime is not None:
+            return self._step_exptime
+        else:
+            print(f'WARNING: exact "step_exptime" values are unknown!' 
+                 +f' Estimating based on the observation start and end times.'
+                 +f' Actual exposure times will be shorter due to on-board'
+                 +f' processing and the specific observation plan.',
+                  file=sys.stderr)
+            total_time = self.duration.to('s').value
+            est_avg_exptime = total_time / self.meta['naxis1']
+            return np.zeros(self.meta['naxis1']) + est_avg_exptime
+
     @classmethod
     def is_datasource_for(cls, data, header, **kwargs):
         """

eispac/core/export_fits.py

@@ -84,13 +84,9 @@ def export_fits(fit_result, save_dir=None, verbose=False):
 
     # Fetch index from the meta structure, cut out spectral data and update
     hdr_dict = copy.deepcopy(fit_result.meta['mod_index'])
-    void = hdr_dict.pop('cname3', None)
-    void = hdr_dict.pop('crval3', None)
-    void = hdr_dict.pop('crpix3', None)
-    void = hdr_dict.pop('cdelt3', None)
-    void = hdr_dict.pop('ctype3', None)
-    void = hdr_dict.pop('cunit3', None)
-    void = hdr_dict.pop('naxis3', None)
+    for KEY in ['cname3', 'crval3', 'crpix3', 
+                'cdelt3', 'ctype3', 'cunit3', 'naxis3']:
+        void = hdr_dict.pop(KEY, None)
     hdr_dict['naxis'] = 2
     code_ver = fit_result.eispac_version
     date_fit = fit_result.date_fit
@@ -112,7 +108,21 @@ def export_fits(fit_result, save_dir=None, verbose=False):
     err_hdr['tdim1'] = err_tdim
     err_hdr['tform1'] = err_tform
 
-    # Loop over all of the parameters and crewate the fits files
+    # Create the bintable table and header for step_date_obs and step_exptime 
+    col_date_obs = fits.Column(name='step_date_obs', format='24A', coord_type='UTC',
+                              array=fit_result.meta['date_obs'].astype('<U24'))
+    col_exptime = fits.Column(name='step_exptime', format='E', unit='s',
+                              array=fit_result.meta['duration'])
+    extra_hdu = fits.BinTableHDU.from_columns([col_date_obs, col_exptime])
+
+    # Add duplicate timestamps (to conform with the FITS-4 standard)
+    hdr_dict['date-obs'] = hdr_dict['date_obs']
+    hdr_dict['date-beg'] = hdr_dict['date_beg']
+    hdr_dict['date-avg'] = hdr_dict['date_avg']
+    hdr_dict['date-end'] = hdr_dict['date_end']
+    data_hdr = fits.Header(hdr_dict) # update primary header ONLY
+
+    # Loop over all of the parameters and create the fits files
     # If there are multiple line_ids, output each line to its own set of files
     # and return a list of filepaths (consistent with old IDL workflow)
     params = ['int', 'vel', 'wid']
@@ -162,7 +172,7 @@ def export_fits(fit_result, save_dir=None, verbose=False):
             err_col = fits.Column(name='errors', format=err_tform, dim=err_tdim,
                                   unit=err_hdr['tunit1'], array=err_array)
             err_hdu = fits.BinTableHDU.from_columns([err_col], header=err_hdr)
-            hdu_list = fits.HDUList([main_hdu, err_hdu])
+            hdu_list = fits.HDUList([main_hdu, err_hdu, extra_hdu])
             # hdu_list.writeto(output_files[-1], output_verify='silentfix',
             hdu_list.writeto(output_files[-1], output_verify='fix',
                              overwrite=True)