Repository: gammapy/gammapy

Issue:
plotting `ConstantSpatialModel`
`Models.plot_regions()/plot_positions()` fails if it contains any `ConstantSpatialModel` (or `ConstantFluxSpatialModel`), with a rather complicated error message. I suppose this is because there is no `position` defined for these models. 

It would be nice to have this fixed before the 1.0rc release.

A quick solution can be to simply   add an exception for `ValueError ` as well here https://github.com/gammapy/gammapy/blob/a68be851d6523742627ed293aa1122c795680f8e/gammapy/modeling/models/core.py#L918.

Discussion:

Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0e69a2253..d29f58780 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -131,7 +131,7 @@ jobs:
           make docs-sphinx
       - name: test code in docs and docstrings
         run: |
-          pytest --doctest-glob="*.rst" docs/
+          pytest --doctest-glob="*.rst" --ignore-glob=docs/development/pigs/*.rst docs/
           pytest --doctest-modules --ignore-glob=*/tests gammapy
       - name: test without datasets
         run: |
diff --git a/docs/tutorials/api/models.ipynb b/docs/tutorials/api/models.ipynb
index 612bea666..71f5e5281 100644
--- a/docs/tutorials/api/models.ipynb
+++ b/docs/tutorials/api/models.ipynb
@@ -632,17 +632,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from regions import write_ds9\n",
+    "from regions import Regions\n",
     "\n",
-    "regions = [gauss.to_region(), gauss_elongated.to_region()]\n",
+    "regions = Regions([gauss.to_region(), gauss_elongated.to_region()])\n",
     "\n",
     "filename = \"regions.reg\"\n",
-    "write_ds9(\n",
-    "    regions,\n",
+    "regions.write(\n",
     "    filename,\n",
-    "    coordsys=\"galactic\",\n",
-    "    fmt=\".4f\",\n",
-    "    radunit=\"deg\",\n",
+    "    format=\"ds9\",\n",
     "    overwrite=True,\n",
     ")"
    ]
@@ -1511,7 +1508,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.12"
+   "version": "3.9.12"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
diff --git a/environment-dev.yml b/environment-dev.yml
index dc8f5bb33..45fafc097 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -24,7 +24,7 @@ dependencies:
   - numpy
   - pydantic
   - pyyaml
-  - regions=0.5
+  - regions>=0.5
   - matplotlib>=3.4
   - scipy
   # test dependencies
diff --git a/gammapy/maps/region/geom.py b/gammapy/maps/region/geom.py
index a152e66a1..52dee0d2a 100644
--- a/gammapy/maps/region/geom.py
+++ b/gammapy/maps/region/geom.py
@@ -4,14 +4,21 @@ import numpy as np
 from astropy import units as u
 from astropy.coordinates import Angle, SkyCoord
 from astropy.io import fits
-from astropy.table import Table
+from astropy.table import QTable, Table
 from astropy.utils import lazyproperty
 from astropy.wcs.utils import (
     proj_plane_pixel_area,
     proj_plane_pixel_scales,
     wcs_to_celestial_frame,
 )
