Harrison 3063/#36 - PMCL/KJON - Maps spx: Added uncertainty for scalar coefficient from spectral fitting

Author: pleasant-menlo

imap_l3_processing/cdf/config/imap_hi_l3_variable_attrs.yaml

@@ -400,4 +400,29 @@ ena_spectral_scalar:
   LABL_PTR_1: energy_label
   LABL_PTR_2: longitude_label
   LABL_PTR_3: latitude_label
+  SCALETYP: log
+  DELTA_PLUS_VAR: ena_spectral_scalar_stat_uncert
+  DELTA_MINUS_VAR: ena_spectral_scalar_stat_uncert
+ena_spectral_scalar_stat_uncert:
+  NAME: ena_spectral_scalar_stat_uncert
+  DATA_TYPE: CDF_REAL4
+  CATDESC: ENA spectral index scalar coefficient statistical uncertainty.
+  DEPEND_0: epoch
+  DEPEND_1: energy
+  DEPEND_2: longitude
+  DEPEND_3: latitude
+  VAR_TYPE: data
+  RECORD_VARYING: RV
+  DISPLAY_TYPE: map_image
+  VARIABLE_PURPOSE: ' '
+  FIELDNAM: Spectral scalar stat unc
+  FORMAT: F12.6
+  LABLAXIS: Statistical unc.
+  UNITS: cm -2 s -1 sr -1 keV -1
+  VALIDMIN: 0.0
+  VALIDMAX: 1.0E+20
+  FILLVAL: -1.00E+31
+  LABL_PTR_1: energy_label
+  LABL_PTR_2: longitude_label
+  LABL_PTR_3: latitude_label
   SCALETYP: log

imap_l3_processing/cdf/config/imap_lo_l3_variable_attrs.yaml

@@ -663,4 +663,29 @@ ena_spectral_scalar:
   LABL_PTR_1: energy_label
   LABL_PTR_2: longitude_label
   LABL_PTR_3: latitude_label
+  SCALETYP: log
+  DELTA_PLUS_VAR: ena_spectral_scalar_stat_uncert
+  DELTA_MINUS_VAR: ena_spectral_scalar_stat_uncert
+ena_spectral_scalar_stat_uncert:
+  NAME: ena_spectral_scalar_stat_uncert
+  DATA_TYPE: CDF_REAL4
+  CATDESC: ENA spectral index scalar coefficient statistical uncertainty.
+  DEPEND_0: epoch
+  DEPEND_1: energy
+  DEPEND_2: longitude
+  DEPEND_3: latitude
+  VAR_TYPE: data
+  RECORD_VARYING: RV
+  DISPLAY_TYPE: map_image
+  VARIABLE_PURPOSE: ' '
+  FIELDNAM: Spectral scalar stat unc
+  FORMAT: F12.6
+  LABLAXIS: Statistical unc.
+  UNITS: cm -2 s -1 sr -1 keV -1
+  VALIDMIN: 0.0
+  VALIDMAX: 1.0E+20
+  FILLVAL: -1.00E+31
+  LABL_PTR_1: energy_label
+  LABL_PTR_2: longitude_label
+  LABL_PTR_3: latitude_label
   SCALETYP: log

imap_l3_processing/cdf/config/imap_ultra_l3_variable_attrs.yaml

@@ -400,4 +400,29 @@ ena_spectral_scalar:
   LABL_PTR_1: energy_label
   LABL_PTR_2: longitude_label
   LABL_PTR_3: latitude_label
