\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.17?urlpath=lab/tree/exclusion_mask.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[exclusion_mask.ipynb](../_static/notebooks/exclusion_mask.ipynb) |\n",
"[exclusion_mask.py](../_static/notebooks/exclusion_mask.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Creating exclusion masks\n",
"\n",
"## Introduction\n",
"\n",
"### Prerequisites\n",
"\n",
"- Understanding of basic analyses in 1D or 3D.\n",
"- Usage of `~regions` and catalogs, see the [catalog notebook](catalog.ipynb). \n",
"\n",
"### Context\n",
"\n",
"Background templates stored in the DL3 IRF are often not reliable enough to be used without some corrections. A set of common techniques to perform background or normalisation from the data is implemented in gammapy: reflected regions for 1D spectrum analysis, field-of-view (FoV) background or ring background for 2D and 3D analyses.\n",
"\n",
"To avoid contamination of the background estimate from gamma-ray bright regions these methods require to exclude those regions from the data used for the estimation. To do so, we use exclusion masks. They are maps containing boolean values where excluded pixels are stored as False. \n",
"\n",
"**Objective: Build an exclusion mask around the Crab nebula excluding gamma-ray sources in the region.**\n",
"\n",
"### Proposed approach\n",
"\n",
"Here we have to build a `Map` object, we must first define its geometry and then we can determine which pixels to exclude.\n",
"\n",
"We can rely on known sources positions and properties to build a list of regions (here `~regions.SkyRegions`) enclosing most of the signal that our detector would see from these objects. We show below how to build this list manually or from an existing catalog. \n",
"\n",
"Finally, we show how to build the mask from a `MapDataset`, finding pixels which contain statistically significant signal. To do so, we use the `ExcessMapEstimator`\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from regions import read_ds9, CircleSkyRegion\n",
"from gammapy.maps import Map, WcsGeom\n",
"from gammapy.utils.regions import make_region\n",
"from gammapy.catalog import SOURCE_CATALOGS\n",
"from gammapy.datasets import Datasets\n",
"from gammapy.estimators import ExcessMapEstimator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create the mask from a list of regions\n",
"\n",
"One can build an exclusion mask from regions. We show here how to proceed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define the geometry\n",
"\n",
"Exclusions masks are stored in `Map` objects. One has therefore to define the geometry to use. Here we consider a region at the Galactic anticentre around the crab nebula."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"position = SkyCoord(83.633083, 22.0145, unit=\"deg\", frame=\"icrs\")\n",
"geom = WcsGeom.create(\n",
" skydir=position, width=\"5 deg\", binsz=0.02, frame=\"galactic\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create the list of regions\n",
"\n",
"A useful function to create region objects is `~gammapy.utils.regions.make_region`. It can take strings defining regions following the \"ds9\" format and convert them to `regions`. \n",
"\n",
"Here we use a region enclosing the Crab nebula with 0.3 degrees. The actual region size should depend on the expected PSF of the data used. We also add another region with a different shape as en example."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[, width=1.0 deg, height=0.5 deg, angle=45.0 deg)>, , radius=0.3 deg)>]\n"
]
}
],
"source": [
"some_region = make_region(\"galactic;box(185,-4,1.0,0.5, 45)\")\n",
"crab_region = make_region(\"icrs;circle(83.633083, 22.0145, 0.3)\")\n",
"regions = [some_region, crab_region]\n",
"print(regions)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Equivalently the regions can be read from a ds9 file, this time using `regions.read_ds9`. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# regions = read_ds9('ds9.reg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create the mask map \n",
"\n",
"We can now create the map. We use the `WcsGeom.region_mask` method putting all pixels inside the regions to False."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"mask_data = geom.region_mask(regions, inside=False)\n",
"mask_map = Map.from_geom(geom, data=mask_data)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(