Merge pull request #25 from Cosmoglobe/astropy_healpix

Replace healpy with astropy-healpi. Add new API to directly support SkyCoord objects. Move freq and weights arguments to Zodipy object initialization and away from methods

.github/workflows/draft-pdf.yml

@@ -1,4 +1,9 @@
-on: [push, workflow_dispatch]
+on:
+  push:
+    paths: ["paper/**"]
+  pull_request:
+    paths: ["paper/**"]
+  workflow_dispatch:
 
 jobs:
   paper:
@@ -20,4 +25,4 @@ jobs:
           # This is the output path where Pandoc will write the compiled
           # PDF. Note, this should be the same directory as the input
           # paper.md
-          path: paper/paper.pdf
+          path: paper/paper.pdf

.pre-commit-config.yaml

@@ -7,12 +7,17 @@ repos:
       langauge_version: python3
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    # Ruff version.
     rev: v0.3.7
     hooks:
       - id: ruff
       - id: ruff-format
 
+  - repo: https://github.com/pre-commit/mirrors-mypy
+    rev: v1.10.0
+    hooks:
+    -   id: mypy
+        exclude: tests/
+        additional_dependencies: [numpy>=1.21.0]
 hooks:
   - id: poetry-export
     args: ["--dev", "-f", "requirements.txt", "-o", "requirements-dev.txt"]

README.md