+  SCALETYP: log
+  DELTA_PLUS_VAR: ena_spectral_scalar_stat_uncert
+  DELTA_MINUS_VAR: ena_spectral_scalar_stat_uncert
+ena_spectral_scalar_stat_uncert:
+  NAME: ena_spectral_scalar_stat_uncert
+  DATA_TYPE: CDF_REAL4
+  CATDESC: ENA spectral index scalar coefficient statistical uncertainty.
+  DEPEND_0: epoch
+  DEPEND_1: energy
+  DEPEND_2: longitude
+  DEPEND_3: latitude
+  VAR_TYPE: data
+  RECORD_VARYING: RV
+  DISPLAY_TYPE: map_image
+  VARIABLE_PURPOSE: ' '
+  FIELDNAM: Spectral scalar stat unc
+  FORMAT: F12.6
+  LABLAXIS: Statistical unc.
+  UNITS: cm -2 s -1 sr -1 keV -1
+  VALIDMIN: 0.0
+  VALIDMAX: 1.0E+20
+  FILLVAL: -1.00E+31
+  LABL_PTR_1: energy_label
+  LABL_PTR_2: longitude_label
+  LABL_PTR_3: latitude_label
   SCALETYP: log

imap_l3_processing/maps/map_models.py

@@ -38,6 +38,7 @@
 ENA_SPECTRAL_INDEX_VAR_NAME = "ena_spectral_index"
 ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME = "ena_spectral_index_stat_uncert"
 ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME = "ena_spectral_scalar"
+ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME = "ena_spectral_scalar_stat_uncert"
 
 ENA_INTENSITY_VAR_NAME = "ena_intensity"
 ENA_INTENSITY_STAT_UNCERT_VAR_NAME = "ena_intensity_stat_uncert"
@@ -123,6 +124,7 @@ class SpectralIndexMapData(MapData):
     ena_spectral_index: np.ndarray
     ena_spectral_index_stat_uncert: np.ndarray
     ena_spectral_index_scalar_coefficient: np.ndarray
+    ena_spectral_index_scalar_coefficient_stat_uncert: np.ndarray
 
 
 @dataclass
@@ -397,6 +399,8 @@ def to_healpix_skymap(self) -> HealpixSkyMap:
                 "ena_spectral_index_stat_uncert": (
                     full_shape, self.spectral_index_map_data.ena_spectral_index_stat_uncert),
                 "ena_spectral_index_scalar_coefficient": (full_shape, self.spectral_index_map_data.ena_spectral_index_scalar_coefficient),
