Merge pull request #2 from Cosmoglobe/interp

Add new argument `interp_kind` to the initialization of `Zodipy`.

docs/index.md

@@ -7,7 +7,7 @@
 [![codecov](https://codecov.io/gh/Cosmoglobe/zodipy/branch/main/graph/badge.svg?token=VZP9L79EUJ)](https://codecov.io/gh/Cosmoglobe/zodipy)
 [![arXiv Paper](https://img.shields.io/badge/arXiv-2205.12962-green)](https://arxiv.org/abs/2205.12962)
 
-ZodiPy simulates zodiacal emission in intensity for arbitrary Solar System observers in the form of timestreams or full-sky maps
+ZodiPy simulates zodiacal emission in intensity for arbitrary solar system observers in the form of timestreams or HEALPix maps.
 ![ZodiPy Logo](img/zodipy_map.png)
 
 
@@ -36,12 +36,12 @@ print(emission)
 
 What's going on here:
 
-- We start by initializing the [`Zodipy`][zodipy.zodipy.Zodipy] class where we specify that we want to use the DIRBE interplanetary dust model.
-- We use the [`get_emission_ang`][zodipy.zodipy.Zodipy.get_emission_ang] method which is a method to compute simulated emission from angular sky coordinates. See the [reference](reference.md) for other available methods.
-- The first argument to the [`get_emission_ang`][zodipy.zodipy.Zodipy.get_emission_ang] method, `25 * u.micron`, specifies the wavelength (or frequency) of the simulated observation. Note that we use Astropy units for many of the input arguments.
-- `theta` and `phi` represent the pointing of the observation (co-latitude and longitude). In this example we observe three sky coordinates.
-- `obs_time` represents the time of observation which is used internally to compute the position of the observer and all other required solar system bodies.
+- We start by initializing the [`Zodipy`][zodipy.zodipy.Zodipy] class, which is our interface, where we specify that we want to use the DIRBE interplanetary dust model.
+- We use the [`get_emission_ang`][zodipy.zodipy.Zodipy.get_emission_ang] method which is a method to simulate emission from angular sky coordinates (see the [reference](reference.md) for other available simulation methods).
+- The first argument to the [`get_emission_ang`][zodipy.zodipy.Zodipy.get_emission_ang] method, `25 * u.micron`, specifies the wavelength of the simulated observation. Note that we use Astropy units for many of the input arguments.
+- `theta` and `phi` represent the pointing of the observation (co-latitude and longitude, following the healpy convention). In this example we observe three sky coordinates.
+- `obs_time` represents the time of observation, which we need to compute the position of the observer and all other required solar system bodies.
 - `obs` represents the observer, and must be an solar system observer supported by the [Astropy ephemeris](https://docs.astropy.org/en/stable/coordinates/solarsystem.html) used internally. If we wish to be more specific about the observer position, we can use the `obs_pos` keyword instead of `obs`, which takes in a heliocentric cartesian position in units of AU.
-- `lonlat` is a boolean which converts the convention of `theta` and `phi` from co-latitude and longitude to longitude and latitude.
+- Finally, `lonlat` is a boolean which converts the convention of `theta` and `phi` from co-latitude and longitude to longitude and latitude.
 
 For more information on using ZodiPy, see [the usage section](usage.md).

docs/introduction.md

@@ -1,17 +1,15 @@
 # Introduction
 
-ZodiPy simulates the zodiacal emission that a solar system observer is predicted to see given an interplanetary dust model. The user selects between a set of built in models, for which the emission can be computed either in the form of timestreams or binned HEALPix maps. 
-
-ZodiPy attempts to make zodiacal emission simulations more accessible by providing the community with a simple Python interface to existing models. For other zodiacal emission tools, see [Zodiacal Light Models on LAMBDA](https://lambda.gsfc.nasa.gov/product/foreground/fg_models.html). ZodiPy is an open source project and all contributions are welcome.
+ZodiPy is an open source Python tool for simulating the zodiacal emission that a solar system observer is predicted to see given an interplanetary dust model. We attempts to make zodiacal emission simulations more accessible by providing the community with a simple interface to existing models. For other zodiacal emission tools, see [Zodiacal Light Models on LAMBDA](https://lambda.gsfc.nasa.gov/product/foreground/fg_models.html). All contributions are most welcome.
 
 
 ## Interplanetary Dust Models
-Currently, ZodiPy supports the following interplanetary dust models:
+ZodiPy supports the following interplanetary dust models:
 
 **1.25-240 $\boldsymbol{\mu}$m**
 
 - DIRBE ([Kelsall et al. 1998](https://ui.adsabs.harvard.edu/abs/1998ApJ...508...44K/abstract))
-- RRM (experimental version in development) ([Rowan-Robinson and May 2013](https://ui.adsabs.harvard.edu/abs/2013MNRAS.429.2894R/abstract))
+- RRM (experimental) ([Rowan-Robinson and May 2013](https://ui.adsabs.harvard.edu/abs/2013MNRAS.429.2894R/abstract))
 
 **100-857 GHz**
 
@@ -22,10 +20,12 @@ Currently, ZodiPy supports the following interplanetary dust models:
 
 !!! info
     The Planck and Odegard models extend the DIRBE interplanetary dust model to CMB frequencies by fitting the blackbody emissivity of the dust in the respective DIRBE interplanetary dust components to Planck HFI data.
+    The distribution of the interplanetary dust is exactly the same as in the DIRBE model.
 
+If you see a missing model, please feel free to contact us by opening an issue on GitHub. 
 
 
 ## Scientific Paper
-For an overview of the ZodiPy model approach and other information regarding zodiacal emission and interplanetary dust modeling we refer to the scientific paper on ZodiPy:
+For an overview of the modeling approach used in ZodiPy and other information regarding zodiacal emission and interplanetary dust modeling we refer to the scientific paper on ZodiPy:
 
 - [Cosmoglobe: Simulating zodiacal emission with ZodiPy](https://arxiv.org/abs/2205.12962)

pyproject.toml

@@ -1,7 +1,7 @@
 [tool.poetry]
 name = "zodipy"
 homepage = "https://github.com/Cosmoglobe/zodipy"
-version = "0.8.3"
+version = "0.8.4"
 description = "Software for simulating zodiacal emission"
 authors = ["Metin San <[email protected]>"]
 readme = "README.md"

tests/test_get_emission.py

@@ -332,7 +332,14 @@ def test_multiprocessing() -> None:
     assert np.allclose(emission_binned_ang.value, emission_binned_ang_parallel.value)
 
 
-@given(model(), time(), nside(), angles(), random_freq(bandpass=True), data())
+@given(
+    model(extrapolate=True),
+    time(),
+    nside(),
+    angles(),
+    random_freq(bandpass=True),
+    data(),
+)
 @settings(deadline=None)
 def test_bandpass_integration(
     model: Zodipy,
@@ -344,7 +351,6 @@ def test_bandpass_integration(
 ) -> None:
     """Property test for bandpass integrations."""
     theta, phi = angles
-    model.extrapolate = True
     observer = data.draw(obs(model, time))
     bp_weights = data.draw(weights(freqs))
     emission_binned = model.get_binned_emission_ang(
@@ -359,7 +365,14 @@ def test_bandpass_integration(
     assert emission_binned.shape == (hp.nside2npix(nside),)
 
 
-@given(model(), time(), nside(), angles(), random_freq(bandpass=True), data())
+@given(
+    model(extrapolate=True),
+    time(),
+    nside(),
+    angles(),
+    random_freq(bandpass=True),
+    data(),
+)
 @settings(deadline=None)
 def test_weights(
     model: Zodipy,
@@ -372,7 +385,6 @@ def test_weights(
     """Property test for bandpass weights."""
 
     theta, phi = angles
-    model.extrapolate = True
     observer = data.draw(obs(model, time))
     bp_weights = data.draw(weights(freqs))
 

tests/test_validators.py

@@ -1,4 +1,5 @@
 import astropy.units as u
+from astropy.time import Time
 import numpy as np
 import pytest
 from hypothesis import given
@@ -13,6 +14,7 @@
 BANDPASS_FREQUENCIES = np.linspace(95, 105, 11) * u.GHz
 BANDPASS_WAVELENGTHS = np.linspace(20, 25, 11) * u.micron
 BANDPASS_WEIGHTS = np.array([2, 3, 5, 9, 11, 12, 11, 9, 5, 3, 2])
+OBS_TIME = Time("2021-01-01T00:00:00")
 
 
 @given(model(extrapolate=False))
@@ -102,3 +104,36 @@ def test_validate_weights_shape() -> None:
     )
     assert weights.size == 1
     assert weights == np.array([1.0], dtype=np.float64)
+
+
+def test_extrapolate_raises_error() -> None:
+    with pytest.raises(ValueError):
+        model = Zodipy("dirbe")
+        model.get_emission_pix(
+            400 * u.micron, pixels=[1, 4, 5], nside=32, obs_time=OBS_TIME
+        )
+
+    model = Zodipy("dirbe", extrapolate=True)
+    model.get_emission_pix(
+        400 * u.micron, pixels=[1, 4, 5], nside=32, obs_time=OBS_TIME
+    )
+
+
+def test_interp_kind() -> None:
+    model = Zodipy("dirbe", interp_kind="linear")
+    linear = model.get_emission_pix(
+        27 * u.micron, pixels=[1, 4, 5], nside=32, obs_time=OBS_TIME
+    )
+
+    model = Zodipy("dirbe", interp_kind="quadratic")
+    quadratic = model.get_emission_pix(
+        27 * u.micron, pixels=[1, 4, 5], nside=32, obs_time=OBS_TIME
+    )
+
+    assert not np.allclose(linear, quadratic)
+
+    with pytest.raises(NotImplementedError):
+        model = Zodipy("dirbe", interp_kind="sdfs")
+        model.get_emission_pix(
+            27 * u.micron, pixels=[1, 4, 5], nside=32, obs_time=OBS_TIME
+        )

zodipy/_interpolate_source.py

@@ -5,7 +5,6 @@
 
 import astropy.units as u
 import numpy as np
-from scipy.interpolate import interp1d
 
 from zodipy._bandpass import Bandpass
 from zodipy._constants import SPECIFIC_INTENSITY_UNITS
@@ -19,12 +18,13 @@
 """Return the source parameters for a given bandpass and model. 
 Must match arguments in the emission fns."""
 GetSourceParametersFn = Callable[
-    [Bandpass, InterplanetaryDustModelT], Dict[Union[ComponentLabel, str], Any]
+    [Bandpass, InterplanetaryDustModelT, Callable],
+    Dict[Union[ComponentLabel, str], Any],
 ]
 
 
 def get_source_parameters_kelsall_comp(
-    bandpass: Bandpass, model: Kelsall
+    bandpass: Bandpass, model: Kelsall, interpolator: Callable
 ) -> dict[ComponentLabel | str, dict[str, Any]]:
     if not bandpass.frequencies.unit.is_equivalent(model.spectrum.unit):
         bandpass.switch_convention()
@@ -35,7 +35,7 @@ def get_source_parameters_kelsall_comp(
         else model.spectrum.to_value(u.micron)
     )
 
-    interpolator = partial(interp1d, x=spectrum, fill_value="extrapolate")
+    interpolator = partial(interpolator, x=spectrum)
 
     source_parameters: dict[ComponentLabel | str, dict[str, Any]] = {}
     for comp_label in model.comps:
@@ -92,7 +92,7 @@ def get_source_parameters_kelsall_comp(
 
 
 def get_source_parameters_rmm(
-    bandpass: Bandpass, model: RRM
+    bandpass: Bandpass, model: RRM, interpolator: Callable
 ) -> dict[ComponentLabel | str, dict[str, Any]]:
     if not bandpass.frequencies.unit.is_equivalent(model.spectrum.unit):
         bandpass.switch_convention()
@@ -104,7 +104,7 @@ def get_source_parameters_rmm(
     )
 
     source_parameters: dict[ComponentLabel | str, dict[str, Any]] = {}
-    calibration = interp1d(x=spectrum, y=model.calibration, fill_value="extrapolate")(
+    calibration = interpolator(x=spectrum, y=model.calibration)(
         bandpass.frequencies.value
     )
     calibration = u.Quantity(calibration, u.MJy / u.AU).to_value(u.Jy / u.cm)

zodipy/zodipy.py

@@ -9,6 +9,7 @@
 import healpy as hp
 import numpy as np
 from astropy.coordinates import solar_system_ephemeris
+from scipy.interpolate import interp1d
 
 from zodipy._bandpass import get_bandpass_interpolation_table, validate_and_get_bandpass
 from zodipy._constants import SPECIFIC_INTENSITY_UNITS
@@ -44,14 +45,18 @@ class Zodipy:
             Defaults to DIRBE.
         gauss_quad_degree (int): Order of the Gaussian-Legendre quadrature used to evaluate
             the line-of-sight integral in the simulations. Default is 50 points.
+        interp_kind (str): Interpolation kind used to interpolate relevant model parameters.
+            Defaults to 'linear'. For more information on available interpolation methods,
+            please visit the [Scipy documentation](
+            https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html).
         extrapolate (bool): If `True` all spectral quantities in the selected model are
             extrapolated to the requested frequencies or wavelengths. If `False`, an
             exception is raised on requested frequencies/wavelengths outside of the
             valid model range. Default is `False`.
         ephemeris (str): Ephemeris used to compute the positions of the observer and the
             Earth. Defaults to 'de432s', which requires downloading (and caching) a ~10MB
-            file. For more information on available ephemeridis, please visit
-            https://docs.astropy.org/en/stable/coordinates/solarsystem.html
+            file. For more information on available ephemeridis, please visit the [Astropy
+            documentation](https://docs.astropy.org/en/stable/coordinates/solarsystem.html)
         solar_cut (u.Quantity[u.deg]): Cutoff angle from the sun in degrees. The emission
             for all the pointing with angular distance between the sun smaller than
             `solar_cutoff` are masked. Defaults to `None`.
@@ -69,6 +74,7 @@ def __init__(
         model: str = "dirbe",
         gauss_quad_degree: int = 50,
         extrapolate: bool = False,
+        interp_kind: str = "linear",
         ephemeris: str = "de432s",
         solar_cut: u.Quantity[u.deg] | None = None,
         solar_cut_fill_value: float = np.nan,
@@ -78,12 +84,18 @@ def __init__(
         self.model = model
         self.gauss_quad_degree = gauss_quad_degree
         self.extrapolate = extrapolate
+        self.interp_kind = interp_kind
         self.ephemeris = ephemeris
         self.solar_cut = solar_cut.to(u.rad) if solar_cut is not None else solar_cut
         self.solar_cut_fill_value = solar_cut_fill_value
         self.parallel = parallel
         self.n_proc = n_proc
 
+        self._interpolator = partial(
+            interp1d,
+            kind=self.interp_kind,
+            fill_value="extrapolate" if self.extrapolate else np.nan,
+        )
         self._ipd_model = model_registry.get_model(model)
         self._gauss_points_and_weights = np.polynomial.legendre.leggauss(
             gauss_quad_degree
@@ -414,7 +426,7 @@ def _compute_emission(
         # Get model parameters, some of which have been interpolated to the given
         # frequency or bandpass.
         source_parameters = SOURCE_PARAMS_MAPPING[type(self._ipd_model)](
-            bandpass, self._ipd_model
+            bandpass, self._ipd_model, self._interpolator
         )
 
         observer_position, earth_position = get_obs_and_earth_positions(