From b075b114545b59e81d2c6451393925529e5d3cb9 Mon Sep 17 00:00:00 2001 From: Metin San Date: Mon, 6 Feb 2023 15:30:23 +0100 Subject: [PATCH 1/4] Add new keyword `interp_kind` to `Zodipy` for controlling interpolation kind. Previously only linear interpolation was performed for the relevant source and spectral parameters in the models. Now the user may select between the various supported methods in https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html. --- tests/test_get_emission.py | 20 ++++++++++++++++---- tests/test_validators.py | 35 +++++++++++++++++++++++++++++++++++ zodipy/_interpolate_source.py | 12 ++++++------ zodipy/zodipy.py | 14 +++++++++++++- 4 files changed, 70 insertions(+), 11 deletions(-) diff --git a/tests/test_get_emission.py b/tests/test_get_emission.py index 7997952..6e276e4 100644 --- a/tests/test_get_emission.py +++ b/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)) diff --git a/tests/test_validators.py b/tests/test_validators.py index 8533b21..b4fbc5a 100644 --- a/tests/test_validators.py +++ b/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 + ) diff --git a/zodipy/_interpolate_source.py b/zodipy/_interpolate_source.py index 0902def..ee584d1 100644 --- a/zodipy/_interpolate_source.py +++ b/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) diff --git a/zodipy/zodipy.py b/zodipy/zodipy.py index 29b887b..0082994 100644 --- a/zodipy/zodipy.py +++ b/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,6 +45,10 @@ 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 the model. Defaults to + 'linear'. For more information on available interpolation methods, please + visit + 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 @@ -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( From 56fdefe48b7739a7674b6e6d77a805a34f81d983 Mon Sep 17 00:00:00 2001 From: Metin San Date: Mon, 6 Feb 2023 15:42:47 +0100 Subject: [PATCH 2/4] Add hyperlinks to docstrings --- zodipy/zodipy.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/zodipy/zodipy.py b/zodipy/zodipy.py index 0082994..92bf78b 100644 --- a/zodipy/zodipy.py +++ b/zodipy/zodipy.py @@ -45,18 +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 the model. Defaults to - 'linear'. For more information on available interpolation methods, please - visit - https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html + 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`. From 15a9c0450faac6e8d62e167e4288c80bad6e5040 Mon Sep 17 00:00:00 2001 From: Metin San Date: Mon, 6 Feb 2023 16:23:50 +0100 Subject: [PATCH 3/4] Bump version number --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index a147008..71b802f 100644 --- a/pyproject.toml +++ b/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 "] readme = "README.md" From 8a054791651b299ca18c4dc0a59829e28a0f0c67 Mon Sep 17 00:00:00 2001 From: Metin San Date: Mon, 6 Feb 2023 17:05:17 +0100 Subject: [PATCH 4/4] Update documentation --- docs/index.md | 14 +++++++------- docs/introduction.md | 12 ++++++------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/index.md b/docs/index.md index 56b80b2..e7f9aa8 100644 --- a/docs/index.md +++ b/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). diff --git a/docs/introduction.md b/docs/introduction.md index ab2f198..b69ae36 100644 --- a/docs/introduction.md +++ b/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) \ No newline at end of file