Skip to content

Update Ultra maps to rectangular projection (for all but intensity) #96

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pleasant-menlo wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Update Ultra maps to rectangular projection (for all but intensity) #96

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
876a3b0
Harrison 3158 - LFAR - using rectangular L2 ultra maps for everything…
pleasant-menlo May 6, 2026
d3aceef
Harrison 3158 - LFAR - ultra processing happening using rectangular i…
pleasant-menlo May 7, 2026
7bc6106
Merge branch 'main' into 3158-rectangular
pleasant-menlo May 7, 2026
39ab7f1
Harrison 3158 - LFAR - ultra processing, using converted rectangular …
pleasant-menlo May 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 30 additions & 14 deletions imap_l3_processing/ultra/science/ultra_combined.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,32 @@
from imap_l3_processing.maps.map_models import HealPixIntensityMapData
from imap_l3_processing.maps.map_models import HealPixIntensityMapData, RectangularIntensityMapData, IntensityMapData, \
HealPixCoords

def map_combined_rectangular_quantities_to_healpix_intensity_map(healpix_map: HealPixIntensityMapData,
rectangular_map: RectangularIntensityMapData) -> HealPixIntensityMapData:
healpix_map_intensities = healpix_map.intensity_map_data
rectangular_intensity_data = rectangular_map.intensity_map_data

def combine_sensors(u45: HealPixIntensityMapData, u90: HealPixIntensityMapData) -> HealPixIntensityMapData:
intensity = (((u45.intensity_map_data.ena_intensity * u45.intensity_map_data.exposure_factor) +
(u90.intensity_map_data.ena_intensity * u90.intensity_map_data.exposure_factor)) /
(u45.intensity_map_data.exposure_factor + u90.intensity_map_data.exposure_factor))

intensity_systematic_error = (((u45.intensity_map_data.ena_intensity_sys_err * u45.intensity_map_data.exposure_factor) +
(u90.intensity_map_data.ena_intensity_sys_err * u90.intensity_map_data.exposure_factor)) /
(u45.intensity_map_data.exposure_factor + u90.intensity_map_data.exposure_factor))

intensity_static_error = (
((u45.intensity_map_data.ena_intensity_stat_uncert * u45.intensity_map_data.exposure_factor) +
(u90.intensity_map_data.ena_intensity_stat_uncert * u90.intensity_map_data.exposure_factor)) /
(u45.intensity_map_data.exposure_factor + u90.intensity_map_data.exposure_factor))
healpix_map_data = HealPixIntensityMapData(
intensity_map_data=IntensityMapData(
ena_intensity_stat_uncert=healpix_map_intensities.ena_intensity_stat_uncert,
ena_intensity_sys_err=healpix_map_intensities.ena_intensity_sys_err,
ena_intensity=healpix_map_intensities.ena_intensity,
epoch=rectangular_intensity_data.epoch,
epoch_delta=rectangular_intensity_data.epoch_delta,
energy=rectangular_intensity_data.energy,
energy_delta_plus=rectangular_intensity_data.energy_delta_plus,
energy_delta_minus=rectangular_intensity_data.energy_delta_minus,
energy_label=rectangular_intensity_data.energy_label,
latitude=rectangular_intensity_data.latitude,
longitude=rectangular_intensity_data.longitude,
exposure_factor=rectangular_intensity_data.exposure_factor,
obs_date=rectangular_intensity_data.obs_date,
obs_date_range=rectangular_intensity_data.obs_date_range,
solid_angle=rectangular_intensity_data.solid_angle,
),
coords=HealPixCoords(
pixel_index=healpix_map.coords.pixel_index,
pixel_index_label=healpix_map.coords.pixel_index_label,
)
)
return healpix_map_data
18 changes: 13 additions & 5 deletions imap_l3_processing/ultra/ultra_l3_dependencies.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,8 @@

