Merge branch 'gammapy:main' into variability_tutorial

docs/api-reference/utils.rst

@@ -8,6 +8,10 @@ utils - Utilities
 
 .. currentmodule:: gammapy.utils
 
+.. automodapi:: gammapy.utils.cluster
+    :no-inheritance-diagram:
+    :include-all-objects:
+
 .. automodapi:: gammapy.utils.coordinates
     :no-inheritance-diagram:
     :include-all-objects:
@@ -32,6 +36,10 @@ utils - Utilities
     :no-inheritance-diagram:
     :include-all-objects:
 
+.. automodapi:: gammapy.utils.parallel
+    :no-inheritance-diagram:
+    :include-all-objects:
+
 .. automodapi:: gammapy.utils.scripts
     :no-inheritance-diagram:
     :include-all-objects:

docs/development/pigs/pig-021.rst

@@ -97,7 +97,7 @@ An example of how this can be achieved with a custom implemented model is given
 .. code::
 
     from gammapy.modeling.models import SpatialModel
-    from astropy.coordinates.angle_utilities import angular_separation
+    from astropy.coordinates import angular_separation
 
     class MyCustomGaussianModel(SpatialModel):
         """My custom gaussian model.

docs/user-guide/estimators.rst

@@ -27,8 +27,9 @@ In general the flux can be estimated using two methods:
 #. **Based on model fitting:** given a (global) best fit model with multiple model components,
    the flux of the component of interest is re-fitted in the chosen energy, time or spatial
    region. The new flux is given as a ``norm`` with respect to the global reference model.
-   Optionally other component parameters in the global model can be re-optimised. This method
-   is also named **forward folding**.
+   Optionally the free parameters of the other models can be re-optimised
+   (but the other parameters of the source of interest are always kept frozen).
+   This method is also named **forward folding**.
 
 #. **Based on excess:** in the case of having one energy bin, neglecting the PSF and
    not re-optimising other parameters, one can estimate the significance based on the

docs/user-guide/howto.rst

@@ -41,8 +41,8 @@ or also select observations based on other information available using the `~gam
 
 .. accordion-header::
     :id: collapseHowToThree
-    :title: Check IRFs
-    :link: ../tutorials/data/cta.html#irfs
+    :title: Make observation duration maps
+    :link: ../tutorials/api/makers.html#observation-duration-and-effective-livetime
 
 Gammapy offers a number of methods to explore the content of the various IRFs
 contained in an observation. This is usually done thanks to their ``peek()``
@@ -373,3 +373,14 @@ using a dummy phase column.
     events_new.write("events.fits.gz", gti=obs.gti, overwrite=True)
 
 .. accordion-footer::
+
+.. accordion-header::
+    :id: collapseHowTo_ObservationMaps
+    :title: Make observation duration maps
+    :link: ../tutorials/api/makers.html#observation-duration-and-effective-livetime
+
+Compute the absolute and acceptance corrected observation duration in the field of view
+
+.. accordion-footer::
+
+

environment-dev.yml

@@ -48,7 +48,7 @@ dependencies:
   - reproject
   - sherpa
   # dev dependencies
-  - black
+  - black=22.6.0
   - codespell
   - flake8
   - isort

examples/models/spectral/plot_absorbed.py

@@ -16,7 +16,7 @@
 redshift z of the source, and :math:`\alpha` is a scale factor
 (default: 1) for the optical depth.
 
-The available EBL models are defined in `~gammapy.modeling.models.spectral.EBL_DATA_BUILTIN`.
+The available EBL models are defined in `~gammapy.modeling.models.EBL_DATA_BUILTIN`.
 """
 
 # %%
@@ -27,12 +27,12 @@
 from astropy import units as u
 import matplotlib.pyplot as plt
 from gammapy.modeling.models import (
+    EBL_DATA_BUILTIN,
     EBLAbsorptionNormSpectralModel,
     Models,
     PowerLawSpectralModel,
     SkyModel,
 )
-from gammapy.modeling.models.spectral import EBL_DATA_BUILTIN
 
 # Print the available EBL models
 print(EBL_DATA_BUILTIN.keys())

examples/tutorials/analysis-3d/cta_data_analysis.py

@@ -2,7 +2,7 @@
 Basic image exploration and fitting
 ===================================
 