-from regions import CompoundSkyRegion, PixCoord, PointSkyRegion, Regions, SkyRegion
+from regions import (
+    CompoundSkyRegion,
+    PixCoord,
+    PointSkyRegion,
+    Regions,
+    RectanglePixelRegion,
+    SkyRegion
+)
 from gammapy.utils.regions import (
     compound_region_center,
     compound_region_to_regions,
@@ -121,7 +128,12 @@ class RegionGeom(Geom):
         for region_pix in regions_pix[1:]:
             bbox = bbox.union(region_pix.bounding_box)
 
-        rectangle_pix = bbox.to_region()
+        try:
+            rectangle_pix = bbox.to_region()
+        except ValueError:
+            rectangle_pix = RectanglePixelRegion(
+                center=PixCoord(*bbox.center[::-1]), width=1, height=1
+            )
         return rectangle_pix.to_sky(self.wcs)
 
     @property
@@ -583,11 +595,13 @@ class RegionGeom(Geom):
         if self.region is None:
             raise ValueError("Region definition required.")
 
-        # TODO: make this a to_hdulist() method
         region_list = compound_region_to_regions(self.region)
+
         pixel_region_list = []
+
         for reg in region_list:
             pixel_region_list.append(reg.to_pixel(self.wcs))
+
         table = Regions(pixel_region_list).serialize(format="fits")
 
         header = WcsGeom(wcs=self.wcs, npix=self.wcs.array_shape).to_header()
@@ -683,16 +697,22 @@ class RegionGeom(Geom):
             region_hdu = hdu + "_" + region_hdu
 
         if region_hdu in hdulist:
-            region_table = Table.read(hdulist[region_hdu])
-            wcs = WcsGeom.from_header(region_table.meta).wcs
+            try:
+                region_table = QTable.read(hdulist[region_hdu])
+                regions_pix = Regions.parse(data=region_table, format="fits")
+            except TypeError:
+                region_table = Table.read(hdulist[region_hdu])
+                regions_pix = Regions.parse(data=region_table, format="fits")
 
+            wcs = WcsGeom.from_header(region_table.meta).wcs
             regions = []
 
-            for reg in Regions.parse(data=region_table, format="fits"):
+            for region_pix in regions_pix:
                 # TODO: remove workaround once regions issue with fits serialization is sorted out
                 # see https://github.com/astropy/regions/issues/400
-                reg.meta["include"] = True
-                regions.append(reg.to_sky(wcs))
+                region_pix.meta["include"] = True
+                regions.append(region_pix.to_sky(wcs))
+
             region = regions_to_compound_region(regions)
         else:
             region, wcs = None, None
@@ -719,7 +739,7 @@ class RegionGeom(Geom):
             else:
                 self._region = other.region
 
-    def plot_region(self, ax=None, kwargs_point=None, **kwargs):
+    def plot_region(self, ax=None, kwargs_point=None, path_effect=None, **kwargs):
         """Plot region in the sky.
 
         Parameters
@@ -728,6 +748,11 @@ class RegionGeom(Geom):
             Axes to plot on. If no axes are given,
             the region is shown using the minimal
             equivalent WCS geometry.
+        kwargs_point : dict
+            Keyword arguments passed to `~matplotlib.lines.Line2D` for plotting
+            of point sources
+        path_effect : `~matplotlib.patheffects.PathEffect`
+            Path effect applied to artists and lines.
         **kwargs : dict
             Keyword arguments forwarded to `~regions.PixelRegion.as_artist`
 
@@ -751,8 +776,9 @@ class RegionGeom(Geom):
                 ax = m.plot(add_cbar=False)
 
         kwargs.setdefault("fc", "None")
+        kwargs.setdefault("ec", "tab:blue")
         kwargs_point.setdefault("color", kwargs.get("edgecolor"))
-        kwargs_point.setdefault("markeredgecolor", "None")
+        kwargs_point.setdefault("marker", "*")
 
         for region in compound_region_to_regions(self.region):
             region_pix = region.to_pixel(wcs=ax.wcs)
@@ -762,6 +788,9 @@ class RegionGeom(Geom):
             else:
                 artist = region_pix.as_artist(**kwargs)
 
+            if path_effect:
+                artist.add_path_effect(path_effect)
+
             ax.add_artist(artist)
 
         return ax
diff --git a/gammapy/modeling/models/core.py b/gammapy/modeling/models/core.py
index 9526f94f5..96dd05da9 100644
--- a/gammapy/modeling/models/core.py
+++ b/gammapy/modeling/models/core.py
@@ -7,7 +7,6 @@ import numpy as np
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from astropy.table import Table
-from regions import PointSkyRegion
 import yaml
 from gammapy.maps import Map, RegionGeom
 from gammapy.modeling import Covariance, Parameter, Parameters
@@ -952,32 +951,14 @@ class DatasetModels(collections.abc.Sequence):
         ax : `~astropy.visualization.WcsAxes`
             WCS axes
         """
-        from astropy.visualization.wcsaxes import WCSAxes
-
-        kwargs_point = kwargs_point or {}
-
-        if ax is None or not isinstance(ax, WCSAxes):
-            ax = Map.from_geom(self.wcs_geom).plot()
-
-        kwargs.setdefault("color", "tab:blue")
-        kwargs.setdefault("fc", "None")
-        kwargs_point.setdefault("marker", "*")
-        kwargs_point.setdefault("markersize", 10)
-        kwargs_point.setdefault("markeredgecolor", "None")
-        kwargs_point.setdefault("color", kwargs["color"])
-
-        for region in self.to_regions():
-            if isinstance(region, PointSkyRegion):
-                artist = region.to_pixel(ax.wcs).as_artist(**kwargs_point)
-            else:
-                artist = region.to_pixel(ax.wcs).as_artist(**kwargs)
-
-            if path_effect:
-                artist.set_path_effects([path_effect])
-
-            ax.add_artist(artist)
-
-        return ax
+        regions = self.to_regions()
+        geom = RegionGeom.from_regions(regions=regions)
+        return geom.plot_region(
+            ax=ax,
+            kwargs_point=kwargs_point,
+            path_effect=path_effect,
+            **kwargs
+        )
 
     def plot_positions(self, ax=None, **kwargs):
         """ "Plot the centers of the spatial models on a given wcs axis
diff --git a/gammapy/modeling/models/spatial.py b/gammapy/modeling/models/spatial.py
index dced2968b..a46a8e454 100644
--- a/gammapy/modeling/models/spatial.py
+++ b/gammapy/modeling/models/spatial.py
@@ -42,7 +42,6 @@ log = logging.getLogger(__name__)
 MAX_OVERSAMPLING = 200
 
 
-
 def compute_sigma_eff(lon_0, lat_0, lon, lat, phi, major_axis, e):
     """Effective radius, used for the evaluation of elongated models"""
     phi_0 = position_angle(lon_0, lat_0, lon, lat)
@@ -125,12 +124,7 @@ class SpatialModel(ModelBase):
     def position_error(self):
         """Get 95% containment position error as (`~regions.EllipseSkyRegion`)"""
         if self.covariance is None:
-            return EllipseSkyRegion(
-                center=self.position,
-                height=np.nan * u.deg,
-                width=np.nan * u.deg,
-                angle=np.nan * u.deg,
-            )
+            raise ValueError("No position error information available.")
 
         pars = self.parameters
         sub_covar = self.covariance.get_subcovariance(["lon_0", "lat_0"]).data.copy()
@@ -960,14 +954,12 @@ class ConstantSpatialModel(SpatialModel):
         """Evaluate model."""
         return value
 
-    @staticmethod
-    def to_region(**kwargs):
-        """Model outline (`~regions.EllipseSkyRegion`)."""
-        return EllipseSkyRegion(
-            center=SkyCoord(np.nan * u.deg, np.nan * u.deg),
-            height=np.nan * u.deg,
-            width=np.nan * u.deg,
-            angle=np.nan * u.deg,
+    def to_region(self, **kwargs):
+        """Model outline (`~regions.RectangleSkyRegion`)."""
+        return RectangleSkyRegion(
+            center=SkyCoord(0 * u.deg, 0 * u.deg, frame=self.frame),
+            height=180 * u.deg,
+            width=360 * u.deg,
             **kwargs,
         )
 
@@ -1007,14 +999,12 @@ class ConstantFluxSpatialModel(SpatialModel):
         """Evaluate model."""
         return Map.from_geom(geom=geom, data=1)
 
-    @staticmethod
-    def to_region(**kwargs):
-        """Model outline (`~regions.EllipseSkyRegion`)."""
-        return EllipseSkyRegion(
-            center=SkyCoord(np.nan * u.deg, np.nan * u.deg),
-            height=np.nan * u.deg,
-            width=np.nan * u.deg,
-            angle=np.nan * u.deg,
+    def to_region(self, **kwargs):
+        """Model outline (`~regions.RectangleSkyRegion`)."""
+        return RectangleSkyRegion(
+            center=SkyCoord(0 * u.deg, 0 * u.deg, frame=self.frame),
+            height=180 * u.deg,
+            width=360 * u.deg,
             **kwargs,
         )
 
@@ -1083,11 +1073,11 @@ class TemplateSpatialModel(SpatialModel):
             self._map = map.copy(data=map.data)
 
         self.meta = {} if meta is None else meta
-        
+
         interp_kwargs = {} if interp_kwargs is None else interp_kwargs
         interp_kwargs.setdefault("method", "linear")
         interp_kwargs.setdefault("fill_value", 0)
-        
+
         self._interp_kwargs = interp_kwargs
         self.filename = filename
         super().__init__()
diff --git a/gammapy/utils/regions.py b/gammapy/utils/regions.py
index ee9d4b7f1..232dca475 100644
--- a/gammapy/utils/regions.py
+++ b/gammapy/utils/regions.py
@@ -197,7 +197,7 @@ def make_orthogonal_rectangle_sky_regions(start_pos, end_pos, wcs, height, nbin=
     return regions
 
 
-def make_concentric_annulus_sky_regions(center, radius_max, nbin=11):
+def make_concentric_annulus_sky_regions(center, radius_max, radius_min=1e-5 * u.deg, nbin=11):
     """Make a list of concentric annulus regions.
 
     Parameters
@@ -206,6 +206,8 @@ def make_concentric_annulus_sky_regions(center, radius_max, nbin=11):
         Center coordinate
     radius_max : `~astropy.units.Quantity`
         Maximum radius.
+    radius_min : `~astropy.units.Quantity`
+        Minimum radius.
     nbin : int
         Number of boxes along the line
 
@@ -216,7 +218,7 @@ def make_concentric_annulus_sky_regions(center, radius_max, nbin=11):
     """
     regions = []
 
-    edges = np.linspace(0 * u.deg, u.Quantity(radius_max), nbin)
+    edges = np.linspace(radius_min, u.Quantity(radius_max), nbin)
 
     for r_in, r_out in zip(edges[:-1], edges[1:]):
         region = CircleAnnulusSkyRegion(
diff --git a/setup.cfg b/setup.cfg
index 983320bf1..deeaec13c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -33,7 +33,7 @@ install_requires =
     scipy>=1.4
     astropy>=5.0
     scipy>=1.2
-    regions>=0.5.0,<0.6.0a0
+    regions>=0.5.0
     pyyaml>=5.1
     click>=7.0
     pydantic>=1.4


Test patch:
diff --git a/gammapy/makers/tests/test_reduce.py b/gammapy/makers/tests/test_reduce.py
index 5abac9089..169c1d081 100644
--- a/gammapy/makers/tests/test_reduce.py
+++ b/gammapy/makers/tests/test_reduce.py
@@ -75,7 +75,7 @@ def get_spectrumdataset(name):
     )
 
 
-def get_spectrumdataset_rad_max(name, e_min=0.005 * u.TeV):
+def get_spectrum_dataset_rad_max(name, e_min=0.005 * u.TeV):
     """get the spectrum dataset maker for the energy-dependent spectrum extraction"""
     target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
     on_center = PointSkyRegion(target_position)
@@ -154,7 +154,7 @@ def makers_spectrum(exclusion_mask):
     ],
 )
 @requires_data()
-def test_datasetsmaker_map(pars, observations_cta, makers_map):
+def test_datasets_maker_map(pars, observations_cta, makers_map):
     makers = DatasetsMaker(
         makers_map,
         stack_datasets=pars["stack_datasets"],
@@ -185,7 +185,7 @@ def test_datasetsmaker_map(pars, observations_cta, makers_map):
 
 
 @requires_data()
-def test_datasetsmaker_map_cutout_width(observations_cta, makers_map, tmp_path):
+def test_datasets_maker_map_cutout_width(observations_cta, makers_map, tmp_path):
     makers = DatasetsMaker(
         makers_map,
         stack_datasets=True,
@@ -206,7 +206,7 @@ def test_datasetsmaker_map_cutout_width(observations_cta, makers_map, tmp_path):
 
 
 @requires_data()
-def test_datasetsmaker_map_2steps(observations_cta, makers_map, tmp_path):
+def test_datasets_maker_map_2steps(observations_cta, makers_map, tmp_path):
 
     makers = DatasetsMaker(
         [MapDatasetMaker()],
@@ -266,7 +266,7 @@ def test_dataset_maker_spectrum_rad_max(observations_magic_rad_max):
     maker = SpectrumDatasetMaker(
         containment_correction=False, selection=["counts", "exposure", "edisp"]
     )
-    dataset = maker.run(get_spectrumdataset_rad_max("spec"), observation)
+    dataset = maker.run(get_spectrum_dataset_rad_max("spec"), observation)
 
     finder = WobbleRegionsFinder(n_off_regions=1)
     bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=finder)
@@ -294,7 +294,7 @@ def test_dataset_maker_spectrum_global_rad_max():
     maker = SpectrumDatasetMaker(
         containment_correction=False, selection=["counts", "exposure", "edisp"]
     )
-    dataset = maker.run(get_spectrumdataset_rad_max("spec"), observation)
+    dataset = maker.run(get_spectrum_dataset_rad_max("spec"), observation)
 
     finder = WobbleRegionsFinder(n_off_regions=3)
     bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=finder)
@@ -322,7 +322,7 @@ def test_dataset_maker_spectrum_rad_max_overlapping(observations_magic_rad_max,
     bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=finder)
 
     with caplog.at_level(logging.WARNING):
-        dataset = maker.run(get_spectrumdataset_rad_max("spec"), observation)
+        dataset = maker.run(get_spectrum_dataset_rad_max("spec"), observation)
         dataset_on_off = bkg_maker.run(dataset, observation)
 
     assert len(caplog.record_tuples) == 3
@@ -340,7 +340,7 @@ def test_dataset_maker_spectrum_rad_max_overlapping(observations_magic_rad_max,
     # rad max, allowing more off regions
     caplog.clear()
     with caplog.at_level(logging.WARNING):
-        dataset = maker.run(get_spectrumdataset_rad_max("spec", e_min=250 * u.GeV), observation)
+        dataset = maker.run(get_spectrum_dataset_rad_max("spec", e_min=250 * u.GeV), observation)
         dataset_on_off = bkg_maker.run(dataset, observation)
         assert dataset_on_off.counts_off is not None
 
@@ -356,7 +356,7 @@ def test_dataset_maker_spectrum_rad_max_all_excluded(observations_magic_rad_max,
     maker = SpectrumDatasetMaker(
         containment_correction=False, selection=["counts", "exposure", "edisp"]
     )
-    dataset = maker.run(get_spectrumdataset_rad_max("spec"), observation)
+    dataset = maker.run(get_spectrum_dataset_rad_max("spec"), observation)
 
     # excludes all possible off regions
     exclusion_region = CircleSkyRegion(
diff --git a/gammapy/modeling/models/tests/test_spatial.py b/gammapy/modeling/models/tests/test_spatial.py
index 8bcb8d298..714a07015 100644
--- a/gammapy/modeling/models/tests/test_spatial.py
+++ b/gammapy/modeling/models/tests/test_spatial.py
@@ -267,7 +267,7 @@ def test_sky_diffuse_constant():
     assert_allclose(val.value, 42)
     radius = model.evaluation_radius
     assert radius is None
-    assert isinstance(model.to_region(), EllipseSkyRegion)
+    assert isinstance(model.to_region(), RectangleSkyRegion)
 
 
 @requires_dependency("matplotlib")