Problems found with Gammapy v0.19

[monicava]

Here are the first 4 problems I found using v0.19

1. Plotting edips gives empty plots

   for ids in obs_ids:
   obs = datastore.obs(ids)
    obs.edisp.peek()

See difference between 0.18.2 and 0.19 versions:

2. Theta2 empty plot

position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
theta2_axis = MapAxis.from_bounds(0, 0.2, nbin=20, interp="lin", unit="deg2")
theta2_table = make_theta_squared_table(
    observations=observations,
    position=position,
    theta_squared_axis=theta2_axis,
)

3. Masking gives the error AttributeError: 'WcsNDMap' object has no attribute 'shape'

exclusion_region = CircleSkyRegion(
    center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
    radius=0.5 * u.deg,
)
skydir = target_position.galactic
exclusion_mask = Map.create(
    npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
)
mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)
exclusion_mask.data = mask
exclusion_mask.plot();

---

AttributeError Traceback (most recent call last)
/var/folders/y2/ty2m_5q50bsdx5hq513k71pr0000gn/T/ipykernel_21610/2571132427.py in
10
11 mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)
---> 12 exclusion_mask.data = mask
13 exclusion_mask.plot();

~/opt/anaconda3/envs/gammapy-0.19/lib/python3.9/site-packages/gammapy/maps/core.py in data(self, value)
93 raise TypeError("Map data must be a Numpy array. Set unit separately")
94
---> 95 if not value.shape == self.geom.data_shape:
96 value = value.reshape(self.geom.data_shape)
97

AttributeError: 'WcsNDMap' object has no attribute 'shape'

4. Data reduction first step gives an error (see below)

e_reco = MapAxis.from_energy_bounds(0.1, 40, 40, unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(
    0.05, 100, 80, unit="TeV", name="energy_true"
)
dataset_empty = SpectrumDataset.create(
    e_reco=e_reco, e_true=e_true, region=on_region
)

---

TypeError Traceback (most recent call last)
/var/folders/y2/ty2m_5q50bsdx5hq513k71pr0000gn/T/ipykernel_21610/786087746.py in
3 0.05, 100, 80, unit="TeV", name="energy_true"
4 )
----> 5 dataset_empty = SpectrumDataset.create(
6 e_reco=e_reco, e_true=e_true, region=on_region
7 )

TypeError: create() missing 1 required positional argument: 'geom'

[registerrier]

Hi @monicava ,

3. Masking gives the error AttributeError: 'WcsNDMap' object has no attribute 'shape'

exclusion_region = CircleSkyRegion(
    center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
    radius=0.5 * u.deg,
)
skydir = target_position.galactic
exclusion_mask = Map.create(
    npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
)
mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)
exclusion_mask.data = mask
exclusion_mask.plot();

AttributeError Traceback (most recent call last) /var/folders/y2/ty2m_5q50bsdx5hq513k71pr0000gn/T/ipykernel_21610/2571132427.py in 10 11 mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False) ---> 12 exclusion_mask.data = mask 13 exclusion_mask.plot();

~/opt/anaconda3/envs/gammapy-0.19/lib/python3.9/site-packages/gammapy/maps/core.py in data(self, value) 93 raise TypeError("Map data must be a Numpy array. Set unit separately") 94 ---> 95 if not value.shape == self.geom.data_shape: 96 value = value.reshape(self.geom.data_shape) 97

AttributeError: 'WcsNDMap' object has no attribute 'shape'

This issue is due to the fact that now geom.region_mask() now returns a Map object. Therefore you code snippet should simplify in v0.19 to:

from gammapy.maps import WcsGeom
exclusion_region = CircleSkyRegion(
        center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
        radius=0.5 * u.deg,
)
geom = WcsGeom.create(
        npix=(150, 150), binsz=0.05, skydir=exclusion_region.center, proj="TAN", frame="icrs"
)
exclusion_mask =  geom.region_mask([exclusion_region], inside=False)
exclusion_mask.plot()

See also the docs on exclusion masks here:
https://docs.gammapy.org/0.19/tutorials/api/mask_maps.html#Creating-an-exclusion-mask

[monicava]