-Detect sources, produce a sky image and a spectrum using CTA 1DC data.
+Detect sources, produce a sky image and a spectrum using CTA-1DC data.
 
 Introduction
 ------------
@@ -11,7 +11,7 @@
 for simulated CTA data with Gammapy.**
 
 The dataset we will use is three observation runs on the Galactic
-center. This is a tiny (and thus quick to process and play with and
+Center. This is a tiny (and thus quick to process and play with and
 learn) subset of the simulated CTA dataset that was produced for the
 first data challenge in August 2017.
 
@@ -32,8 +32,6 @@
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from regions import CircleSkyRegion
-
-# %matplotlib inline
 import matplotlib.pyplot as plt
 from IPython.display import display
 from gammapy.data import DataStore
@@ -74,8 +72,8 @@
 # A Gammapy analysis usually starts by creating a
 # `~gammapy.data.DataStore` and selecting observations.
 #
-# This is shown in detail in the other notebook, here we just pick three
-# observations near the galactic center.
+# This is shown in detail in other notebooks (see e.g. the :doc:`/tutorials/starting/analysis_2` tutorial),
+# here we choose three observations near the Galactic Center.
 #
 
 data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps")
@@ -108,8 +106,19 @@
 # analysis
 #
 
-axis = MapAxis.from_edges(
-    np.logspace(-1.0, 1.0, 10), unit="TeV", name="energy", interp="log"
+axis = MapAxis.from_energy_bounds(
+    0.1,
+    10,
+    nbin=10,
+    unit="TeV",
+    name="energy",
+)
+axis_true = MapAxis.from_energy_bounds(
+    0.05,
+    20,
+    nbin=20,
+    name="energy_true",
+    unit="TeV",
 )
 geom = WcsGeom.create(
     skydir=(0, 0), npix=(500, 400), binsz=0.02, frame="galactic", axes=[axis]
@@ -123,8 +132,7 @@
 #
 
 # %%time
-stacked = MapDataset.create(geom=geom)
-stacked.edisp = None
+stacked = MapDataset.create(geom=geom, energy_axis_true=axis_true)
 maker = MapDatasetMaker(selection=["counts", "background", "exposure", "psf"])
 maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=2.5 * u.deg)
 
@@ -329,7 +337,7 @@
 
 
 ######################################################################
-# Call `plot_npred_signal` to plot the predicted counts.
+# Call `~gammapy.visualization.plot_npred_signal` to plot the predicted counts.
 #
 
 
@@ -361,7 +369,7 @@
 #
 # Let’s plot the spectral model and points. You could do it directly, but
 # for convenience we bundle the model and the flux points in a
-# `FluxPointDataset`:
+# `~gammapy.datasets.FluxPointsDataset`:
 #
 
 flux_points_dataset = FluxPointsDataset(data=flux_points, models=model)

examples/tutorials/api/datasets.py

@@ -145,7 +145,7 @@
 
 ######################################################################
 # To further explore the contents of a `Dataset`, you can use
-# e.g. `~gammapy.datasets.MapDataset.info_dict()`
+# e.g. `~gammapy.datasets.MapDataset.info_dict()`
 #
 
 ######################################################################
@@ -171,7 +171,7 @@
 
 
 ######################################################################
-# Of course you can also access IRF related maps, e.g. the psf as
+# Of course you can also access IRF related maps, e.g. the psf as
 # `~gammapy.irf.PSFMap`:
 #
 
@@ -235,7 +235,7 @@
 # use:
 #
 
-npred_source = dataset_cta.npred_signal(model_name="gc")
+npred_source = dataset_cta.npred_signal(model_names=["gc"])
 npred_source.sum_over_axes().plot()
 plt.show()
 
@@ -263,7 +263,7 @@
 #    according to the specified selection cuts, and should not be changed
 #    by the user.
 # -  During modelling and fitting, the user might want to additionally
-#    ignore some parts of a reduced dataset, e.g. to restrict the fit to a
+#    ignore some parts of a reduced dataset, e.g. to restrict the fit to a
 #    specific energy range or to ignore parts of the region of interest.
 #    This should be done by applying the `~gammapy.datasets.MapDataset.mask_fit`. To see details of
 #    applying masks, please refer to :ref:`masks-for-fitting`.

examples/tutorials/api/makers.py

@@ -20,6 +20,7 @@
 
 import numpy as np
 from astropy import units as u
+from astropy.coordinates import SkyCoord
 from regions import CircleSkyRegion
 import matplotlib.pyplot as plt
 from IPython.display import display
@@ -33,6 +34,7 @@
     SafeMaskMaker,
     SpectrumDatasetMaker,
 )
+from gammapy.makers.utils import make_effective_livetime_map, make_observation_time_map
 from gammapy.maps import MapAxis, RegionGeom, WcsGeom
 
 ######################################################################
@@ -234,7 +236,7 @@
 # This method is only used for 1D spectral analysis.
 #
 # For more details and usage, see
-# the :doc:`Reflected background </user-guide/makers/ring>`
+# the :doc:`Reflected background </user-guide/makers/reflected>`
 #
 # Data reduction loop
 # -------------------
@@ -342,3 +344,95 @@
 print(datasets)
 
 plt.show()
+
+######################################################################
+# Observation duration and effective livetime
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# It can often be useful to know the total number of hours spent
+# in the given field of view (without correcting for the acceptance
+# variation). This can be computed using `make_observation_time_map`
+# as shown below
+#
+
+# Get the observations
+obs_id = data_store.obs_table["OBS_ID"][data_store.obs_table["OBJECT"] == "MSH 15-5-02"]
+observations = data_store.get_observations(obs_id)
+print("No. of observations: ", len(observations))
+
+# Define an energy range
+energy_min = 100 * u.GeV
+energy_max = 10.0 * u.TeV
+
+# Define an offset cut (the camera field of view)
+offset_max = 2.5 * u.deg
+
+# Define the geom
+source_pos = SkyCoord(228.32, -59.08, unit="deg")
+energy_axis_true = MapAxis.from_energy_bounds(
+    energy_min, energy_max, nbin=2, name="energy_true"
+)
+geom = WcsGeom.create(
+    skydir=source_pos,
+    binsz=0.02,
+    width=(6, 6),
+    frame="icrs",
+    proj="CAR",
+    axes=[energy_axis_true],
+)
+
+total_obstime = make_observation_time_map(observations, geom, offset_max=offset_max)
+
+
+plt.figure(figsize=(5, 5))
+ax = total_obstime.plot(add_cbar=True)
+# Add the pointing position on top
+for obs in observations:
+    ax.plot(
+        obs.get_pointing_icrs(obs.tmid).to_pixel(wcs=ax.wcs)[0],
+        obs.get_pointing_icrs(obs.tmid).to_pixel(wcs=ax.wcs)[1],
+        "+",
+        color="black",
+    )
+ax.set_title("Total observation time")
+plt.show()
+
+######################################################################
+# As the acceptance of IACT cameras vary within the field of
+# view, it can also be interesting to plot the on-axis equivalent
+# number of hours.
+#
+
+effective_livetime = make_effective_livetime_map(
+    observations, geom, offset_max=offset_max
+)
+
+
+axs = effective_livetime.plot_grid(add_cbar=True)
+# Add the pointing position on top
+for ax in axs:
+    for obs in observations:
+        ax.plot(
+            obs.get_pointing_icrs(obs.tmid).to_pixel(wcs=ax.wcs)[0],
+            obs.get_pointing_icrs(obs.tmid).to_pixel(wcs=ax.wcs)[1],
+            "+",
+            color="black",
+        )
+plt.show()
+
+######################################################################
+# To get the value of the observation time at a particular position,
+# use `get_by_coord`
+
+obs_time_src = total_obstime.get_by_coord(source_pos)
+effective_times_src = effective_livetime.get_by_coord(
+    (source_pos, energy_axis_true.center)
+)
+
+print(f"Time spent on position {source_pos}")
+print(f"Total observation time: {obs_time_src}* {total_obstime.unit}")
+print(
+    f"Effective livetime at {energy_axis_true.center[0]}: {effective_times_src[0]} * {effective_livetime.unit}"
+)
+print(
+    f"Effective livetime at {energy_axis_true.center[1]}: {effective_times_src[1]} * {effective_livetime.unit}"
+)

examples/tutorials/api/models.py

@@ -760,7 +760,7 @@ def evaluate(energy, index, amplitude, reference, mean, width):
 # models.
 #
 
-from astropy.coordinates.angle_utilities import angular_separation
+from astropy.coordinates import angular_separation
 from gammapy.modeling.models import SpatialModel