+                "ena_spectral_index_scalar_coefficient_stat_uncert": (full_shape,
+                                                          self.spectral_index_map_data.ena_spectral_index_scalar_coefficient_stat_uncert),
             },
             coords={
                 CoordNames.TIME.value: self.spectral_index_map_data.epoch,
@@ -506,6 +510,7 @@ def _spectral_index_data_variables(data: SpectralIndexMapData) -> list[DataProdu
         DataProductVariable(ENA_SPECTRAL_INDEX_VAR_NAME, data.ena_spectral_index),
         DataProductVariable(ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME, data.ena_spectral_index_stat_uncert),
         DataProductVariable(ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, data.ena_spectral_index_scalar_coefficient),
+        DataProductVariable(ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, data.ena_spectral_index_scalar_coefficient_stat_uncert),
     ]
 
 

imap_l3_processing/maps/spectral_fit.py

@@ -32,6 +32,7 @@ def calculate_spectral_index_for_multiple_ranges(intensity_data: IntensityMapDat
         ena_spectral_index=np.concat([m.ena_spectral_index for m in spectral_maps], axis=1),
         ena_spectral_index_stat_uncert=np.concat([m.ena_spectral_index_stat_uncert for m in spectral_maps], axis=1),
         ena_spectral_index_scalar_coefficient=np.concat([m.ena_spectral_index_scalar_coefficient for m in spectral_maps], axis=1),
+        ena_spectral_index_scalar_coefficient_stat_uncert=np.concat([m.ena_spectral_index_scalar_coefficient_stat_uncert for m in spectral_maps], axis=1)
     )
 
 
@@ -81,7 +82,7 @@ def fit_spectral_index_map(intensity_data: IntensityMapData) -> SpectralIndexMap
     mean_energy = np.sqrt(min_energy * max_energy)
     new_energy_label = f"{min_energy} - {max_energy}"
 
-    output_scalar_coefficients, output_gammas, output_gamma_errors = fit_arrays_to_power_law(fluxes, uncertainty, energy)
+    output_scalar_coefficients, output_scalar_errors, output_gammas, output_gamma_errors = fit_arrays_to_power_law(fluxes, uncertainty, energy)
     mean_obs_date = calculate_datetime_weighted_average(intensity_data.obs_date,
                                                         weights=intensity_data.exposure_factor,
                                                         axis=1, keepdims=True)
@@ -110,12 +111,13 @@ def fit_spectral_index_map(intensity_data: IntensityMapData) -> SpectralIndexMap
         solid_angle=intensity_data.solid_angle,
         ena_spectral_index=output_gammas,
         ena_spectral_index_stat_uncert=output_gamma_errors,
-        ena_spectral_index_scalar_coefficient=output_scalar_coefficients
+        ena_spectral_index_scalar_coefficient=output_scalar_coefficients,
+        ena_spectral_index_scalar_coefficient_stat_uncert=output_scalar_errors,
     )
 
 
 def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energy: np.ndarray) -> tuple[
-    np.ndarray, np.ndarray, np.ndarray]:
+    np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
     par_info = [
         {'limits': [0.0, 1000.0]},
         {'limits': [0.0, 1000.0]},
@@ -125,6 +127,7 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
     output_gammas = np.full(output_shape, np.nan, dtype=float)
     output_gamma_errors = np.full_like(output_gammas, np.nan)
     output_scalar_coefficients = np.full_like(output_gammas, np.nan)
+    output_scalar_coefficients_errors = np.full_like(output_gammas, np.nan)
 
     for epoch in range(fluxes.shape[0]):
         intensity = fluxes[epoch].reshape((fluxes[epoch].shape[0], -1))
@@ -133,6 +136,7 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
         gammas = np.full(intensity.shape[-1], np.nan, dtype=float)
         gamma_errors = np.full_like(gammas, np.nan)
         scalar_coefficients = np.full(intensity.shape[-1], np.nan, dtype=float)
+        scalar_coefficient_errors = np.full(intensity.shape[-1], np.nan, dtype=float)
 
         for i in range(intensity.shape[-1]):
             flux = intensity[:, i]
@@ -158,10 +162,12 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
                     gammas[i] = gamma
                     scalar_coefficients[i] = a
                     gamma_errors[i] = gamma_error
+                    scalar_coefficient_errors[i] = a_error
         output_gammas[epoch, 0] = gammas.reshape(fluxes.shape[2:])
         output_gamma_errors[epoch, 0] = gamma_errors.reshape(fluxes.shape[2:])
         output_scalar_coefficients[epoch, 0] = scalar_coefficients.reshape(fluxes.shape[2:])
-    return output_scalar_coefficients, output_gammas, output_gamma_errors
+        output_scalar_coefficients_errors[epoch, 0] = scalar_coefficient_errors.reshape(fluxes.shape[2:])
+    return output_scalar_coefficients, output_scalar_coefficients_errors, output_gammas, output_gamma_errors
 
 
 def power_law(params, **kwargs):

tests/maps/test_map_models.py

@@ -52,6 +52,7 @@ def test_rectangular_spectral_index_to_data_product_variables(self):
                     ena_spectral_index=sentinel.ena_spectral_index,
                     ena_spectral_index_stat_uncert=sentinel.ena_spectral_index_stat_uncert,
                     ena_spectral_index_scalar_coefficient=sentinel.ena_spectral_index_scalar_coefficient,
+                    ena_spectral_index_scalar_coefficient_stat_uncert=sentinel.ena_spectral_index_scalar_coefficient_stat_uncert,
                 ),
                 coords=RectangularCoords(
                     latitude_delta=sentinel.latitude_delta,
@@ -82,6 +83,7 @@ def test_rectangular_spectral_index_to_data_product_variables(self):
             DataProductVariable(map_models.ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME,
                                 sentinel.ena_spectral_index_stat_uncert),
             DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient),
