Harrison 3098 - EYEA - Maps: Use default SP values of 1 where there is no GLOWS data

Author: Menlo Innovations - CAVA Project

imap_l3_processing/maps/rectangular_survival_probability.py

@@ -1,3 +1,5 @@
+from typing import Optional
+
 import numpy as np
 import xarray as xr
 from imap_processing.ena_maps.ena_maps import RectangularSkyMap, PointingSet
@@ -23,7 +25,7 @@ def interpolate_angular_data_to_nearest_neighbor(input_azimuths: np.array, glows
 
 class RectangularSurvivalProbabilityPointingSet(PointingSet):
     def __init__(self, l1c_dataset: InputRectangularPointingSet, sensor: Sensor, spin_phase: SpinPhase,
-                 glows_dataset: GlowsL3eRectangularMapInputData, energies: np.ndarray, cg_corrected: bool = False):
+                 glows_dataset: GlowsL3eRectangularMapInputData | None, energies: np.ndarray, cg_corrected: bool = False):
         num_spin_angle_bins = l1c_dataset.exposure_times.shape[-1]
 
         deg_spacing = 360 / num_spin_angle_bins
@@ -135,7 +137,7 @@ def __init__(self, l1c_dataset: InputRectangularPointingSet, sensor: Sensor, spi
                     self.azimuths, glows_dataset.spin_angle, sp_interpolated_to_hi_energies)
                 sp_final = sp_interpolated_to_pset_angles
         else:
-            sp_final = np.full((1, len(energies), 3600), np.nan)
+            sp_final = np.ones((1, len(energies), 3600))
 
         dataset["survival_probability_times_exposure"] = xr.DataArray(
             sp_final * exposure,

imap_l3_processing/maps/sp_map_initializer.py

@@ -34,11 +34,10 @@ def get_maps_that_can_be_produced(self, l3_descriptor: str) -> list[PossibleMapT
         map_duration = get_duration_from_map_descriptor(l3_descriptor_parts)
 
         glows_file_by_repointing = self._collect_glows_psets_by_repoint(l3_descriptor_parts)
+
         if len(glows_file_by_repointing) == 0:
-            logger.info(f"No GLOWS data available for descriptor {l3_descriptor}, no maps will be produced!")
-            return []
-        glows_start_dates = [ScienceFilePath(glows_path).start_date for glows_path in glows_file_by_repointing.values()]
-        glows_data_end_date = datetime.strptime(max(glows_start_dates), "%Y%m%d")
+            logger.info(f"No GLOWS data available for descriptor {l3_descriptor}. "
+                        f"Processing will use default survival probability of 1.0.")
 
         l2_descriptors = self._get_l2_dependencies(l3_descriptor_parts)
         l2_descriptor_strs = [map_descriptor_parts_to_string(parts) for parts in l2_descriptors]
@@ -51,9 +50,6 @@ def get_maps_that_can_be_produced(self, l3_descriptor: str) -> list[PossibleMapT
         possible_maps = []
         for str_start_date in sorted(list(possible_start_dates)):
             start_date = datetime.strptime(str_start_date, "%Y%m%d")
-            if start_date + map_duration > glows_data_end_date:
-                logger.info(f"Not enough GLOWS data to produce map {l3_descriptor} {str_start_date}")
-                continue
 
             l2_file_paths = []
             for l2_descriptor in l2_descriptor_strs:
@@ -83,14 +79,13 @@ def get_maps_that_can_be_produced(self, l3_descriptor: str) -> list[PossibleMapT
             glows_files = [glows_file_by_repointing[repoint] for repoint in l1c_repointings if
                            repoint in glows_file_by_repointing]
 
-            if len(glows_files) > 0:
-                input_metadata = InputMetadata(instrument=self.instrument, data_level='l3', start_date=start_date,
-                                               end_date=start_date + map_duration, version='v001',
-                                               descriptor=l3_descriptor)
+            input_metadata = InputMetadata(instrument=self.instrument, data_level='l3', start_date=start_date,
+                                           end_date=start_date + map_duration, version='v001',
+                                           descriptor=l3_descriptor)
 
-                possible_map_to_produce = PossibleMapToProduce(
-                    input_files=set(l2_file_paths + glows_files + l1c_names + self._get_ancillary_files()),
-                    input_metadata=input_metadata
-                )
-                possible_maps.append(possible_map_to_produce)
+            possible_map_to_produce = PossibleMapToProduce(
+                input_files=set(l2_file_paths + glows_files + l1c_names + self._get_ancillary_files()),
+                input_metadata=input_metadata
+            )
+            possible_maps.append(possible_map_to_produce)
         return possible_maps

imap_l3_processing/ultra/science/ultra_survival_probability.py

@@ -15,35 +15,41 @@
 
 
 class UltraSurvivalProbability(UltraPointingSet):
-    def __init__(self, l1c_pset: UltraL1CPSet, l3e_glows: UltraGlowsL3eData, bin_groups: Optional[np.ndarray] = None):
+    def __init__(self, l1c_pset: UltraL1CPSet, l3e_glows: Optional[UltraGlowsL3eData],
+                 bin_groups: Optional[np.ndarray] = None):
         coarse_bins = bin_pset_energy_bins(l1c_pset.to_xarray(), bin_groups)
         super().__init__(coarse_bins, geometry.SpiceFrame.IMAP_DPS)
 
-        l1c_epoch_in_et = spiceypy.unitim(self.data.coords[CoordNames.TIME.value].values[0] / ONE_SECOND_IN_NANOSECONDS,
-                                          "TT", "ET")
-        rotated_az_el_points = geometry.frame_transform_az_el(l1c_epoch_in_et, self.az_el_points,
-                                                              geometry.SpiceFrame.IMAP_DPS,
-                                                              geometry.SpiceFrame.ECLIPJ2000)
-        glows_nside = npix_to_nside(len(l3e_glows.healpix_index))
-        glows_healpix = HEALPix(nside=glows_nside)
+        if l3e_glows is None:
+            num_pixels = len(l1c_pset.healpix_index)
+            num_energies = len(coarse_bins.energy_bin_geometric_mean)
+            energy_interpolated_sp = np.ones((num_energies, num_pixels))
+        else:
+            l1c_epoch_in_et = spiceypy.unitim(
+                self.data.coords[CoordNames.TIME.value].values[0] / ONE_SECOND_IN_NANOSECONDS, "TT", "ET")
+            rotated_az_el_points = geometry.frame_transform_az_el(l1c_epoch_in_et, self.az_el_points,
+                                                                  geometry.SpiceFrame.IMAP_DPS,
+                                                                  geometry.SpiceFrame.ECLIPJ2000)
+            glows_nside = npix_to_nside(len(l3e_glows.healpix_index))
+            glows_healpix = HEALPix(nside=glows_nside)
 
-        npixels = len(l1c_pset.healpix_index)
+            npixels = len(l1c_pset.healpix_index)
 
-        spatially_interpolated_sp = np.zeros((len(l3e_glows.energy), npixels))
+            spatially_interpolated_sp = np.zeros((len(l3e_glows.energy), npixels))
 
-        for energy_index in range(len(l3e_glows.energy)):
-            spatially_interpolated_sp[energy_index, :] = glows_healpix.interpolate_bilinear_lonlat(
-                Quantity(rotated_az_el_points[:, 0], unit='deg'),
-                Quantity(rotated_az_el_points[:, 1], unit='deg'),
-                l3e_glows.survival_probability[0, energy_index, :])
+            for energy_index in range(len(l3e_glows.energy)):
+                spatially_interpolated_sp[energy_index, :] = glows_healpix.interpolate_bilinear_lonlat(
+                    Quantity(rotated_az_el_points[:, 0], unit='deg'),
+                    Quantity(rotated_az_el_points[:, 1], unit='deg'),
+                    l3e_glows.survival_probability[0, energy_index, :])
 
-        energy_interpolated_sp = np.zeros((len(coarse_bins.energy_bin_geometric_mean), npixels))
-        for healpix_index in range(npixels):
-            energy_interpolated_sp[:, healpix_index] = np.interp(
-                np.log10(coarse_bins.energy_bin_geometric_mean),
-                np.log10(l3e_glows.energy),
-                spatially_interpolated_sp[:, healpix_index]
-            )
+            energy_interpolated_sp = np.zeros((len(coarse_bins.energy_bin_geometric_mean), npixels))
+            for healpix_index in range(npixels):
+                energy_interpolated_sp[:, healpix_index] = np.interp(
+                    np.log10(coarse_bins.energy_bin_geometric_mean),
+                    np.log10(l3e_glows.energy),
+                    spatially_interpolated_sp[:, healpix_index]
+                )
 
         self.data["survival_probability_times_exposure"] = (
             [

imap_l3_processing/utils.py

@@ -228,8 +228,8 @@ def combine_glows_l3e_with_l1c_pointing(glows_l3e_data: list[GlowsL3eData], l1c_
     l1c_by_repoint = {l1c.repointing: l1c for l1c in l1c_data}
     glows_by_repoint = {l3e.repointing: l3e for l3e in glows_l3e_data}
 
-    return [(l1c_by_repoint[repoint], glows_by_repoint[repoint]) for repoint in l1c_by_repoint.keys() if
-            repoint in glows_by_repoint]
+    return [(l1c_by_repoint[repoint], glows_by_repoint.get(repoint, None))
+            for repoint in l1c_by_repoint.keys()]
 
 
 def get_dependency_paths_by_descriptor(deps: ProcessingInputCollection, descriptors: list[str]) -> dict[

tests/hi/test_hi_sp_initializer.py

@@ -83,7 +83,40 @@ def test_get_maps_that_can_be_produced(self, mock_read_cdf_parents):
                     version="v001",
                     descriptor='h45-ena-h-sf-sp-anti-hae-4deg-3mo',
                 )
-            )
+            ),
+            PossibleMapToProduce(
+                input_files={
+                    'imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100701_v001.cdf',
+                    'imap_glows_l3e_survival-probability-hi-45_20100702-repoint00202_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20100401-repoint00201_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20100402-repoint00202_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20100403-repoint00203_v001.cdf',
+                },
+                input_metadata=InputMetadata(
+                    instrument="hi",
+                    data_level="l3",
+                    start_date=datetime(2010, 7, 1),
+                    end_date=datetime(2010, 9, 30, 7, 30),
+                    version="v001",
+                    descriptor='h45-ena-h-sf-sp-anti-hae-4deg-3mo',
+                )
+            ),
+            PossibleMapToProduce(
+                input_files={
+                    'imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20101001_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20101001-repoint00301_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20101002-repoint00302_v001.cdf',
+                    'imap_hi_l1c_45sensor-pset_20101003-repoint00303_v001.cdf',
+                },
+                input_metadata=InputMetadata(
+                    instrument="hi",
+                    data_level="l3",
+                    start_date=datetime(2010, 10, 1),
+                    end_date=datetime(2010, 12, 31, 7, 30),
+                    version="v001",
+                    descriptor='h45-ena-h-sf-sp-anti-hae-4deg-3mo',
+                )
+            ),
         ]
 
         initializer = HiSPInitializer()
@@ -100,10 +133,39 @@ def test_get_maps_that_can_be_produced(self, mock_read_cdf_parents):
         mock_read_cdf_parents.assert_has_calls([
             call('imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100101_v001.cdf'),
             call('imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100401_v001.cdf'),
+            call('imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100701_v001.cdf'),
+            call('imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20101001_v001.cdf'),
         ], any_order=True)
 
         self.assertEqual(expected_possible_maps, actual_possible_maps)
 
+    @patch('imap_l3_processing.maps.map_initializer.read_cdf_parents')
+    def test_get_maps_that_can_be_produced_with_no_glows_data(self, mock_read_cdf_parents):
+        self.mock_query.side_effect = [
+            create_mock_query_results([]),
+            create_mock_query_results([]),
+            create_mock_query_results([
+                'imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100101_v001.cdf',
+                'imap_hi_l2_h45-ena-h-sf-nsp-anti-hae-4deg-3mo_20100401_v001.cdf',
+            ]),
+            create_mock_query_results([])
+        ]
+
+        mock_read_cdf_parents.side_effect = self.create_fake_read_cdf_parents("45")
+
+        initializer = HiSPInitializer()
+
+        actual_possible_maps = initializer.get_maps_that_can_be_produced('h45-ena-h-sf-sp-anti-hae-4deg-3mo')
+
+        self.assertEqual(2, len(actual_possible_maps))
+        for possible_map in actual_possible_maps:
+            glows_files = [f for f in possible_map.input_files if 'glows' in f]
+            self.assertEqual([], glows_files)
+            l2_files = [f for f in possible_map.input_files if '_l2_' in f]
+            self.assertGreater(len(l2_files), 0)
+            l1c_files = [f for f in possible_map.input_files if '_l1c_' in f]
+            self.assertGreater(len(l1c_files), 0)
+
     @patch('imap_l3_processing.maps.map_initializer.read_cdf_parents')
     def test_get_maps_that_can_be_produced_full_spin_descriptor(self, mock_read_cdf_parents):
         self.mock_query.side_effect = [
@@ -275,7 +337,40 @@ def test_get_maps_that_should_be_produced(self, mock_read_cdf_parents):
                     version="v002",
                     descriptor='h90-ena-h-sf-sp-anti-hae-4deg-3mo',
                 )
-            )
+            ),
+            PossibleMapToProduce(
+                input_files={
+                    'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20100701_v001.cdf',
+                    'imap_glows_l3e_survival-probability-hi-90_20100702-repoint00201_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20100401-repoint00201_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20100402-repoint00202_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20100403-repoint00203_v001.cdf',
+                },
+                input_metadata=InputMetadata(
+                    instrument="hi",
+                    data_level="l3",
+                    start_date=datetime(2010, 7, 1),
+                    end_date=datetime(2010, 9, 30, 7, 30),
+                    version="v001",
+                    descriptor='h90-ena-h-sf-sp-anti-hae-4deg-3mo',
+                )
+            ),
+            PossibleMapToProduce(
+                input_files={
+                    'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20101001_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20101001-repoint00301_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20101002-repoint00302_v001.cdf',
+                    'imap_hi_l1c_90sensor-pset_20101003-repoint00303_v001.cdf',
+                },
+                input_metadata=InputMetadata(
+                    instrument="hi",
+                    data_level="l3",
+                    start_date=datetime(2010, 10, 1),
+                    end_date=datetime(2010, 12, 31, 7, 30),
+                    version="v001",
+                    descriptor='h90-ena-h-sf-sp-anti-hae-4deg-3mo',
+                )
+            ),
         ]
 
         initializer = HiSPInitializer()
@@ -292,6 +387,8 @@ def test_get_maps_that_should_be_produced(self, mock_read_cdf_parents):
         mock_read_cdf_parents.assert_has_calls([
             call(f'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20100101_v001.cdf'),
             call(f'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20100401_v001.cdf'),
+            call(f'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20100701_v001.cdf'),
+            call(f'imap_hi_l2_h90-ena-h-sf-nsp-anti-hae-4deg-3mo_20101001_v001.cdf'),
 
             call(f'imap_hi_l3_h90-ena-h-sf-sp-anti-hae-4deg-3mo_20100101_v001.cdf'),
             call(f'imap_hi_l3_h90-ena-h-sf-sp-anti-hae-4deg-3mo_20100401_v001.cdf'),

tests/lo/l3/test_lo_sp_initializer.py

@@ -55,6 +55,16 @@ def test_get_maps_that_should_be_produced(self,  mock_read_cdf_parents):
                 'imap_lo_l1c_pset_20100402-repoint00102_v001.cdf',
                 'imap_lo_l1c_pset_20100403-repoint00103_v001.cdf'
             ],
+            [
+                'imap_lo_l1c_pset_20100701-repoint00201_v001.cdf',
+                'imap_lo_l1c_pset_20100702-repoint00202_v001.cdf',
+                'imap_lo_l1c_pset_20100703-repoint00203_v001.cdf'
+            ],
+            [
+                'imap_lo_l1c_pset_20101001-repoint00301_v001.cdf',
+                'imap_lo_l1c_pset_20101002-repoint00302_v001.cdf',
+                'imap_lo_l1c_pset_20101003-repoint00303_v001.cdf'
+            ],
             [
                 'imap_glows_l3e_survival-probability-lo_20100101-repoint00001_v001.cdf',
                 'imap_glows_l3e_survival-probability-lo_20100102-repoint00002_v001.cdf',
@@ -87,6 +97,8 @@ def test_get_maps_that_should_be_produced(self,  mock_read_cdf_parents):
         mock_read_cdf_parents.assert_has_calls([
             call('imap_lo_l2_l090-ena-h-sf-nsp-ram-hae-6deg-12mo_20100101_v001.cdf'),
             call('imap_lo_l2_l090-ena-h-sf-nsp-ram-hae-6deg-12mo_20100401_v001.cdf'),
+            call('imap_lo_l2_l090-ena-h-sf-nsp-ram-hae-6deg-12mo_20100701_v001.cdf'),
+            call('imap_lo_l2_l090-ena-h-sf-nsp-ram-hae-6deg-12mo_20101001_v001.cdf'),
 
             call('imap_lo_l3_l090-ena-h-sf-sp-ram-hae-6deg-12mo_20100101_v001.cdf'),
             call('imap_lo_l3_l090-ena-h-sf-sp-ram-hae-6deg-12mo_20100401_v001.cdf')
@@ -112,7 +124,8 @@ def test_get_maps_that_should_be_produced(self,  mock_read_cdf_parents):
             }
         )
 
-        self.assertEqual([expected_possible_map_to_produce], actual_maps_to_produce)
+        self.assertIn(expected_possible_map_to_produce, actual_maps_to_produce)
+        self.assertEqual(3, len(actual_maps_to_produce))
 
     @patch('imap_l3_processing.lo.l3.lo_sp_initializer.furnish_spice_metakernel')
     def test_furnish_spice_dependencies(self, mock_furnish_metakernel):

tests/maps/test_rectangular_survival_probability.py

@@ -365,6 +365,15 @@ def test_exposure_weighting_with_interpolated_survival_probabilities(self):
                     expected_mask
                 )
 
+    def test_uses_default_survival_probability_of_one_when_glows_is_none(self):
+        pointing_set = RectangularSurvivalProbabilityPointingSet(
+            self.l1c_hi_dataset, Sensor.Hi90, SpinPhase.RamOnly,
+            glows_dataset=None, energies=self.hi_energies)
+
+        sp_times_exposure = pointing_set.data["survival_probability_times_exposure"].values
+        expected = np.ones((1, self.num_energies, 3600)) * self.l1c_hi_dataset.exposure_times
+        np.testing.assert_array_equal(sp_times_exposure, expected)
+
     def test_exposure_weighted_survivals_are_repeated_to_match_l1c_shape(self):
         pointing_set = RectangularSurvivalProbabilityPointingSet(self.l1c_hi_dataset, Sensor.Hi90, SpinPhase.RamOnly,
                                                                  self.glows_data,

tests/maps/test_survival_probability_processing.py

@@ -145,12 +145,13 @@ def test_integration_uses_fill_values_for_missing_l3e_data(self):
 
         descriptor = parse_map_descriptor("h90-ena-h-sf-nsp-ram-hae-4deg-3mo")
         survival_dependencies = HiLoL3SurvivalDependencies(l2_intensity_map, l1c_psets, l3e_psets, descriptor,
-                                                           l2_intensity_map.intensity_map_data.energy)
+                                                           Mock())
 
         expected_sky_strip = np.full(45, 2.0)
+        default_sp_sky_strip = np.full(45, 1.0)
 
         output_map = process_survival_probabilities(survival_dependencies, SpiceFrame.IMAP_HAE)
         np.testing.assert_equal(output_map.intensity_map_data.ena_intensity[0, 0, 76, :], expected_sky_strip)
         np.testing.assert_equal(output_map.intensity_map_data.ena_intensity[0, 0, 78, :], expected_sky_strip)
-        np.testing.assert_equal(output_map.intensity_map_data.ena_intensity[0, 0, 79, :], np.full(45, np.nan))
+        np.testing.assert_equal(output_map.intensity_map_data.ena_intensity[0, 0, 79, :], default_sp_sky_strip)
         np.testing.assert_equal(output_map.intensity_map_data.ena_intensity[0, 0, 80, :], expected_sky_strip)

tests/test_utils.py

@@ -479,7 +479,9 @@ def test_combine_glows_l3e_with_l1c_pointing(self):
 
         expected = [
             (hi_l1c_data[0], glows_l3e_data[0]),
+            (hi_l1c_data[1], None),
             (hi_l1c_data[2], glows_l3e_data[3]),
+            (hi_l1c_data[3], None),
         ]
 
         actual = combine_glows_l3e_with_l1c_pointing(glows_l3e_data, hi_l1c_data)