@@ -57,8 +57,8 @@ ZodiPy is installed using `pip install zodipy`.
 ## Dependencies
 ZodiPy supports all Python versions >= 3.9, and has the following dependencies:
 - [Astropy](https://www.astropy.org/) (>=5.0.1)
+- [Astropy-healpix](https://astropy-healpix.readthedocs.io/en/latest/)
 - [NumPy](https://numpy.org/)
-- [healpy](https://healpy.readthedocs.io/en/latest/)
 - [jplephem](https://pypi.org/project/jplephem/)
 - [SciPy](https://scipy.org/)
 
@@ -67,6 +67,7 @@ Contributing developers will need to download the following additional dependenc
 - pytest
 - pytest-cov
 - hypothesis
+- healpy
 - coverage
 - ruff
 - mypy

docs/examples/get_bandpass_integrated_emission.py

@@ -1,3 +1,5 @@
+from multiprocessing import cpu_count
+
 import astropy.units as u
 import healpy as hp
 import matplotlib.pyplot as plt
@@ -15,16 +17,15 @@
 plt.plot(freqs, weights)
 plt.xlabel("Frequency [GHz]")
 plt.ylabel("Weights")
-plt.savefig("../img/bandpass.png", dpi=300)
 
-model = Zodipy(model="planck18")
+model = Zodipy(model="planck18", n_proc=cpu_count())
 
 emission_central_freq = model.get_binned_emission_pix(
     freq=center_freq,
     pixels=np.arange(hp.nside2npix(nside)),
     nside=nside,
     obs_time=Time("2022-03-10"),
-    obs="SEMB-L2",
+    obs_pos="SEMB-L2",
 )
 
 emission_bandpass_integrated = model.get_binned_emission_pix(
@@ -33,7 +34,7 @@
     pixels=np.arange(hp.nside2npix(nside)),
     nside=nside,
     obs_time=Time("2022-03-10"),
-    obs="SEMB-L2",
+    obs_pos="SEMB-L2",
 )
 
 hp.mollview(
@@ -43,7 +44,6 @@
     cmap="afmhot",
     norm="log",
 )
-plt.savefig("../img/center_freq.png", dpi=300)
 
 hp.mollview(
     emission_bandpass_integrated,
@@ -52,5 +52,4 @@
     cmap="afmhot",
     norm="log",
 )
-plt.savefig("../img/bandpass_integrated.png", dpi=300)
 plt.show()

docs/examples/get_binned_emission.py

@@ -1,3 +1,5 @@
+from multiprocessing import cpu_count
+
 import astropy.units as u
 import healpy as hp
 import matplotlib.pyplot as plt
@@ -6,15 +8,15 @@
 
 from zodipy import Zodipy
 
-model = Zodipy("planck18")
+model = Zodipy("planck18", n_proc=cpu_count())
 nside = 256
 
 binned_emission = model.get_binned_emission_pix(
-    857 * u.GHz,
-    pixels=np.arange(hp.nside2npix(nside)),
+    np.arange(hp.nside2npix(nside)),
+    freq=857 * u.GHz,
     nside=nside,
     obs_time=Time("2022-06-14"),
-    obs="earth",
+    obs_pos="earth",
 )
 
 hp.mollview(
@@ -25,5 +27,4 @@
     max=1,
     cmap="afmhot",
 )
-plt.savefig("../img/binned.png", dpi=300)
 plt.show()

docs/examples/get_binned_emission_solar_cutoff.py

@@ -1,3 +1,5 @@
+from multiprocessing import cpu_count
+
 import astropy.units as u
 import healpy as hp
 import matplotlib.pyplot as plt
@@ -6,15 +8,16 @@
 
 from zodipy import Zodipy
 
-model = Zodipy("dirbe", solar_cut=60 * u.deg)
+model = Zodipy("dirbe", n_proc=cpu_count())
 nside = 256
 
 binned_emission = model.get_binned_emission_pix(
-    25 * u.micron,
-    pixels=np.arange(hp.nside2npix(nside)),
+    np.arange(hp.nside2npix(nside)),
     nside=nside,
+    freq=25 * u.micron,
     obs_time=Time("2020-01-01"),
-    obs="earth",
+    obs_pos="earth",
+    solar_cut=60 * u.deg,
 )
 
 hp.mollview(
@@ -25,5 +28,4 @@
     coord="E",
     cmap="afmhot",
 )
-plt.savefig("../img/binned_solar_cutoff.png", dpi=300)
 plt.show()

docs/examples/get_binned_gal_emission.py

@@ -1,3 +1,5 @@
+from multiprocessing import cpu_count
+
 import astropy.units as u
 import healpy as hp
 import matplotlib.pyplot as plt
@@ -6,16 +8,16 @@
 
 from zodipy import Zodipy
 
-model = Zodipy("planck18")
+model = Zodipy("planck18", n_proc=cpu_count())
 nside = 256
 
 binned_emission = model.get_binned_emission_pix(
-    857 * u.GHz,
-    pixels=np.arange(hp.nside2npix(nside)),
+    np.arange(hp.nside2npix(nside)),
+    freq=857 * u.GHz,
     nside=nside,
     obs_time=Time("2022-02-20"),
-    obs="earth",
-    coord_in="G",  # Coordinates of the input pointing
+    obs_pos="earth",
+    frame="galactic",  # Coordinates of the input pointing
 )
 
 hp.mollview(
@@ -26,5 +28,4 @@
     min=0,
     max=1,
 )
-plt.savefig("../img/binned_gal.png", dpi=300)
 plt.show()

docs/examples/get_comp_binned_emission.py

@@ -1,3 +1,5 @@
+from multiprocessing import cpu_count
+
 import astropy.units as u
 import healpy as hp
 import matplotlib.pyplot as plt
@@ -6,15 +8,15 @@
 
 from zodipy import Zodipy
 
-model = Zodipy("dirbe")
+model = Zodipy("dirbe", n_proc=cpu_count())
 nside = 256
 
 binned_emission = model.get_binned_emission_pix(
-    25 * u.micron,
-    pixels=np.arange(hp.nside2npix(nside)),
+    np.arange(hp.nside2npix(nside)),
+    freq=25 * u.micron,
     nside=nside,
     obs_time=Time("2022-01-01"),
-    obs="earth",
+    obs_pos="earth",
     return_comps=True,
 )
 fig = plt.figure(figsize=(8, 6.5), constrained_layout=True)
@@ -29,5 +31,4 @@
         sub=(3, 2, idx + 1),
         fig=fig,
     )
-# plt.savefig("../img/binned_comp.png", dpi=300)
 plt.show()

docs/examples/get_density_contour.py

@@ -26,5 +26,4 @@
 plt.title("Cross section of the interplanetary dust density (yz-plane)")
 plt.xlabel("x [AU]")
 plt.ylabel("z [AU]")
-# plt.savefig("../img/density_grid.png", dpi=300)
 plt.show()

docs/examples/get_emission_ang.py

@@ -1,27 +1,30 @@
 import astropy.units as u
 import matplotlib.pyplot as plt
 import numpy as np
+from astropy.coordinates import BarycentricMeanEcliptic, SkyCoord
 from astropy.time import Time
 
 from zodipy import Zodipy
 
-model = Zodipy("dirbe")
+model = Zodipy()
 
-latitudes = np.linspace(-90, 90, 10000) * u.deg
-longitudes = np.zeros_like(latitudes)
+# Longitude and Latitude values corresponding to a scan through the eclitpic plane
+lats = np.linspace(-90, 90, 100) * u.deg
+lons = np.zeros_like(lats)
 
-emission = model.get_emission_ang(
-    30 * u.micron,
-    theta=longitudes,
-    phi=latitudes,
-    lonlat=True,
-    obs_time=Time("2022-06-14"),
-    obs="earth",
+obs_time = Time("2022-06-14")
+
+# The SkyCoord object needs to include the coordinate frame and time of observation
+coords = SkyCoord(
+    lons,
+    lats,
+    frame=BarycentricMeanEcliptic,
+    obstime=obs_time,
 )
 
+emission = model.get_emission_skycoord(coords, freq=30 * u.micron)
 
-plt.plot(latitudes, emission)
+plt.plot(lats, emission)
 plt.xlabel("Latitude [deg]")
 plt.ylabel("Emission [MJy/sr]")
-plt.savefig("../img/timestream.png", dpi=300)
 plt.show()

docs/examples/get_parallel_emission.py

@@ -1,4 +1,5 @@
 import time
+from multiprocessing import cpu_count
 
 import astropy.units as u
 import healpy as hp
@@ -10,15 +11,15 @@
 nside = 256
 pixels = np.arange(hp.nside2npix(nside))
 obs_time = Time("2020-01-01")
-n_proc = 8
+n_proc = cpu_count()
 
 model = Zodipy()
 model_parallel = Zodipy(n_proc=n_proc)
 
 start = time.perf_counter()
 emission = model.get_binned_emission_pix(
-    40 * u.micron,
     pixels=pixels,
+    freq=40 * u.micron,
     nside=nside,
     obs_time=obs_time,
 )
@@ -27,9 +28,9 @@
 
 start = time.perf_counter()
 emission_parallel = model_parallel.get_binned_emission_pix(
-    40 * u.micron,
     pixels=pixels,
     nside=nside,
+    freq=40 * u.micron,
     obs_time=obs_time,
 )
 print(

docs/install.md

@@ -13,7 +13,7 @@ pip install zodipy
 ZodiPy has the following dependencies (these are automatically downloaded alongside ZodiPy):
 
 - [Astropy](https://www.astropy.org) (>= 5.0.1)
+- [Astropy-healpix](https://astropy-healpix.readthedocs.io/en/latest/)
 - [NumPy](https://numpy.org)
-- [healpy](https://healpy.readthedocs.io/en/latest/)
 - [jplehem](https://pypi.org/project/jplephem/)
 - [SciPy](https://scipy.org/)

docs/usage.md

@@ -49,10 +49,9 @@ Note that the bandpass weights must be in spectral radiance units (Jy/sr), even
 ![Center frequency emission](img/center_freq.png)
 ![Bandpass integrated emission](img/bandpass_integrated.png)
 
-### Solar cutoff angle
-Few experiments look directly in towards the Sun. We can initialize `Zodipy` with the `solar_cut` 
-argument to mask all input pointing that looks in towards the sun with an angular distance smaller 
-than the `solar_cut` value.
+### Solar cutoff
+In the case where the user wishes to bin a set of observations into a HEALPix map, using the `*_binned_*` API functions,
+a solar cutoff angle can be specified to mask out the all pixels closer to the Sun by some angle given by `solar_cut`
 
 ```python hl_lines="9"
 {!examples/get_binned_emission_solar_cutoff.py!}