+            DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient_stat_uncert),
             DataProductVariable(map_models.LATITUDE_DELTA_VAR_NAME, sentinel.latitude_delta),
             DataProductVariable(map_models.LATITUDE_LABEL_VAR_NAME, sentinel.latitude_label),
             DataProductVariable(map_models.LONGITUDE_DELTA_VAR_NAME, sentinel.longitude_delta),
@@ -236,6 +238,7 @@ def test_healpix_spectral_index_to_data_product_variables(self):
                     ena_spectral_index=sentinel.ena_spectral_index,
                     ena_spectral_index_stat_uncert=sentinel.ena_spectral_index_stat_uncert,
                     ena_spectral_index_scalar_coefficient=sentinel.ena_spectral_index_scalar_coefficient,
+                    ena_spectral_index_scalar_coefficient_stat_uncert=sentinel.ena_spectral_index_scalar_coefficient_stat_uncert,
                 ),
                 coords=HealPixCoords(
                     pixel_index=sentinel.pixel_index,
@@ -264,6 +267,7 @@ def test_healpix_spectral_index_to_data_product_variables(self):
             DataProductVariable(map_models.ENA_SPECTRAL_INDEX_STAT_UNC_VAR_NAME,
                                 sentinel.ena_spectral_index_stat_uncert),
             DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient),
+            DataProductVariable(map_models.ENA_SPECTRAL_INDEX_SCALAR_COEFFICIENT_STAT_UNCERT_VAR_NAME, sentinel.ena_spectral_index_scalar_coefficient_stat_uncert),
             DataProductVariable(map_models.PIXEL_INDEX_VAR_NAME, sentinel.pixel_index),
             DataProductVariable(map_models.PIXEL_INDEX_LABEL_VAR_NAME, sentinel.pixel_index_label),
         ]
@@ -721,6 +725,7 @@ def test_healpix_spectral_index_map_data_to_skymap(self):
             ena_spectral_index=fake_data_per_energy_per_pixel * 2.2,
             ena_spectral_index_stat_uncert=fake_data_per_energy_per_pixel * 2.3,
             ena_spectral_index_scalar_coefficient=fake_data_per_energy_per_pixel * 2.4,
+            ena_spectral_index_scalar_coefficient_stat_uncert=fake_data_per_energy_per_pixel * 2.5,
         )
 
         healpix_spectral_index_map_data = HealPixSpectralIndexMapData(
@@ -755,9 +760,10 @@ def test_healpix_spectral_index_map_data_to_skymap(self):
         np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_stat_uncert"].values, spectral_index_map_data.ena_spectral_index_stat_uncert)
 
         np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_scalar_coefficient"].values, spectral_index_map_data.ena_spectral_index_scalar_coefficient)
+        np.testing.assert_array_equal(actual_dataset.data_vars["ena_spectral_index_scalar_coefficient_stat_uncert"].values, spectral_index_map_data.ena_spectral_index_scalar_coefficient_stat_uncert)
 
         for key in [ "obs_date", "obs_date_range", "exposure_factor", "ena_spectral_index", "ena_spectral_index_stat_uncert",
-                     "ena_spectral_index_scalar_coefficient" ]:
+                     "ena_spectral_index_scalar_coefficient", "ena_spectral_index_scalar_coefficient_stat_uncert"  ]:
             self.assertEqual((CoordNames.TIME.value, CoordNames.ENERGY_L2.value, CoordNames.GENERIC_PIXEL.value), actual_dataset.data_vars[key].dims)
         for key in [ "latitude", "longitude", "solid_angle" ]:
             self.assertEqual((CoordNames.GENERIC_PIXEL.value,), actual_dataset.data_vars[key].dims)