How to use the regions.CircleSkyRegion function in regions

Returns
    -------
    compound : `~regions.CompoundSkyRegion`
        Compound sky region
    """

    region_union = regions[0]

    for region in regions[1:]:
        region_union = region_union.union(region)

    return region_union


class SphericalCircleSkyRegion(CircleSkyRegion):
    """Spherical circle sky region.

    TODO: is this separate class a good idea?

    - If yes, we could move it to the ``regions`` package?
    - If no, we should implement some other solution.
      Probably the alternative is to add extra methods to
      the ``CircleSkyRegion`` class and have that support
      both planar approximation and spherical case?
      Or we go with the approach to always make a
      TAN WCS and not have true cone select at all?
    """

    def contains(self, skycoord, wcs=None):
        """Defined by spherical distance."""
        separation = self.center.separation(skycoord)

def apply_containment_correction(self, observation, bkg):
        """Apply PSF containment correction.

        Parameters
        ----------
        observation : `~gammapy.data.DataStoreObservation`
            observation
        bkg : `~gammapy.spectrum.BackgroundEstimate`
            background esimate
        """
        if not isinstance(bkg.on_region, CircleSkyRegion):
            raise TypeError(
                "Incorrect region type for containment correction."
                " Should be CircleSkyRegion."
            )

        log.info("Apply containment correction")
        # First need psf
        angles = np.linspace(0.0, 1.5, 150) * u.deg
        offset = observation.pointing_radec.separation(bkg.on_region.center)
        if isinstance(observation.psf, PSF3D):
            psf = observation.psf.to_energy_dependent_table_psf(theta=offset)
        else:
            psf = observation.psf.to_energy_dependent_table_psf(offset, angles)

        new_aeff = apply_containment_fraction(self._aeff, psf, bkg.on_region.radius)

from astropy.coordinates import Angle, SkyCoord
from regions import CircleSkyRegion
import matplotlib.pyplot as plt
from gammapy.makers import ReflectedRegionsFinder
from gammapy.maps import WcsNDMap

# Exclude a rectangular region
exclusion_mask = WcsNDMap.create(npix=(801, 701), binsz=0.01, skydir=(83.6, 23.0))

coords = exclusion_mask.geom.get_coord().skycoord
mask = (Angle("23 deg") < coords.dec) & (coords.dec < Angle("24 deg"))
exclusion_mask.data = np.invert(mask)

pos = SkyCoord(83.633, 22.014, unit="deg")
radius = Angle(0.3, "deg")
on_region = CircleSkyRegion(pos, radius)
center = SkyCoord(83.633, 24, unit="deg")

# One can impose a minimal distance between ON region and first reflected regions
finder = ReflectedRegionsFinder(
    region=on_region,
    center=center,
    exclusion_mask=exclusion_mask,
    min_distance_input="0.2 rad",
)
finder.run()

fig1 = plt.figure(1)
finder.plot(fig=fig1)

# One can impose a minimal distance between two adjacent regions
finder = ReflectedRegionsFinder(

from regions import CircleSkyRegion, make_example_dataset

# load example dataset to get skymap
config = dict(crval=(0, 0),
              crpix=(180, 90),
              cdelt=(-1, 1),
              shape=(180, 360))

dataset = make_example_dataset(data='simulated', config=config)
wcs = dataset.wcs

# remove sources
dataset.image.data = np.zeros_like(dataset.image.data)

# define 2 sky circles
circle1 = CircleSkyRegion(
    center=SkyCoord(20, 0, unit='deg', frame='galactic'),
    radius=Angle('30 deg')
)

circle2 = CircleSkyRegion(
    center=SkyCoord(50, 45, unit='deg', frame='galactic'),
    radius=Angle('30 deg'),
)

# define skycoords
lon = np.arange(-180, 181, 10)
lat = np.arange(-90, 91, 10)
coords = np.array(np.meshgrid(lon, lat)).T.reshape(-1, 2)
skycoords = SkyCoord(coords, unit='deg', frame='galactic')

# get events in AND and XOR

"""Example how to compute and plot reflected regions."""
import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion, EllipseAnnulusSkyRegion, RectangleSkyRegion
import matplotlib.pyplot as plt
from gammapy.maps import WcsNDMap

position = SkyCoord(83.63, 22.01, unit="deg", frame="icrs")

on_circle = CircleSkyRegion(position, 0.3 * u.deg)

on_ellipse_annulus = EllipseAnnulusSkyRegion(
    center=position,
    inner_width=1.5 * u.deg,
    outer_width=2.5 * u.deg,
    inner_height=3 * u.deg,
    outer_height=4 * u.deg,
    angle=130 * u.deg,
)

another_position = SkyCoord(80.3, 22.0, unit="deg")
on_rectangle = RectangleSkyRegion(
    center=another_position, width=2.0 * u.deg, height=4.0 * u.deg, angle=50 * u.deg
)

# Now we plot those regions. We first create an empty map

def from_config(cls, config):
        """Initialize target from config.

        The config dict is stored as attribute for later use by other analysis
        classes.
        """
        obs_id = config['obs']
        if not isinstance(obs_id, list):
            from . import ObservationTable
            obs_table = ObservationTable.read(obs_id)
            obs_id = obs_table['OBS_ID'].data
        # TODO : This should also accept also Galactic coordinates
        pos = SkyCoord(config['ra'], config['dec'], unit='deg')
        on_radius = config['on_size'] * u.deg
        on_region = CircleSkyRegion(pos, on_radius)
        target = cls(pos, on_region, obs_id=obs_id,
                     name=config['name'], tag=config['tag'])
        target.config = config
        return target

cdelt=(-1, 1),
              shape=(180, 360))

dataset = make_example_dataset(data='simulated', config=config)
wcs = dataset.wcs

# remove sources
dataset.image.data = np.zeros_like(dataset.image.data)

# define 2 sky circles
circle1 = CircleSkyRegion(
    center=SkyCoord(20, 0, unit='deg', frame='galactic'),
    radius=Angle('30 deg')
)

circle2 = CircleSkyRegion(
    center=SkyCoord(50, 45, unit='deg', frame='galactic'),
    radius=Angle('30 deg'),
)

# define skycoords
lon = np.arange(-180, 181, 10)
lat = np.arange(-90, 91, 10)
coords = np.array(np.meshgrid(lon, lat)).T.reshape(-1, 2)
skycoords = SkyCoord(coords, unit='deg', frame='galactic')

# get events in AND and XOR
compound_and = circle1 & circle2
compound_xor = circle1 ^ circle2

mask_and = compound_and.contains(skycoords, wcs)
skycoords_and = skycoords[mask_and]

wcs = dataset.wcs

fig = plt.figure()
ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=wcs)

ax.set_xlim(-0.5, dataset.config['shape'][1] - 0.5)
ax.set_ylim(-0.5, dataset.config['shape'][0] - 0.5)

ax.imshow(dataset.image.data, cmap='gray', vmin=0, vmax=1,
          interpolation='nearest', origin='lower')

for source in dataset.source_table:
    # Plot a sky circle around each source
    center = SkyCoord(source['GLON'], source['GLAT'], unit='deg', frame='galactic')
    radius = Angle(20, 'deg')
    region = CircleSkyRegion(center=center, radius=radius)
    pix_region = region.to_pixel(wcs=wcs)

    pix_region.plot(ax=ax, edgecolor='yellow', facecolor='yellow', alpha=0.5, lw=3)

if self.background_model is not None:
            kwargs["background"] = self.background_model.evaluate().get_spectrum(
                on_region, np.sum
            )

        if self.exposure is not None:
            exposure = self.exposure.get_spectrum(on_region, np.mean)
            energy = exposure.geom.axes[0].edges
            kwargs["aeff"] = EffectiveAreaTable(
                energy_lo=energy[:-1],
                energy_hi=energy[1:],
                data=exposure.quantity[:, 0, 0] / kwargs["livetime"],
            )

        if containment_correction:
            if not isinstance(on_region, CircleSkyRegion):
                raise TypeError(
                    "Containement correction is only supported for"
                    " `CircleSkyRegion`."
                )
            elif self.psf is None or isinstance(self.psf, PSFKernel):
                raise ValueError("No PSFMap set. Containement correction impossible")
            else:
                psf = self.psf.get_energy_dependent_table_psf(on_region.center)
                containment = psf.containment(
                    kwargs["aeff"].energy.center, on_region.radius
                )
                kwargs["aeff"].data.data *= containment.squeeze()

        if self.edisp is not None:
            if isinstance(self.edisp, EDispKernel):
                edisp = self.edisp