@dataclass
class UltraL3Dependencies:
ultra_l2_map: HealPixIntensityMapData
ultra_l2_healpix_map: HealPixIntensityMapData
ultra_l2_rectangular_map: RectangularIntensityMapData
ultra_l1c_pset: list[UltraL1CPSet]
glows_l3e_sp: list[UltraGlowsL3eData]
energy_bin_group_sizes: Optional[np.ndarray]
Expand Down Expand Up @@ -58,12 +59,13 @@ def from_file_paths(cls, l2_map_path: Path, l1c_file_paths: list[Path], glows_fi
l1c_paths_dict = {path.stem: path for path in l1c_file_paths}
l2_healpix_datasets = ultra_l2(l1c_paths_dict)
l2_healpix_map_data = HealPixIntensityMapData.read_from_xarray(l2_healpix_datasets[0])
l2_rectangular_map_data = RectangularIntensityMapData.read_from_path(l2_map_path)
energy_bin_group_sizes = None
if energy_bin_path:
paths.append(energy_bin_path)
energy_bin_group_sizes = np.loadtxt(energy_bin_path, delimiter=',', dtype=np.uint8)

return cls(ultra_l2_map=l2_healpix_map_data, ultra_l1c_pset=ultra_l1c_data, glows_l3e_sp=glows_l3e_data,
return cls(ultra_l2_healpix_map=l2_healpix_map_data, ultra_l2_rectangular_map=l2_rectangular_map_data, ultra_l1c_pset=ultra_l1c_data, glows_l3e_sp=glows_l3e_data,
dependency_file_paths=paths, energy_bin_group_sizes=energy_bin_group_sizes)


Expand Down Expand Up @@ -98,8 +100,10 @@ def get_fit_energy_ranges(self) -> np.ndarray:

@dataclass
class UltraL3CombinedDependencies:
u45_l2_map: HealPixIntensityMapData
u90_l2_map: HealPixIntensityMapData
u45_l2_healpix_map: HealPixIntensityMapData
u90_l2_healpix_map: HealPixIntensityMapData
u45_l2_rectangular_map: RectangularIntensityMapData
u90_l2_rectangular_map: RectangularIntensityMapData
u45_l1c_psets: list[UltraL1CPSet]
u90_l1c_psets: list[UltraL1CPSet]
glows_l3e_psets: list[UltraGlowsL3eData]
Expand Down Expand Up @@ -159,6 +163,9 @@ def from_file_paths(cls, u45_pset_paths: list[Path], u90_pset_paths: list[Path],
l2_u90_maps = ultra_l2(l1c_u90_paths_dict)
l2_u90_healpix_map_data = HealPixIntensityMapData.read_from_xarray(l2_u90_maps[0])

l2_u45_rectangular_map_data = RectangularIntensityMapData.read_from_path(u45_map_path)
l2_u90_rectangular_map_data = RectangularIntensityMapData.read_from_path(u90_map_path)

dependency_paths = [*u45_pset_paths, *u90_pset_paths, *glows_l3e_pset_paths, u45_map_path,
u90_map_path]
if energy_bin_group_sizes_path:
Expand All @@ -167,7 +174,8 @@ def from_file_paths(cls, u45_pset_paths: list[Path], u90_pset_paths: list[Path],
else:
energy_bin_group_sizes = None

return cls(u45_l2_map=l2_u45_healpix_map_data, u90_l2_map=l2_u90_healpix_map_data,
return cls(u45_l2_healpix_map=l2_u45_healpix_map_data, u90_l2_healpix_map=l2_u90_healpix_map_data,
u45_l2_rectangular_map=l2_u45_rectangular_map_data, u90_l2_rectangular_map=l2_u90_rectangular_map_data,
u45_l1c_psets=u45_l1c_psets, u90_l1c_psets=u90_l1c_psets,
glows_l3e_psets=survival_probability_ul_pset, energy_bin_group_sizes=energy_bin_group_sizes,
dependency_file_paths=dependency_paths)
76 changes: 39 additions & 37 deletions imap_l3_processing/ultra/ultra_processor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,8 @@
from datetime import datetime, timedelta
from typing import Optional

import numpy as np
from imap_processing.spice.geometry import SpiceFrame

from imap_l3_processing.constants import TT2000_EPOCH
from imap_l3_processing.maps.map_combination import ExposureWeightedCombination, UncertaintyWeightedCombination
from imap_l3_processing.maps.map_descriptors import MapDescriptorParts, MapQuantity, SurvivalCorrection, \
parse_map_descriptor, PixelSize, Sensor
Expand All @@ -16,6 +14,7 @@
from imap_l3_processing.maps.map_processor import MapProcessor
from imap_l3_processing.maps.spectral_fit import calculate_spectral_index_for_multiple_ranges, \
slice_energy_range_by_bin, fit_spectral_index_map
from imap_l3_processing.ultra.science.ultra_combined import map_combined_rectangular_quantities_to_healpix_intensity_map
from imap_l3_processing.ultra.science.ultra_survival_probability import UltraSurvivalProbabilitySkyMap, \
UltraSurvivalProbability
from imap_l3_processing.ultra.ultra_l3_dependencies import UltraL3Dependencies, UltraL3SpectralIndexDependencies, \
Expand Down Expand Up @@ -56,7 +55,8 @@ def process(self, spice_frame_name: SpiceFrame = SpiceFrame.ECLIPJ2000):
combined_deps = UltraL3CombinedDependencies.fetch_dependencies(self.dependencies)
combined_data = self._process_combined_survival_probability(combined_deps, spice_frame_name)

data_product = self._process_healpix_intensity_to_rectangular(combined_data, parsed_descriptor.grid,
data_product = self._process_healpix_intensity_to_rectangular(combined_data,
parsed_descriptor.grid,
spice_frame_name=spice_frame_name)
data_product.add_paths_to_parents(combined_deps.dependency_file_paths)

Expand All @@ -67,9 +67,17 @@ def process(self, spice_frame_name: SpiceFrame = SpiceFrame.ECLIPJ2000):

combination_strategy = ExposureWeightedCombination()

healpix_intensity_map_data = combination_strategy.combine_healpix_intensity_map_data(
[deps.u45_l2_map, deps.u90_l2_map])
data_product = self._process_healpix_intensity_to_rectangular(healpix_intensity_map_data,
combined_healpix_intensity_map_data = combination_strategy.combine_healpix_intensity_map_data(
[deps.u45_l2_healpix_map, deps.u90_l2_healpix_map])

combined_rectangular_intensity_map_data = combination_strategy.combine_rectangular_intensity_map_data(
[deps.u45_l2_rectangular_map, deps.u90_l2_rectangular_map]
)

healpix_intensity_map_rectangular_quantities = map_combined_rectangular_quantities_to_healpix_intensity_map(
combined_healpix_intensity_map_data, combined_rectangular_intensity_map_data)

data_product = self._process_healpix_intensity_to_rectangular(healpix_intensity_map_rectangular_quantities,
parsed_descriptor.grid,
spice_frame_name=spice_frame_name)
data_product.add_paths_to_parents(deps.dependency_file_paths)
Expand All @@ -81,14 +89,16 @@ def process(self, spice_frame_name: SpiceFrame = SpiceFrame.ECLIPJ2000):

def _process_combined_survival_probability(self, deps: UltraL3CombinedDependencies, spice_frame_name: SpiceFrame):
u45_dep = UltraL3Dependencies(
ultra_l2_map=deps.u45_l2_map,
ultra_l2_healpix_map=deps.u45_l2_healpix_map,
ultra_l2_rectangular_map=deps.u45_l2_rectangular_map,
ultra_l1c_pset=deps.u45_l1c_psets,
glows_l3e_sp=deps.glows_l3e_psets,
dependency_file_paths=deps.dependency_file_paths,
energy_bin_group_sizes=deps.energy_bin_group_sizes,
)
u90_dep = UltraL3Dependencies(
ultra_l2_map=deps.u90_l2_map,
ultra_l2_healpix_map=deps.u90_l2_healpix_map,
ultra_l2_rectangular_map=deps.u90_l2_rectangular_map,
ultra_l1c_pset=deps.u90_l1c_psets,
glows_l3e_sp=deps.glows_l3e_psets,
dependency_file_paths=deps.dependency_file_paths,
Expand All @@ -108,32 +118,33 @@ def _process_survival_probability(self, deps: UltraL3Dependencies,
_l1c, _l3e in
combined_psets]

intensity_data = deps.ultra_l2_map.intensity_map_data
coords = deps.ultra_l2_map.coords
healpix_intensity_data = deps.ultra_l2_healpix_map.intensity_map_data
rectangular_intensity_data = deps.ultra_l2_rectangular_map.intensity_map_data
coords = deps.ultra_l2_healpix_map.coords
corrected_skymap = UltraSurvivalProbabilitySkyMap(survival_probability_psets, spice_frame_name, coords.nside)
survival_probability_map = corrected_skymap.to_dataset()["exposure_weighted_survival_probabilities"].values

corrected_intensity = intensity_data.ena_intensity / survival_probability_map
corrected_stat_uncert = intensity_data.ena_intensity_stat_uncert / survival_probability_map
corrected_sys_unc = intensity_data.ena_intensity_sys_err / survival_probability_map
corrected_intensity = healpix_intensity_data.ena_intensity / survival_probability_map
corrected_stat_uncert = healpix_intensity_data.ena_intensity_stat_uncert / survival_probability_map
corrected_sys_unc = healpix_intensity_data.ena_intensity_sys_err / survival_probability_map

healpix_map_data = HealPixIntensityMapData(
intensity_map_data=IntensityMapData(
ena_intensity_stat_uncert=corrected_stat_uncert,
ena_intensity_sys_err=corrected_sys_unc,
ena_intensity=corrected_intensity,
epoch=intensity_data.epoch,
epoch_delta=intensity_data.epoch_delta,
energy=intensity_data.energy,
energy_delta_plus=intensity_data.energy_delta_plus,
energy_delta_minus=intensity_data.energy_delta_minus,
energy_label=intensity_data.energy_label,
latitude=intensity_data.latitude,
longitude=intensity_data.longitude,
exposure_factor=intensity_data.exposure_factor,
obs_date=intensity_data.obs_date,
obs_date_range=intensity_data.obs_date_range,
solid_angle=intensity_data.solid_angle,
epoch=rectangular_intensity_data.epoch,
epoch_delta=rectangular_intensity_data.epoch_delta,
energy=rectangular_intensity_data.energy,
energy_delta_plus=rectangular_intensity_data.energy_delta_plus,
energy_delta_minus=rectangular_intensity_data.energy_delta_minus,
energy_label=rectangular_intensity_data.energy_label,
latitude=rectangular_intensity_data.latitude,
longitude=rectangular_intensity_data.longitude,
exposure_factor=rectangular_intensity_data.exposure_factor,
obs_date=rectangular_intensity_data.obs_date,
obs_date_range=rectangular_intensity_data.obs_date_range,
solid_angle=rectangular_intensity_data.solid_angle,
),
coords=HealPixCoords(
pixel_index=coords.pixel_index,
Expand Down Expand Up @@ -164,22 +175,14 @@ def _process_healpix_intensity_to_rectangular(self, healpix_map_data: HealPixInt
spacing_deg: int,
spice_frame_name: SpiceFrame) -> RectangularIntensityDataProduct:
variables_to_convert_to_rectangular = [
"exposure_factor",
"ena_intensity",
"ena_intensity_stat_uncert",
"ena_intensity_sys_err",
"obs_date",
"obs_date_range",
]
healpix_map = healpix_map_data.to_healpix_skymap()
rectangular_map, _ = healpix_map.to_rectangular_skymap(spacing_deg, variables_to_convert_to_rectangular)
rectangular_map_xarray_dataset = rectangular_map.to_dataset()

obs_date_seconds = np.ma.masked_invalid(
rectangular_map_xarray_dataset["obs_date"].values / 1e9)
obs_date_seconds.data[np.ma.getmask(obs_date_seconds)] = 0
obs_date = TT2000_EPOCH + timedelta(seconds=1) * obs_date_seconds

input_map_intensity_data = healpix_map_data.intensity_map_data
intensity_map_data = IntensityMapData(
epoch=input_map_intensity_data.epoch,
Expand All @@ -190,10 +193,10 @@ def _process_healpix_intensity_to_rectangular(self, healpix_map_data: HealPixInt
energy_label=input_map_intensity_data.energy_label,
latitude=rectangular_map.sky_grid.el_bin_midpoints,
longitude=rectangular_map.sky_grid.az_bin_midpoints,
obs_date=obs_date,
obs_date_range=rectangular_map_xarray_dataset["obs_date_range"].values,
obs_date=input_map_intensity_data.obs_date,
obs_date_range=input_map_intensity_data.obs_date_range,
solid_angle=rectangular_map.solid_angle_grid.T,
exposure_factor=rectangular_map_xarray_dataset["exposure_factor"].values,
exposure_factor=input_map_intensity_data.exposure_factor,
ena_intensity=rectangular_map_xarray_dataset["ena_intensity"].values,
ena_intensity_stat_uncert=rectangular_map_xarray_dataset["ena_intensity_stat_uncert"].values,
ena_intensity_sys_err=rectangular_map_xarray_dataset["ena_intensity_sys_err"].values,
Expand Down Expand Up @@ -262,7 +265,6 @@ def _process_healpix_spectral_index_to_rectangular(self, healpix_map_data: HealP
)
)


@dataclass
class UltraMapDescriptorParts:
grid_size: int
37 changes: 37 additions & 0 deletions tests/ultra/test_ultra_combined.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
import unittest

import numpy as np

from imap_l3_processing.ultra.science.ultra_combined import map_combined_rectangular_quantities_to_healpix_intensity_map
from tests.maps.test_builders import create_rectangular_intensity_map_data
from tests.ultra.test_ultra_processor import _create_ultra_l2_healpix_data


class TestUltraSurvivalProbability(unittest.TestCase):

def test_map_combined_rectangular_quantities_to_healpix_intensity_map_data(self):
healpix_input = _create_ultra_l2_healpix_data()
rectangular_input = create_rectangular_intensity_map_data()
healpix_input_intensity_map_data = healpix_input.intensity_map_data
rectangular_input_intensity_map_data = rectangular_input.intensity_map_data

result_map = map_combined_rectangular_quantities_to_healpix_intensity_map(healpix_input, rectangular_input)
result_intensity_data = result_map.intensity_map_data

np.testing.assert_array_equal(result_intensity_data.ena_intensity, healpix_input_intensity_map_data.ena_intensity)
np.testing.assert_array_equal(result_intensity_data.ena_intensity_stat_uncert, healpix_input_intensity_map_data.ena_intensity_stat_uncert)
np.testing.assert_array_equal(result_intensity_data.ena_intensity_sys_err, healpix_input_intensity_map_data.ena_intensity_sys_err)

np.testing.assert_array_equal(result_intensity_data.epoch, rectangular_input_intensity_map_data.epoch)
np.testing.assert_array_equal(result_intensity_data.epoch_delta, rectangular_input_intensity_map_data.epoch_delta)
np.testing.assert_array_equal(result_intensity_data.energy, rectangular_input_intensity_map_data.energy)
np.testing.assert_array_equal(result_intensity_data.energy_delta_plus, rectangular_input_intensity_map_data.energy_delta_plus)
np.testing.assert_array_equal(result_intensity_data.energy_delta_minus, rectangular_input_intensity_map_data.energy_delta_minus)
np.testing.assert_array_equal(result_intensity_data.energy_label, rectangular_input_intensity_map_data.energy_label)
np.testing.assert_array_equal(result_intensity_data.latitude, rectangular_input_intensity_map_data.latitude)
np.testing.assert_array_equal(result_intensity_data.longitude, rectangular_input_intensity_map_data.longitude)
np.testing.assert_array_equal(result_intensity_data.obs_date, rectangular_input_intensity_map_data.obs_date)
np.testing.assert_array_equal(result_intensity_data.obs_date_range, rectangular_input_intensity_map_data.obs_date_range)
np.testing.assert_array_equal(result_intensity_data.solid_angle, rectangular_input_intensity_map_data.solid_angle)
np.testing.assert_array_equal(result_intensity_data.exposure_factor, rectangular_input_intensity_map_data.exposure_factor)
self.assertEqual(result_map.coords, healpix_input.coords)
Loading