Center pointings

litebird_sim/observations.py

@@ -6,6 +6,7 @@
 import astropy.time
 import numpy as np
 import numpy.typing as npt
+from ducc0.healpix import Healpix_Base
 
 from .coordinates import DEFAULT_TIME_SCALE
 from .detectors import DetectorInfo, InstrumentInfo
@@ -14,7 +15,7 @@
 from .input_sky import SkyGenerationParams
 from .maps_and_harmonics import HealpixMap, SphericalHarmonics
 from .mpi import MPI_COMM_GRID, _SerialMpiCommunicator
-from .pointings import PointingProvider, DEFAULT_INTERNAL_BUFFER_SIZE_FOR_POINTINGS_MB
+from .pointings import DEFAULT_INTERNAL_BUFFER_SIZE_FOR_POINTINGS_MB, PointingProvider
 from .scanning import RotQuaternion
 from .units import Units
 
@@ -1232,6 +1233,40 @@ def precompute_pointings(
         self.pointing_matrix = pointing_matrix
         self.hwp_angle = hwp_angle
 
+    def center_pointings(
+        self,
+        nside,
+        pointings_dtype=np.float64,
+    ) -> None:
+        """Force the pointings to the center of pixels for a given nside.
+
+        This method can be useful to ensure that the pixels used when scanning
+        are exactly the same as the ones used at the map-making step.
+
+        Args:
+            nside (int):
+                The nside of the map to whose pixels it centers the pointings.
+            pointings_dtype (data-type, optional):
+                Data type to use for the computed arrays. Defaults to `np.float64`.
+
+        Raises:
+            AssertionError:
+                If `precompute_pointings()` has not been called prior to this method,
+                i.e., if `self.pointing_matrix is not defined.
+
+        Notes:
+            This method must be called after `precompute_pointings()`.
+        """
+
+        assert "pointing_matrix" in dir(self), (
+            "you must call precompute_pointings() on a set of observations "
+            "before calling center_pointings()"
+        )
+
+        hpx = Healpix_Base(nside, "RING")
+        pixel_ind = hpx.ang2pix(self.pointing_matrix[:, :, 0:2])
+        self.pointing_matrix[:, :, 0:2] = hpx.pix2ang(pixel_ind)
+
     def _set_mpi_subcommunicators(self):
         """This function splits the global MPI communicator into three kinds of
         sub-communicators:

litebird_sim/pointings_in_obs.py

@@ -1,20 +1,17 @@
 from collections.abc import Callable
 
+import astropy.time
 import numpy as np
 import numpy.typing as npt
-import astropy.time
-
 from deprecated import deprecated
-
 from ducc0.healpix import Healpix_Base
 
+from .coordinates import CoordinateSystem, rotate_coordinates_e2g
 from .detectors import InstrumentInfo
 from .hwp import HWP
 from .observations import Observation
 from .scanning import RotQuaternion
 
-from .coordinates import CoordinateSystem, rotate_coordinates_e2g
-
 
 def prepare_pointings(
     observations: Observation | list[Observation],
@@ -112,6 +109,53 @@ def precompute_pointings(
         cur_obs.precompute_pointings(pointings_dtype=pointings_dtype)
 
 
+def center_pointings(
+    observations: Observation | list[Observation],
+    nside: int,
+    pointings_dtype=np.float64,
+) -> None:
+    """Force the pointings to the center of pixels for a given nside.
+
+    This function triggers the change of the full time-domain pointing matrix
+    and, The results are stored internally in each observation's
+    ``pointing_matrix`` field.
+
+    This method can be useful to ensure that the pixels used when scanning
+    are exactly the same as the ones used at the map-making step.
+
+    Args:
+        observations (Observation or list[Observation]):
+            A single observation or a list of observations for which pointings should be precomputed.
+
+        nside (int):
+            The nside of the map to whose pixels it centers the pointings.
+
+        pointings_dtype (data-type, optional):
+            Data type to use when allocating the pointing and HWP arrays.
+            Defaults to `np.float64`.
+
+    Returns:
+        None
+
+    Raises:
+        AssertionError:
+            If any observation does not have a pointing matrix defined.
+            Make sure to call :func:`precompute_pointings()` beforehand.
+
+    Notes:
+        - This function must be called after pointing precomputation (i.e., after `precompute_pointings()`).
+        - Output arrays are stored in `Observation.pointing_matrix`.
+    """
+
+    if isinstance(observations, Observation):
+        obs_list = [observations]
+    else:
+        obs_list = observations
+
+    for cur_obs in obs_list:
+        cur_obs.center_pointings(nside=nside, pointings_dtype=pointings_dtype)
+
+
 @deprecated(
     version="0.15.0",
     reason="This function adds the HWP angle to the orientation, but this is logically wrong "

litebird_sim/simulations.py

@@ -59,7 +59,7 @@
 from .noise import add_noise_to_observations
 from .non_linearity import NonLinParams, apply_quadratic_nonlin_to_observations
 from .observations import Observation, TodDescription
-from .pointings_in_obs import precompute_pointings, prepare_pointings
+from .pointings_in_obs import precompute_pointings, prepare_pointings, center_pointings
 from .profiler import TimeProfiler, profile_list_to_speedscope
 from .scan_map import scan_map_in_observations
 from .scanning import ScanningStrategy, SpinningScanningStrategy
@@ -1667,6 +1667,14 @@ def precompute_pointings(self, pointings_dtype=np.float64) -> None:
             observations=self.observations, pointings_dtype=pointings_dtype
         )
 
+    def center_pointings(self, nside, pointings_dtype=np.float64) -> None:
+        """Force the pointings to the center of pixels for a given nside.
+        Changes the values in the field ``pointing matrix
+        """
+        center_pointings(
+            observations=self.observations, nside=nside, pointings_dtype=pointings_dtype
+        )
+
     @_profile
     def compute_pos_and_vel(
         self,

test/test_simulations.py

@@ -3,15 +3,15 @@
 import os
 import pathlib
 from pathlib import Path
-from tempfile import TemporaryDirectory, NamedTemporaryFile
+from tempfile import NamedTemporaryFile, TemporaryDirectory
+from unittest.mock import MagicMock, patch
 from uuid import UUID
 
 import astropy
+import litebird_sim as lbs
 import numpy as np
 import pytest
-from unittest.mock import MagicMock, patch
-
-import litebird_sim as lbs
+from ducc0.healpix import Healpix_Base
 from litebird_sim.detectors import DetectorInfo, FreqChannelInfo
 
 
@@ -21,7 +21,7 @@ def savefig(*args, **kwargs):
 
 
 def test_healpix_map_write(tmp_path):
-    sim = lbs.Simulation(base_path=tmp_path / "simulation_dir", random_seed=12345)
+    sim = lbs.Simulation(base_path="./simulation_dir", random_seed=12345)
     output_file = sim.write_healpix_map(filename="test.fits.gz", pixels=np.zeros(12))
 
     assert isinstance(output_file, pathlib.Path)
@@ -40,7 +40,7 @@ def test_healpix_map_write(tmp_path):
 
 def test_markdown_report(tmp_path):
     sim = lbs.Simulation(
-        base_path=tmp_path / "simulation_dir",
+        base_path="./simulation_dir",
         name="My simulation",
         description="Lorem ipsum",
         start_time=1.0,
@@ -98,7 +98,7 @@ def test_imo_in_report(tmp_path):
     imo = lbs.Imo(flatfile_location=curpath / "test_imo")
 
     sim = lbs.Simulation(
-        base_path=tmp_path / "simulation_dir",
+        base_path="./simulation_dir",
         name="My simulation",
         description="Lorem ipsum",
         imo=imo,
@@ -428,7 +428,7 @@ def _configure_simulation_for_pointings(
     )
 
     sim = lbs.Simulation(
-        base_path=tmp_path / "simulation_dir",
+        base_path="./simulation_dir",
         start_time=0.0,
         duration_s=61.0,
         random_seed=12345,
@@ -887,3 +887,38 @@ def det_side_effect(imo, url):
         # Asking for 5 detectors when only 2 exist in the channel
         with pytest.raises(ValueError, match="Expected 5 detectors, but got 2"):
             sim.set_detectors(channels=["url_ch_a"], detectors=5)
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_center_pointings(tmp_path, dtype):
+    sim = _configure_simulation_for_pointings(
+        tmp_path,
+        include_hwp=False,
+        store_full_pointings=True,
+        dtype=dtype,
+    )
+
+    sim.center_pointings(nside=16, pointings_dtype=dtype)
+
+    for cur_obs in sim.observations:
+        assert "pointing_matrix" in dir(cur_obs)
+        assert cur_obs.pointing_matrix.dtype == dtype
+        assert cur_obs.pointing_matrix.shape == (
+            cur_obs.n_detectors,
+            cur_obs.n_samples,
+            3,
+        )
+
+        # confirming that the pointings are centered
+        # by confirming they are exactly the same after
+        # doing ang2pix and pix2ang in sequence
+        hpx = Healpix_Base(16, "RING")
+        aux_pointings = cur_obs.pointing_matrix.copy()
+        aux_pointings[:, :, 0:2] = hpx.pix2ang(
+            hpx.ang2pix(cur_obs.pointing_matrix[:, :, 0:2])
+        )
+
+        np.testing.assert_allclose(
+            cur_obs.pointing_matrix,
+            aux_pointings,
+        )