The code works but the x-axis range moves, funny ...

[registerrier]

I suppose this is because I used the exclusion_region.center rather than your target_position.galactic as the center of the geom.

[registerrier]

Regarding problem 4:

4. Data reduction first step gives an error (see below)

e_reco = MapAxis.from_energy_bounds(0.1, 40, 40, unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(
    0.05, 100, 80, unit="TeV", name="energy_true"
)
dataset_empty = SpectrumDataset.create(
    e_reco=e_reco, e_true=e_true, region=on_region
)

TypeError Traceback (most recent call last) /var/folders/y2/ty2m_5q50bsdx5hq513k71pr0000gn/T/ipykernel_21610/786087746.py in 3 0.05, 100, 80, unit="TeV", name="energy_true" 4 ) ----> 5 dataset_empty = SpectrumDataset.create( 6 e_reco=e_reco, e_true=e_true, region=on_region 7 )

TypeError: create() missing 1 required positional argument: 'geom'

Again this is an API change, now the SpectrumDataset.create() method takes a geom in argument for consistency with MapDataset, this translates the code this way:

from gammapy.maps import RegionGeom

e_reco = MapAxis.from_energy_bounds(0.1, 40, 40, unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(
    0.05, 100, 80, unit="TeV", name="energy_true"
)
geom = RegionGeom.create(region=on_region, axes=[energy_axis]) 

This wiki (in progress) also refers to this change:
https://github.com/gammapy/gammapy/wiki/Main-API-changes-from-v0.18.2-to-v0.19#creating-a-spectrumdataset

[registerrier]

Hi @monicava

2. Theta2 empty plot

position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
theta2_axis = MapAxis.from_bounds(0, 0.2, nbin=20, interp="lin", unit="deg2")
theta2_table = make_theta_squared_table(
    observations=observations,
    position=position,
    theta_squared_axis=theta2_axis,
)

I don't know if your two first issues are connected or not.
It is possible that this one is due to the observations skipped while reading in the DataStore. Could you check that observations is not empty?

Now DataStore.get_observations() takes a required_irfs argument set to "all" by default. This might be an issue here. Can you try with

observations = datastore.get_observations(obs_ids, required_irfs = ["aeff", "edisp"])

see:
https://docs.gammapy.org/dev/api/gammapy.data.DataStore.html#gammapy.data.DataStore.get_observations

[monicava]

I confirm it is a single bin:

    print(obs.edisp.axes["offset"])
    print(obs.edisp.axes["offset"].center)

MapAxis

name       : offset    
unit       : 'deg'     
nbins      : 1         
node type  : edges     
edges min  : 3.0e-01 deg
edges max  : 5.0e-01 deg
interp     : lin       

[0.4] deg

[adonath]

And how does the plot look like if you use something like:

import matplotlib.pyplot as plt
from astropy import units as u

offset = 0.4 * u.deg
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=figsize)
edisp.plot_bias(ax=axes[0], offset=offset)
edisp.plot_migration(ax=axes[1], offset=offset)
edisp = self.to_edisp_kernel(offset=offset)
edisp.plot_matrix(ax=axes[2])

[monicava]

I commented the lines below to avoid the "self" not defined

    offset = 0.4 * u.deg
    fig, axes = plt.subplots(nrows=1, ncols=3,figsize=(10, 4))
    obs.edisp.plot_bias(ax=axes[0], offset=offset)
    obs.edisp.plot_migration(ax=axes[1], offset=offset)
    #obs.edisp = self.to_edisp_kernel(offset=offset)
    #obs.edisp.plot_matrix(ax=axes[2])

---> 12 obs.edisp = self.to_edisp_kernel(offset=offset)
13 obs.edisp.plot_matrix(ax=axes[2])
NameError: name 'self' is not defined

I get the following:

[monicava]

I was focusing on issue 1) (Edisp empty plots) but I confirm that using:
observations = datastore.get_observations(obs_ids, required_irf = ["aeff", "edisp"])
the issue 2 (theta2 plot) now appears filled with v0.19 :)

[registerrier]

So here we have a documentation issue, and we need to remove the default option as required_irfs=None.