{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\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.19?urlpath=lab/tree/tutorials/api/mask_maps.ipynb)\n", "- You may download all the notebooks in the documentation as a\n", "[tar file](../../_downloads/notebooks-0.19.tar).\n", "- **Source files:**\n", "[mask_maps.ipynb](../../_static/notebooks/mask_maps.ipynb) |\n", "[mask_maps.py](../../_static/notebooks/mask_maps.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mask maps\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", "There are two main categories of masks in Gammapy for different use cases.\n", "- Fitting often requires to 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 that the user does not want to model, or both. \n", "Gammapy's `Datasets` therefore contain a `mask_fit` sharing the same geometry as the data (i.e. `counts`).\n", "- During data reduction, some background makers will normalize the background model template\n", "on the data themselves. To limit contamination by real photons, one has to exclude parts of the \n", "field-of-view where signal is expected to be large. To do so, one needs to provide an exclusion mask. \n", "The latter can be provided in a different geometry as it will be reprojected by the `Makers`.\n", "\n", "We explain in more details these two types of masks below:\n", "\n", "### Masks for fitting: `mask_fit`\n", "\n", "The region of interest used for the fit can defined through the dataset `mask_fit` attribute.\n", "The `mask_fit` is a map containing boolean values where pixels used in the fit are stored as True. \n", "\n", "A spectral fit (1D or 3D) can be restricted to a specific energy range where e.g. the background\n", "is well estimated or where the number of counts is large enough.\n", "Similarly, 2D and 3D analyses usually require to work with a wider map than the region of interest so sources laying outside but reconstructed inside because of the PSF are correctly taken into account. Then the `mask_fit` have to include a margin that take into account the PSF width. We will show an example in the boundary mask sub-section.\n", "\n", "The `mask_fit` also can be used to exclude sources or complex regions for which we don't have good enough models. In that case the masking is an extra security, it is prefereable to include the available models even if the sources are masked and frozen.\n", "\n", "Note that a dataset contains also a `mask_safe` attribute that is created and filled during \n", "data reduction. It is not to be modified directly by users. The `mask_safe` is defined only from \n", "the options passed to the `~gammapy.makers.SafeMaskMaker` (More details [here](https://docs.gammapy.org/dev/makers/index.html#safe-data-range-handling)).\n", "\n", "### Exclusion masks\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", "\n", "## Proposed approach\n", "\n", "Even if the use cases for exclusion masks and fit masks are different, the way to create these masks is exactly the same, so in the following we show how to work with masks in general:\n", "- Creating masks from scratch\n", "- Combining multiple masks\n", "- Extending and reducing an existing mask\n", "- Reading and writing masks\n", "\n", "\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.240600Z", "iopub.status.busy": "2021-11-22T21:09:28.235980Z", "iopub.status.idle": "2021-11-22T21:09:28.410758Z", "shell.execute_reply": "2021-11-22T21:09:28.410934Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.413301Z", "iopub.status.busy": "2021-11-22T21:09:28.412930Z", "iopub.status.idle": "2021-11-22T21:09:28.932760Z", "shell.execute_reply": "2021-11-22T21:09:28.932955Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from astropy.coordinates import SkyCoord, Angle\n", "import astropy.units as u\n", "from regions import CircleSkyRegion, Regions\n", "from gammapy.maps import Map, WcsGeom, MapAxis\n", "from gammapy.catalog import CATALOG_REGISTRY\n", "from gammapy.datasets import Datasets\n", "from gammapy.estimators import ExcessMapEstimator\n", "from gammapy.modeling.models import FoVBackgroundModel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a mask for fitting\n", "\n", "One can build a `mask_fit` to restrict the energy range of pixels used to\n", "fit a `Dataset`. The mask being a `Map` it needs to use the same geometry\n", "(i.e. a `Geom` object) as the `Dataset` it will be applied to. \n", "\n", "We show here how to proceed on a `MapDataset` taken from Fermi data used in the 3FHL catalog. The dataset is already in the form of a `Datasets` object. We read it from disk.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.935100Z", "iopub.status.busy": "2021-11-22T21:09:28.934813Z", "iopub.status.idle": "2021-11-22T21:09:28.946171Z", "shell.execute_reply": "2021-11-22T21:09:28.946358Z" } }, "outputs": [], "source": [ "filename = \"$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml\"\n", "datasets = Datasets.read(filename=filename)\n", "dataset = datasets[\"Fermi-LAT\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check the default energy range of the dataset. In the absence of a `mask_fit` it is equal to the safe energy range." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.948216Z", "iopub.status.busy": "2021-11-22T21:09:28.947909Z", "iopub.status.idle": "2021-11-22T21:09:28.989634Z", "shell.execute_reply": "2021-11-22T21:09:28.989871Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "HDU 'MASK_FIT' not found\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Fit range : (WcsNDMap\n", "\n", "\tgeom : WcsGeom \n", " \taxes : ['lon', 'lat']\n", "\tshape : (50, 40)\n", "\tndim : 2\n", "\tunit : TeV\n", "\tdtype : float64\n", ", WcsNDMap\n", "\n", "\tgeom : WcsGeom \n", " \taxes : ['lon', 'lat']\n", "\tshape : (50, 40)\n", "\tndim : 2\n", "\tunit : TeV\n", "\tdtype : float64\n", ")\n" ] } ], "source": [ "print(f\"Fit range : {dataset.energy_range}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create a mask in energy\n", "\n", "We show first how to use a simple helper function `~gammapy.maps.Geom.energy_range()`.\n", "\n", "We obtain the `Geom` that is stored on the `counts` map inside the `Dataset`\n", "and we can directly create the `Map`.\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.992000Z", "iopub.status.busy": "2021-11-22T21:09:28.991696Z", "iopub.status.idle": "2021-11-22T21:09:28.992913Z", "shell.execute_reply": "2021-11-22T21:09:28.993092Z" } }, "outputs": [], "source": [ "mask_energy = dataset.counts.geom.energy_mask(10 * u.GeV, 700 * u.GeV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now set the dataset `mask_fit` attribute. \n", "\n", "And we check that the total fit range has changed accordingly. The bin edges closest to requested range provide the actual fit range. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:28.995454Z", "iopub.status.busy": "2021-11-22T21:09:28.995129Z", "iopub.status.idle": "2021-11-22T21:09:28.996787Z", "shell.execute_reply": "2021-11-22T21:09:28.997060Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fit range : (WcsNDMap\n", "\n", "\tgeom : WcsGeom \n", " \taxes : ['lon', 'lat']\n", "\tshape : (50, 40)\n", "\tndim : 2\n", "\tunit : TeV\n", "\tdtype : float64\n", ", WcsNDMap\n", "\n", "\tgeom : WcsGeom \n", " \taxes : ['lon', 'lat']\n", "\tshape : (50, 40)\n", "\tndim : 2\n", "\tunit : TeV\n", "\tdtype : float64\n", ")\n" ] } ], "source": [ "dataset.mask_fit = mask_energy\n", "print(f\"Fit range : {dataset.energy_range}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mask some sky regions\n", "\n", "One might also exclude some specific part of the sky for the fit. For instance, if one wants not to model a specific source in the region of interest, or if one want to reduce the region of interest in the dataset `Geom`.\n", "\n", "In the following we restrict the fit region to a square around the Crab nebula. **Note**: the dataset geometry is aligned on the galactic frame, we use the same frame to define the box to ensure a correct alignment. \n", "We can now create the map. We use the `WcsGeom.region_mask` method putting all pixels outside the regions to False (because we only want to consider pixels inside the region. For convenience we can directly pass a ds9 region string to the method:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.007283Z", "iopub.status.busy": "2021-11-22T21:09:29.006912Z", "iopub.status.idle": "2021-11-22T21:09:29.011255Z", "shell.execute_reply": "2021-11-22T21:09:29.011483Z" } }, "outputs": [], "source": [ "regions = \"galactic;box(184.55, -5.78, 3.0, 3.0)\"\n", "mask_map = dataset.counts.geom.region_mask(regions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now combine this mask with the energy mask using the logical and operator" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.013566Z", "iopub.status.busy": "2021-11-22T21:09:29.013279Z", "iopub.status.idle": "2021-11-22T21:09:29.014635Z", "shell.execute_reply": "2021-11-22T21:09:29.014807Z" } }, "outputs": [], "source": [ "dataset.mask_fit &= mask_map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check the result and plot the full mask." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.017125Z", "iopub.status.busy": "2021-11-22T21:09:29.016836Z", "iopub.status.idle": "2021-11-22T21:09:29.289022Z", "shell.execute_reply": "2021-11-22T21:09:29.289267Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.\n", " return super().imshow(X, *args, origin=origin, **kwargs)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAACyCAYAAACA5hTOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbnklEQVR4nO3deZRtZX3m8e/DReEyqoBxilxMjBDTgiLOKKJGMXEFGhMvjizHuDDd6nJOOssYE1uUaEewTWfQG6cYBY2aOMQBRQUVCQgkcWhxtqNBGe4FB+DXf5xdsG9Z59Q5p2qfU+fu72etvWoP797v+9Z+6t56aw8nVYUkSZIk9c1u826AJEmSJM2DgyFJkiRJveRgSJIkSVIvORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSLy3UYCjJ15Ncm2R7azp93u1aTZItST6e5Jok/57koSPKJskrk1zeTKcmyQrljm59D3YkqWXflzuusM99m7L7rrDtX5I8y/4NbbfZu6nchjk3E/bvwU3ZK5N8fUS5DdO/Zv8+ZG/hzs2E/Xt+kkuSXJ3ksiTP3+j9a/bvQ/YW7txM2L8PLGvnT5NcvJH71+zfh+wt1LlJcuskb0/y3Qz+rf50knuPKJ8syu8UVbUwE/B14KEzqGf3dT7eucCfAZuBE4ErgIOGlH0G8CXgDsDtgX8FfneV428Bapx2N8c+edm6XwN+Ahxg/8zeIp2bCft3L+AJwNOBr495fLM3m+wt3LmZsH8vAO4B7A7cBfgGsHUj969H2Vu4czNJ/1bY92zgDzdy//qSvUU7N8CdgOcCtwU2Mfj3+j+BfYaUX5jfKToN2XpPo344gJOBTwGvBn4EXAYc19q+P/DXwPeA7wAvBza19v008Brgh822A4D3AVcBn2/WfaopfwZw2rL63wc8e4V2/UpzYvZtrTtnWCCAzwBPby0/BThvkvCs0teXAB9btv+pwFlTnpNdun9mb+Oem0n71yrzUKb8hdvsdZO9RTs30/avVfbPgddt1P71MXuLcm7W0r+m3dcDh2zU/vU1e4t0bpYd7yrgyCHbFuZ3ioW6TW4M92YwUjyQwTfkr1uX5LYB1wG/DNwd+HXgqcv2/Rpwa+BPGPwA7ABuAzypmWgd66QkuwEkORB4CPD2Fdp0V+BrVXV1a91FzfqV3LXZPk7ZYUb19c3A0UuXHZs+PBb42wnraLd3V+7fuMzeTe3bqNlbD2avm+yth4XIXnNujgYunbBOs9dx9hbk3KzlZ+uJwDlVddmEdZq97v/dW7hzk+QI4ObAV4cUWZzfKdZrdDiLicFfCrYzuOy4ND2tbhrtf7VVdi8Go8vbAL/AYLS+ubX9JODjrX2/2dq2CfgZcJfWuhv/UtAs/xvwsGb+WcA/DWnzE1g2Embww/emIeWvBw5tLd+56UfGGUmv1tdm+SPAS5r5hzG4zHmzKc/JLt0/s7dxz82k/WuVmerqg9nrLnuLdm6m7V9T7o8Y/FKwx0btXx+ztyjnZo39+yrLbivaaP3rcfYW5tw0x9oPuBh48YgyC/M7xSJeGTq+qm7Rmv6yte3/Lc1U1TXN7D7AwcDNgO8luSLJFcBfMPirwJJvteYPYnAivjVkOwxGq49v5h/PYIS6ku0MQtO2H3D1CmVXKr8fsL2aMz2Gcfq6jcFfIWDww/u2qvrZSgdb9oDf48Zo71KbF6J/EzJ7o2307K2V2esue2u1ENlrHup9IvAbVfWTVfrUZvY6zt5GOTcdZu8BDAYL71qlP8uZve6ztyHOzRjZWyq3mcGtgudV1StGtG9D/U4xyu6T7rCgvsVgdHlgVV03pEz75PyAwWW5OwBfbtb94rLybwEuSXI4cBjwniHHvRS4U5J966ZLp4cDbxtR/nDgc62yk1yyH6evZwGvT/Jg4L8Cxww7WFUdt0p9C92/GTB7O5tn9tbK7HWXvbXa8NlL8mTgRcADq+rbqxx/ObPXYfY20rnp8N+9JzF4lmL7KuWWM3vd/7u3Ic7NGNkjyR4M+v8dBi9IGGVD/U4x0lovlc1yYowH6patK+CXm/l/AP4Xg5HpbsAvAQ8ase87GAR4L+BQ4JsrlPln4IvA36zS7vMYPOi3J3ACo99s87sMLsneHrgdg+BM9PaNUX1t7fPG5vt56Tqcl126f2Zv456bCfu3W1PuOAZvjNoTuPlG7l+Psrdw52bC/j2OwV+yD5vg+GZvNtlbuHMzSf+a8pubMscuQv/6kr1FOzcMrtK8j8FgaJy3vS3M7xRrDuwsp6az1zK49LY0vXtEwNs/HPsD/xv4NnAl8C80r88csu9BwD9y09tFXgl8dFmZxzd1PHiMk3t20/Yv0foBZ/Cw5vbWchg8DPjDZjqVEfdXDgnP0L629jmm2eeF63Bedun+mb2Ne24m7N9Sve3p7I3cvx5lb+HOzYT9u4zBcwntc/iGjdy/HmVv4c7NJP1r1p3E4I8MI/893yj960v2Fu3cAA9qjnPNsvNy9JCfrYX5nSLNgbSKJK8EblNVT2qteyCDy6dbquqGuTVOuzSzp3kxe5oXs6d5MXv9s4gvUJiJJIcmuVsG7sXg/ejvbm2/GfDfgb/yB0PryexpXsye5sXsaV7MnhwMDbcvgwezdgB/D5zG4N5FkhzG4B7P2wKvnU/ztAsze5oXs6d5MXuaF7PXc94mJ0mSJKmXvDIkSZIkqZccDEmSJEnqpYX70NUk3tenNbnVrW5Vl19++cR/CDB7WqtpswfmT2tj9jQvZk/zNE7+vDKk3jnkkEMy7zaon8ye5sXsaV7MnuZpnPzNfDCU5Kgk1yd5dGvd1iQXJHn2rNsjSZIkqZ9mOhhKsonBJ/t+aNmmrcBRwH2S7DPLNkmSJEnqp1lfGfo94Ezg+8vWL13Cqta8JEmSJHVmZoOhJLcHTgDesMLms4DzgfOr6upZtUmSJElSf83ybXKvBV5YVdcnO1/8qaptwLYZtkWSJElSz3U6GEpyCvC0ZnF/4O+agdCBwCOTXFdV7xnjODs6a6Q0gtnTPJk/zYvZ07yYPc1ap4OhqjoDOGP5+iRvAt4/zkCoOc7erX1957xmxuxpnsyf5sXsaV7MnmbNzxmSJEmS1EuzfGboRlV18jzqlSRJkqQlXhmSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSL83lBQqSxlPlW0W7svzDn/XzzF93zN9oZq87Zk/amVeGJEmSJPXSWIOhJAcneWgzvznJvt02S5IkSZK6tepgKMnTgHcBf9GsugPwng7bJEmSJEmdG+fK0CnA/YGrAKrqK8Ctu2yUJEmSJHVtnMHQT6rqp0sLSXYHfLJRkiRJ0kIbZzD0iSQvATYneRjwTuB93TZLkiRJkro1zmDoRcAPgIuBZwD/BPxBl42SJEmSpK6t+jlDVXUD8JfNJEmSJEm7hKGDoSQXM+LZoKq6WyctkiRJkqQZGHVl6Debr6c0X9/cfH0ccE1nLZIkSZKkGRg6GKqqbwAkuX9V3b+16UVJPg28rOvGSZIkSVJXxnmBwt5JHrC0kOR+wN7dNUmSJEmSurfqCxSApwB/k2T/ZvkK4MmdtUiSJEmSZmCct8l9ATg8yX5AqurK7pslSZIkSd1adTCU5A+XLQNQVT4zJEmSJGlhjfPM0I7WdD1wHLBl0oqSHJPkyiQXNtMftrZtTXJBkmdPelxJkiRJmsY4t8md1l5O8mrgvVPWd05V/eYK67cCRwFvTbJPVW2f8viSJEmSNJZxrgwttxdwp3VuR5qv1ZqXJEmSpM6M88zQxQwGKQCbgIOAP56yvvsmuQj4LvC8qrq0WX8WcD7wlqq6espjS5IkSdLYxnm1dvu2tuuA/6iq66ao6wLg4KranuSRwHuAOwNU1TZg27Adk+yYoj5pzcye5sn8aV7MnubF7GnWxrlN7uVV9Y1m+k5VXZfkzeMcPMkpSy9MAG58Fqiq/gm4WZIDxzlOVe29NI1TXlovZk/zZP40L2ZP82L2NGvjDIbu2l5Isjtw5DgHr6ozquqIqjoCuCHNe7mT3Kup+/LJmitJkiRJ62PobXJJXgy8BNic5Kql1cBPgf8zRV2PBp6Z5DrgWmBrVdUq+0iSJElSJ4YOhqrqFcArkryiql681oqq6nTg9LUeR5IkSZLWw6grQ4dW1b8D70xyj+Xbq+qCTlsmSZIkSR0a9Ta55wJPB05bYVsBx3bSIkmSJEmagVG3yT29mT2uqn7c3pZkz05bJUmSJEkdG+dtcp8Zc50kSZIkLYxRzwzdBrg9g7fJ3Z3Bm+QA9gP2mkHbJEmSJKkzo54ZejhwMnAH4M9a669m8MptSZIkSVpYo54Z2gZsS3JiVZ05wzZJkiRJUudGXRkCoKrOTPIbwF2BPVvrX9ZlwyRJkiSpS6u+QCHJG4DHAL/H4Lmh3wYO7rhdkiRJktSpcd4md7+qeiLwo6r6I+C+wC922yxJkiRJ6tY4g6Frm6/XJLkd8DPgkO6aJEmSJEndW/WZIeD9SW4BvAq4ACjgL7tslCRJkiR1bZwXKPxxM3tmkvczeInCoZ22SpIkSZI6Ns5tcjeqqp9U1ZXAOztqjyRJkiTNxESDoZasayskSZIkacamHQzVurZCkiRJkmZs6DNDSd7HyoOeAAd01iJJkiRJmoFRL1B49ZTbJEmSJGnDGzoYqqpPzLIhkiRJkjRL0z4zJEmSJEkLbaaDoSTHJLkwyaVJPtFavzXJBUmePcv2SJIkSeqvVT90db0kuQXweuARVfXNJLdubd4KHAW8Nck+VbV9Vu2SJEmS1E+rXhlK8s/NQGZp+ZZJPjRFXY8FzqqqbwJU1ffb1TRfCz/DSJIkSdIMjHOb3IFVdcXSQlX9CLj18OJD/QpwyyRnJ/lCkie2tp0FnA+cX1VXT3FsSZIkSZrIOLfJ3ZDkjktXdJIczHQfuro7cCTwEGAzcG6S86rqy1W1Ddg2bMckO6aoT1ozs6d5Mn+aF7OneTF7mrVxBkO/D3yq9cKDBwJPH+fgSU4BntYs/j3wwaraAexI8kngcODLqx2nqvZuHXOagZg0FbOneTJ/mhezp3kxe5q1VQdDVfXBJPcA7sPgeZ7nVNV/jnPwqjoDOAMgyWHA6Ul2B24O3Bt4zbQNlyRJkqS1GDoYSnJoVf17MxAC+G7z9Y7NbXMXTFJRVf1bkg8CXwRuAP6qqi6ZqtWSJEmStEajrgw9l8HtcKetsK2AYyetrKpeBbxq0v0kSZIkab0NHQxV1dJzQcdV1Y/b25Ls2WmrJEmSJKlj47xa+zNjrpMkSZKkhTHqmaHbALcHNie5Ozd9GOp+wF4zaJskSZIkdWbUM0MPB04G7sDguaGlwdBVwEu6bZYkSZIkdWvUM0PbgG1JTqyqM2fYJkmSJEnq3DjPDB2Z5BZLC0lumeTl3TVJkiRJkro3zmDouKq6Ymmhqn4EPLKzFkmSJEnSDIwzGNqUZI+lhSSbgT1GlJckSZKkDW/UCxSWvAX4aJI3Mviw1ScD2zptlSRJkiR1bNXBUFWdmuRi4CEM3ij3x1X1oc5bJkmSJEkdGufKEFX1AeADHbdFkiRJkmZm1WeGktwnyeeTbE/y0yTXJ7lqFo2TJEmSpK6M8wKF04GTgK8Am4GnAq/rslGSJEmS1LVxb5P7apJNVXU98MYkn+m4XZIkSZLUqXEGQ9ckuTlwYZJTge8Be3fbLEmSJEnq1ji3yT0B2AQ8C9gB/CJwYpeNkiRJkqSujfNq7W80s9cCf9RtcyRJkiRpNoYOhprPFqph26vqbp20SJIkSZJmYNSVod+cWSskSZIkacaGDoZat8etiyTPBx7Xqvcw4KCq+mGSrcALgL+tqteuZ72SJEmStJKZfehqVb2qqo6oqiOAFwOfqKofNpu3AkcB90myz6THliRJkqRJzetDV08C3t5aTvO1WvOSJEmS1JlxBkNU1VeBTVV1fVW9EXjwtBUm2Qt4BHBma/VZwPnA+VV19bTHliRJkqRxzeNDVx8FfLp1ixxVtQ3YtoZjSpIkSdJExv3Q1d2Y4kNXk5yS5MJmul2zeis73yI3znF2LE2T7CetldnTPJk/zYvZ07yYPc3aJB+6+mMm/NDVqjoDOGNpOcn+wIOAx094nBuvRCUZ+tlH0noze5on86d5MXuaF7OnWRt6ZSjJbyU5pbX82SRfa6ZHT1nfCcCHq8rRviRJkqS5GnWb3AuA97aW92Dw+utjgGdOU1lVvamqtk6zryRJkiStp1G3yd28qr7VWv5UVV0OXJ5kLS9QkCRJkqS5G3Vl6Jbthap6VmvxoG6aI0mSJEmzMWow9NkkT1u+MskzgM911yRJkiRJ6t6o2+SeA7wnyWOBC5p1RzJ4duj4jtslSZIkSZ0aOhiqqu8D90tyLHDXZvU/VtXHZtIySZIkSerQOJ8z9DHAAZA0B0nm3QT1mPnTvJg9SbMy6pkhSZIkSdplORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSLzkYkiRJktRLDoYkSZIk9ZKDIUmSJEm95GBIkiRJUi85GJIkSZLUSw6GJEmSJPWSgyFJkiRJveRgSJIkSVIv7T7vBkzqgAMOYMuWLVx22WUccsghM63bOneZ+mqa/Zay1xxjl/4e9aXORckezO/fvj7kYB51mr3x7OrnpQ91mr2NWV9f6hwnf6maOqNzlWRHVe1tnYtf5zz6uFZ9+B71oU6zt/Hq60udZm9j1tmHPs6rzrXow/eoD32cV52r8TY5SZIkSb3kYEiSJElSLy3sbXKSJEmStBZeGZIkSZLUSw6GJEmSJPXShhoMJfmbJN9Pcklr3RFJzktyYZLzk9yrWf+4Zt3SdEOSI5ptxzRlT13POpttd0tybpJLk1ycZM8u60yyJcm1rX6+obVPZ/1stt8xyfYkz5u0zgn7eK9W/y5KcsI0fVyrWefP7Jm9Kdts9szeupl19iats9m2pvz1IXtT9HPu+TN7Zm9e2fs5VbVhJuCBwD2AS1rrPgwc18w/Ejh7hf3+C/C11vI7gM3AacCh61Ung89l+iJweLN8ALCp4zq3tMstO04ndba2nwm8E3jepHVO2Me9gN2b+dsC328tj93HRcuf2TN7Zs/s9S1788hfH7K3iPkze2ZvXtlbPm2oK0NV9Ungh8tXA/s18/sD311h15OAt7eWd2v2uwHIOtb568AXq+qiZt/Lq+r6juscpbM6kxwPfA24dJo6J6mvqq6pquua9Xs25Saqbz3MOn9mz+zd2Dizt7zOUczeOurD/7t9yN6kdW6E/Jm9n6tvFLPXpVmNusadWDYyBg4Dvgl8C/gOcPAK+/xf4Ndayw8HLgBOW886gWcDbwY+1Bz/BTOocwuwA/gX4BPA0TOoc2/gXGAf4KXs/JeCseuc5FwC92bwg7gdOGHaPi5a/sye2TN7Zq9v2ZtH/vqQvUXMn9kze/PK3k5tn2VlU34z/xw4sZn/HeAjy8rfG7h4FnUCzwMuAw5kcJnvXOAhHde5B3BAM39kE6r9Oq7z1cDvNPM7/XB0eS6b9YcBnwP27EP+zJ7ZM3tmr2/Zm0f++pC9Rcyf2TN788reTu2YZ+VjfjOvhBs/DynAVcvKvwZ4ySzqBLYCb2qV+x/A82fRz1a5s4F7dtzPc4CvN9MVDC59PmuGffz4tH1ctPyZPbNn9sxe37I3j/z1IXuLmD+zZ/bmlb32tKGeGRriu8CDmvljga8sbUiyG/DbwN/NqM4PAXdLsleS3Zsy/9plnUkOSrKpmb8TcGcG93Z2VmdVHV1VW6pqC/Ba4E+r6vSu6ktySPP9JMnBwF0Y/GBuBLPOn9kze0vMHmZvTvrw/24fsje0zg2cP7OH2Zu5eY7EVhgdvh34HvAz4NvAU4AHAF8ALgI+CxzZKn8McN6M63w8g/scLwFO7bpO4MSmvosY3Ef5qFn0s7XfS5nisumEfXxC08cLmz4e34f8mT2zZ/bMXt+yN4/89SF7i5g/s2f25pW95dPS5StJkiRJ6pVFuE1OkiRJktadgyFJkiRJveRgSJIkSVIvORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYGiEJL+Q5G1JvpbkC0nOTXLCKvtsSXLJlPWdnOR2reW/SvKrY+57TJL3T1PvmMe/XZJ3NfNHJHnkFMd4aZLnrX/rdj1mb6fjm70ZMns7Hd/szZDZ2+n4Zm+GzN5Ox+9d9hwMDZEkwHuAT1bVnarqSGArcIcOqz0ZuPGHo6qeWlXr9Unva1JV362qRzeLRwAT/3BoPGZvZ2Zvdszezsze7Ji9nZm92TF7O+tj9hwMDXcs8NOqesPSiqr6RlW9Dm78i8A5SS5opvstP8CoMklekOTiJBcl+Z9JHg3cE3hrkguTbE5ydpJ7NuUf0RzjoiQfHbcTSU5q6rkkyStb67cn+ZPmeOcl+YVm/S81y59P8rIk21t9uSTJzYGXAY9p2vmY5X8BaMptaeZ/P8mXknwEuEurzC8l+WDzF5hzkhw6bp96wOyZvXkxe2ZvXsye2ZsXs9f37FWV0woT8N+A14zYvhewZzN/Z+D8Zn4LcMkqZY4DPgPs1Szfqvl6NnDPVh1nM/iBOQj4FnBIu/yy9hwDvH/ZutsB32z23x34GHB8s62ARzXzpwJ/0My/Hzipmf9dYPsK/ToZOL1Vz0uB57WWL2nKHwlc3Hwf9gO+ulQO+Chw52b+3sDH5n3ON8pk9sye2TN7fZvMntkze2ZvXtPuaCxJzgAewOCvB0cBNwNOT3IEcD3wKyvsNqzMQ4E3VtU1AFX1w1Wqvw+Dy7eXjVl+yVHA2VX1g6YPbwUeyOBy8E8Z/CAAfAF4WDN/X+D4Zv5twKvHrGslRwPvXupnkvc2X/cB7ge8M8lS2T3WUM8uzexNxeytA7M3FbO3DszeVMzeOjB7U1no7DkYGu5S4MSlhao6JcmBwPnNqucA/wEczuB2wx+vcIxhZcJgpD6uScu39xvmZ9UM0xn84K4lC9ex8y2Xe7bmV2r3bsAVVXXEGurclZm98Zm99WX2xmf21pfZG5/ZW19mb3y7ZPZ8Zmi4jwF7Jnlma91erfn9ge9V1Q3AE4BNKxxjWJkPA09OshdAkls1668G9l3hOOcCD0pyyLLyq/lss9+BSTYBJwGfWGWf87jpH4WtQ8osb+fXgXs0bbsHcEiz/pPACc39sPsCjwKoqquAy5L8drNPkhw+Zp/6wOyZvXkxe2ZvXsye2ZsXs9fz7DkYGqIZRR/PIFyXJfkcsA14YVPk9cCTkpzH4HLojhUOs2KZqvog8F7g/CQXAksPo70JeEOaB+pabfkB8HTgrCQXAe8Y0uyHJPn20sTgPs4XAx8HLgIuqKp/WKXrzwae2/T3tsCVK5T5OPCrTTsfA5wJ3KrpyzOBLzftvqBp64VNmXNax3gc8JSmP5cCv7VKu3rD7Jm9eTF7Zm9ezJ7ZmxezZ/Zy05UzCZq/XlxbVZVkK4OH6zZccLXrMXuaF7OneTF7mhezdxOfGdJyRzJ4CDDAFcCT59sc9YjZ07yYPc2L2dO8mL2GV4YkSZIk9ZLPDEmSJEnqJQdDkiRJknrJwZAkSZKkXnIwJEmSJKmXHAxJkiRJ6iUHQ5IkSZJ66f8DLHV0HnNRUJAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "_ = dataset.mask_fit.plot_grid(ncols=5, vmin=0, vmax=1, figsize=(14, 3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating a mask manually\n", "\n", "If you are more familiar with the `Geom` and `Map` API, you can also create the mask manually from the coordinates of all pixels in the geometry. Here we simply show how to obtain the same behaviour as the `energy_mask` helper method.\n", "\n", "In practice, this allows to create complex energy dependent masks." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.291691Z", "iopub.status.busy": "2021-11-22T21:09:29.291390Z", "iopub.status.idle": "2021-11-22T21:09:29.293741Z", "shell.execute_reply": "2021-11-22T21:09:29.293967Z" } }, "outputs": [], "source": [ "coords = dataset.counts.geom.get_coord()\n", "mask_data = (coords[\"energy\"] >= 10 * u.GeV) & (coords[\"energy\"] < 700 * u.GeV)\n", "mask_energy = Map.from_geom(dataset.counts.geom, data=mask_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating an exclusion mask\n", "\n", "Exclusion masks are typically used for background estimation to mask out regions where gamma-ray signal is expected. \n", "An exclusion mask is usually a simple 2D boolean `Map` where excluded positions are stored as `False`. Their actual geometries are independent of the target datasets that a user might want to build. The first thing to do is to build the geometry.\n", "\n", "### Define the geometry\n", "\n", "Masks are stored in `Map` objects. We must first define its geometry and then we can determine which pixels to exclude. Here we consider a region at the Galactic anticentre around the crab nebula." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.297841Z", "iopub.status.busy": "2021-11-22T21:09:29.297509Z", "iopub.status.idle": "2021-11-22T21:09:29.298649Z", "shell.execute_reply": "2021-11-22T21:09:29.298895Z" } }, "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 mask from a list of regions\n", "\n", "One can build an exclusion mask from regions. We show here how to proceed.\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.\n", "\n", "A useful function to create region objects is `~regions.regions.parse`. It can take strings defining regions e.g. 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": 12, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.346246Z", "iopub.status.busy": "2021-11-22T21:09:29.345880Z", "iopub.status.idle": "2021-11-22T21:09:29.347126Z", "shell.execute_reply": "2021-11-22T21:09:29.347359Z" } }, "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": [ "regions_ds9 = (\n", " \"galactic;box(185,-4,1.0,0.5, 45);icrs;circle(83.633083, 22.0145, 0.3)\"\n", ")\n", "regions = Regions.parse(regions_ds9, format=\"ds9\")\n", "print(regions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equivalently the regions can be read from a ds9 file, this time using `Regions.read`. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.349070Z", "iopub.status.busy": "2021-11-22T21:09:29.348787Z", "iopub.status.idle": "2021-11-22T21:09:29.350104Z", "shell.execute_reply": "2021-11-22T21:09:29.350335Z" } }, "outputs": [], "source": [ "# regions = Regions.read('ds9.reg', format=\"ds9\")" ] }, { "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": 14, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.352131Z", "iopub.status.busy": "2021-11-22T21:09:29.351850Z", "iopub.status.idle": "2021-11-22T21:09:29.362196Z", "shell.execute_reply": "2021-11-22T21:09:29.362390Z" } }, "outputs": [], "source": [ "# to define the exclusion mask we take the inverse\n", "mask_map = ~geom.region_mask(regions)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.366786Z", "iopub.status.busy": "2021-11-22T21:09:29.366497Z", "iopub.status.idle": "2021-11-22T21:09:29.420168Z", "shell.execute_reply": "2021-11-22T21:09:29.420391Z" }, "nbsphinx-thumbnail": { "tooltip": "Create and apply masks maps." } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEHCAYAAABbS3hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUfElEQVR4nO3deZgkdX3H8feHJbDscsSZxQMPZvTRoCawYWcFQXQDXqA8wgPq4BWCgiFrEvTBAzV58MghSDRx14ckPgqPoibIikgUPJDDyOGw2XUXNWoY8cBHtFdkDxRhv/mjqqG26en59cxUdXX35/U89Uzd9e3Zrs9W1VT9ShGBmVmK3XpdgJn1DweGmSXr68CQtL3XNaTqp1qhv+rtp1qh/+ot6uvAMLNqOTDMLJn67a8ky5Yti7GxMQCmp6cZHx/vbUGJ+qlW6K96+6lW6I96p6eno9FoPOyAYvdeFDMfY2NjTE1N9boMs4E2MTGhduN9SmJmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwY1jVJSG3/TG8DzoFhyVqDwqExfBwYlmSmcHBoDBcHhs1qtlBwaAwPB4Z1lBoGDo3h4MCwBePQGHwODJvRXALAoTHYHBjW1nx2fIfG4HJg2C4W6h4Lh8ZgcmDYgxZ6J3doDB4HhgHl7dwOjcHiwLDSd2qHxuBwYFglHBqDwYEx5KrckR0a/c+BMcS8A1u3+u41AzZ/DgqbKx9hDJlehkW/vTTLHs6BMUQcFjZfDowh4bCwhVB5YEhaKekBSScXxk1KWi/prKrrsXI5LAZLpYEhaRHwPuDqlkmTwErgcEl7V1nTMOjV0YXDYvBUfYTxl8BlwF0t45vf6Cj02wJwWNhCqiwwJD0WOBG4sM3kdcAUMBURW6uqaZD1smVvh8XgqvI+jA8Cb42IB1q/yBFxMXDxTAtK2t7sHxkZKau+geELnDZf09PTu+x3EbEUSg4MSauB0/PB/YBP51/mZcBxku6PiMtnW0+zWICJiQl/IztwWNhCGB8fp9FoLG0dX2pgRMRaYG3reEkXAVemhIX1B4fFcPB9GAPE1yysbD15liQiTu3FdgeZw8Kq4IfP+pyvWViVfErSxxwWVjUHhnXNYTG8HBh9rBc7rsNiuDkw+lyVO7DDwhwYA6CKHdlhYeDAGBhl7tAOC2tyYAyQMnZsh4UVOTAGzELu4A4La+XAGEALsaM7LKwdB8aAms8O77CwmTgwBthcdnyHhXXiwBhw3QSAw8Jm48AYAilB4LCwFA6MIdEpEBwWlsqBMUTaBYPDwrrhwBgyxYBwWFi33IDOEHJQ2Fz5CMPMkjkwzCyZA8PMkjkwzCyZA8PMkjkwzCyZA8PMkjkwzCyZA8PMkjkwzCyZA8PMkjkwzCyZA8PMkiUFhqQDJT03799L0j7llmVmdTRrYEg6HfgM8K/5qMcBl5dYk5nVVMoRxmrgSOAegIj4PvDIMosys3pKCYzfRsR9zQFJuwNdt8AiaZWkX0vakHd/W5g2KWm9pLO6Xa+ZVSelxa3rJL0d2EvS84C/AD4/x+3dEBEvbjN+ElgJXCJp74jYNsf1m1mJUo4w3gb8AtgEvB74AvDOBa5D+c8o9JtZzcx6hBERO4F/z7v5eqakjcCdwNkRcVs+fh0wBXwiIrYuwHbMrASaqUFYSZvocK0iIg7uakPSvsDOiNgm6TjgnyPiyYnLbm/2j4yMLGk0Gt1s2sy6NDo6ypYtW3Y0hyNiKXQ+wmhea1id//x4/vOVwI6Hz/5wklYDp+eDx0XEnfnGvyDpw5KWRcQvZ1tPs1iAiYkJN3ltVrLx8XEajcbS1vEzBkZE3AEg6ciIOLIw6W2S/ht492wbjYi1wNp8PY+WpIgISc8gu37iQwWzPpLyV5Klkp4VEV8HkHQE8LDkSXAycKak+4F7gcnwCzLM+kpKYLwW+Kik/fLhu4HTut1QRKwB1nS7nJnVR8pfSW4FDskvWioifl1+WWZWR7MGRvGOzHwYgIiY9RqGmQ2WlFOS7YX+xWR/PflOOeWYWZ2lnJJcUByW9H7gitIqMrPamksDOkuAJy50IWZWfynXMIp3fC4C9gfeU2ZRZlZPKdcwik+X3g/8PCLuL6keM6uxlFOS90bEHXn304i4X9LHZ1/MzAZNSmA8vTiQN6CzopxyzKzOZgwMSedI2gocLOmevNsK/Bz4XGUVmlltzBgYEfEPEbEPcH5E7Jt3+0TEaEScU2GNZlYTM170lHRQRHwXuFTSoa3TI2J9qZWZWe10+ivJm4AzgAvaTAvg6FIqMrPa6tQexhl577ER8ZviNEmLS63KzGop5a8k30gcZ2YDrtM1jEcDjyV7vcAf81Br3vuS3R5uZkOm0zWMFwCnkr0a8Z8K47cCby+xJhsAzWYQWrmRtf7W6RrGxcDFkk6KiMsqrMn61EwhMdM8Do/+k/J4+2WSXkR2x+fiwng3oGMPSgmLdss4NPpLytOqF5Jds/gT4CNkjfneUnJd1gfmEhKd1uHwqL+Uv5IcERGvAX4VEe8Cngk8vtyyzKyOUgLj3vznDkkHAL8DxssryczqKqU9jCsl/T5wPrCe7C7PhXjPqvWxhTgdabdOn5bUW8pFz2brWpdJupLswudBpVZltVZGWBTX7dCor67a9IyI3+bvJbm0pHrMrMbm0ggwPHTXp5kNkbkGho8Zh1SZpyNVbsPmptOzJJ+nfTAIGC2tIjOrrU4XPd8/x2lmNqA6PUtyXZWFmFn9zfUahpkNIQeGmSVzYFhXqripyjdu1desgSHpy/mt4c3hR0i6utSqzKyWUo4wlkXE3c2BiPgV8MjSKjKz2koJjJ2SntAckHQgc7xxS9IqSRsk3SbpusL4SUnrJZ01l/Vatco8ZfDpSL2lPK36DuDrhR382WTvK+lKflrzYeCFEfEjScWjlElgJXCJpL0jYlu367dqRcSC35HpsKi/lKdVr8rffHY42V2eb4yIX85hW68A1kXEj/L13lWY1vzmBX5Oxay2Or2M+aD856HAE4A7gZ8CT2j36sQETwEeIelaSbdKek1h2jpgCpiKiK1zWLeZVUAzHQZK+reIOEPS19pMjojo6lWJktYAE8AxwF7AjcCLIuJ7Cctub/aPjIwsaTQa3WzaKjDX0xOfhtTT6OgoW7Zs2dEcjoilUPKrEiWtBk7PB/8TuCoitgPbJV0PHALMGhjNYgEmJib8DauhuVzTcFjU1/j4OI1GY2nr+FJflRgRayNieUQsBz4LHCVpd0lLgMOA76Ssx/pDRDzYzWceq6/KXpUYEd+RdBXwLWAn8JGI2Nx9ydYPHAiDKfVViRfwUGDcwxxflRgR55M1JmxmfcivSjSzZCnXMFa0eZbkveWVZGZ1lRIYx7Z5luS40ioys9pKCYxFkvZsDkjaC9izw/xmNqBSniX5BPBVSR8ju3X7NODiUqsys1pKeZbkPEmbyO7QFPCeiHB7GGZDKOUIg4j4IvDFkmsxs5pLaXHrcEnflLRN0n2SHpB0TxXFmVm9pFz0XAOcAnyf7KGx1wEfKrMoM6un1FOSH0haFBEPAB+TlPQsiZkNlpTA2CFpD2CDpPOAnwEPe4rNzAZfyinJq4FFwBuA7cDjgZPKLMrM6inlz6p35L33Au8qtxwzq7NOj7dvokPr4BFxcCkVmVltdTrCeHFlVZhZX+j0ePsdM00zs+HkG7fMLJlv3DKzZL5xy8yS+cYtM0uWeuPWbvjGLbOh182NW7/BN26ZDbVO71Z9Sf7msubwzZJuz7uTqynPzOqk0ynJW4ArCsN7AiuBVcCZJdZkZjXV6ZRkj4j4cWH46xHRABqSfNHTbAh1OsJ4RHEgIt5QGNy/nHLMrM46BcbNkk5vHSnp9cAt5ZVkZnXV6ZTkjcDlkl4BrM/HrSC7lnFCyXWZWQ11evjsLuAISUcDT89H/1dEXFNJZWZWOyn3YVwDOCTMLOlOTzMzwIFhZl1wYJhZMgeGmSWrLDAkvVnShrzbnLfcNZJPm5S0XtJZVdVjZt2rLDAi4vyIWB4Ry4FzgOsiYks+eZLsOZXDJe1dVU1m1p1enZKcAnyqMKz8ZxT6zaxmKg8MSUuAFwKXFUavA6aAqYjYWnVNZpZGETO+q6icDUovB14VEcd3scz2Zv/IyMiSRqNRSm1mlhkdHWXLli07msMRsRRKPsKQtLpwofOAfPQku56OzCoilja78fHxhS/UzHYxPj6+y37XHJ/UavhcRcRaYG1zWNJ+wHOAV5W5XTMrR9XXME4EvhQR22ed08xqp9LAiIiLImKyym2a2cLxnZ5mlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlqyywJC0n6TPS9oo6TZJf1aYNilpvaSzqqrHzLpX5RHGauDbEXEIsAq4QNIe+bRJYCVwuKS9K6zJzLpQZWAEsI8kAXsDW4D782kqzKM2y5pZDVQZGGuApwJ3ApuAv46Infm0dcAUMBURWyusycy6oIioZkPSycCRwJuAJwFfBg6JiHsSlt3e7B8ZGVnSaDRKq9PMYHR0lC1btuxoDkfEUij5CEPSakkbJG0gu4axLjI/AKaBg1LWExFLm934+HiJFZsZwPj4+C77XXN8qYEREWsjYnlELAe+CxwDIOlRwB8At5e5fTNbWLtXuK33ABdJ2kR2YfOtEfHLCrdvZvNUWWBExJ3A86vanpktPN/paWbJHBhmlsyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJKnuWZKEsW7YsxsbGHhyenp6mX24Xd63lcK0Lb3p6OhqNxsMOKPouMFpJ2l68173OXGs5XGt1fEpiZskcGGaWrO9PScysOj7CMLNkDgwzS+bAMLNktQwMSR+VdJekzYVxyyXdlDf5NyXpGfn4VzabAcy7nZKW59NW5fOeV4da82kHS7oxfzfLJkmL61irpDFJ9xZ+rxcWlqlVrYXpT5C0TdLZda1V0jMKv9ONkk6sstZ5i4jadcCzgUOBzYVxXwKOzfuPA65ts9wfAbcXhv8D2Au4ADio17WSNVj0LbLGjwFGgUU1rXWsOF/LempVa2H6ZcClwNl1rRVYAuye9z8GuKswXHqt8+1qeYQREdeTvbdkl9HAvnn/fmSvK2h1CvCpwvBu+XI7Kel9J13W+nzgWxGxMV+2EREP1LTWTmpXq6QTyNqIva3OtUbEjohovo9ncT5fZbXOW68Tq0Nqj7FrYj8V+BHwY+CnwIFtlvk/4A8Lwy8A1gMX1KFW4Czg48DVeV1vqXGtY8B24H+A64CjalzrUuBGshdkncuuRxi1qjWfdhhZsG0DTqy61nl9zl4X0MU/wL8AJ+X9LwO+0jL/YcCmOtcKnE32eoVlZIemNwLH1LTWPYHRvH9F/sXft6a1vh94Wd6/S2DUrdaWZZ4K3AIsrrreOX/OXhfQxT/Ar3noRjMB97TM/wHg7XWulewdshcV5vsb4M11rLXNctcCE3WsFbgB+GHe3U12evCGOtbaZrmvVf17nU9Xy2sYM7gTeE7efzTw/eYESbsBLwU+3YO62pmp1quBgyUtkbR7Ps+3e1BfUdtaJe0vaVHe/0TgyfT+PTJta42IoyJiLCLGgA8Cfx8Ra3pS4UNm+r2O5//2SDqQ7P08P+xFgXPS68SaIXU/BfwM+B3wE+C1wLOAW4GNwM3AisL8q4Cb+qTWV5Gdv24GzqtrrcBJeZ0byc6rj69rrS3LnUvFpyRd/l5fnf9eN+S/1xN68b2da+dnScwsWT+dkphZjzkwzCyZA8PMkjkwzCyZA8PMkjkwzCyZA6OPSHqUpE9Kul3Srflj8ifOssxY8bHrLrd3qqQDCsMfkfS0xGVXSbpyLttNXP8Bkj6T9y+XdNwc1nFu8VF4m50Do09IEnA5cH1EPDEiVpDdav64Ejd7KvBgYETE6yKi13emAhARd0bEyfngcrJHyK1kDoz+cTRwX0Q82JBNRNwRER+CB48kbpC0Pu+OaF1Bp3kkvSVv0GejpH+UdDIwAVySN/ayl6RrJU3k878wX8dGSV9N/RCSTsm3s1nS+wrjt0n6u3x9N0l6VD7+SfnwNyW9W9K2wmfZLGkP4N3Ay/M6X9565JDPN5b3v0PS/0r6Ctlt2RS2c1V+5HaDpINSP9NQ6fWtpu7SOuCvgA90mL6E/KlHsuc+pvL+MfKHojrMcyzwDWBJPjyS/7yWwoNRzWFgf7KnV8eL87fUswq4smXcAWSPfO9P1pjQNeS3RpO1A3F83n8e8M68/0rglLz/z4FtbT7XqcCawnbOZddH3Dfn868ANuW/h32BHzTnA74KPDnvPwy4ptf/5nXsdk9OFqsVSWvJnle4LyJWAr8HrFHWPOEDwFPaLDbTPM8FPhYROwAiorUxmFaHk50aTSfO37SSrOWpX+Sf4RKy1qouB+4jCwfInsF4Xt7/TOCEvP+TZI+yz9VRwGebn1PSFfnPvYEjgEuzMz8ge7zfWjgw+sdtZA+EARARqyUtA6byUW8Efg4cQnaq+Zs265hpHrFry0+z6Xb+4nIz+V3k/72Thdl8vpv3s+vp9uJCf7u6dwPujojl89jmUPA1jP5xDbBY0pmFcUsK/fsBP4uInWRPRC5qs46Z5vkScJqkJQCSRvLxW4F92qznRuA5ksZb5p/Nzflyy/JH508ha82rk5t4KCgnZ5intc4fkrWxiaRDgebbj68HTsyvx+wDHA8QEfcA05Jemi8jSYckfqah4sDoE/n/vieQ7XDTkm4BLgbems/yYeBPJd1Edqqxvc1q2s4TEVcBVwBTkjaQtQwGcBFwYfOiZ6GWXwBnAOskbSRrvLadYyT9pNmRXUc4h6zRmI3A+oj43Cwf/SzgTfnnfQxZwzStvgY8rXnRk6wx4JH8s5wJfC+ve31e64Z8nhsK63gl8Nr889wGvGSWuoaSH2+3WsuPeu6NiJA0SXYB1Dtzj/gahtXdCrILtSJrfu+03pYz3HyEYWbJfA3DzJI5MMwsmQPDzJI5MMwsmQPDzJL9P7AWIM9UmThPAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask_map.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the mask from a catalog of sources\n", "\n", "We can also build our list of regions from a list of catalog sources. Here we use the Fermi 4FGL catalog which we read using `~gammapy.catalog.SourceCatalog`." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.422304Z", "iopub.status.busy": "2021-11-22T21:09:29.421953Z", "iopub.status.idle": "2021-11-22T21:09:29.660892Z", "shell.execute_reply": "2021-11-22T21:09:29.661067Z" } }, "outputs": [], "source": [ "fgl = CATALOG_REGISTRY.get_cls(\"4fgl\")()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now select sources that are contained in the region we are interested in." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.663427Z", "iopub.status.busy": "2021-11-22T21:09:29.663074Z", "iopub.status.idle": "2021-11-22T21:09:29.667961Z", "shell.execute_reply": "2021-11-22T21:09:29.668150Z" } }, "outputs": [], "source": [ "inside_geom = geom.contains(fgl.positions)\n", "positions = fgl.positions[inside_geom]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now create the list of regions using our 0.3 degree radius a priori value. If the sources were extended, one would have to adapt the sizes to account for the larger size." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.670238Z", "iopub.status.busy": "2021-11-22T21:09:29.669961Z", "iopub.status.idle": "2021-11-22T21:09:29.671244Z", "shell.execute_reply": "2021-11-22T21:09:29.671429Z" } }, "outputs": [], "source": [ "exclusion_radius = Angle(\"0.3 deg\")\n", "regions = [\n", " CircleSkyRegion(position, exclusion_radius) for position in positions\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can build the mask map the same way as above." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.673087Z", "iopub.status.busy": "2021-11-22T21:09:29.672767Z", "iopub.status.idle": "2021-11-22T21:09:29.687961Z", "shell.execute_reply": "2021-11-22T21:09:29.688122Z" } }, "outputs": [], "source": [ "mask_map_catalog = ~geom.region_mask(regions)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.693821Z", "iopub.status.busy": "2021-11-22T21:09:29.693479Z", "iopub.status.idle": "2021-11-22T21:09:29.748686Z", "shell.execute_reply": "2021-11-22T21:09:29.748928Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEHCAYAAABbS3hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUcklEQVR4nO3deZRedX3H8feHpBASljoTXHBhRo8WtYWUTARBNAU3UI5wQB3cSlGwNLZFDy6o7cGliyDV1sRDW4/CUdQWiYhUwQVZrCxO0sQEtWoZccEjOhHJgiLk2z/u74Gb4Zlnfs/M3PvcZ/J5nfOcufv9PpN5Prn399zfvYoIzMxy7NHrAsysfzgwzCxbXweGpO29riFXP9UK/VVvP9UK/VdvWV8HhpnVy4FhZtnUb9+SLF26NIaGhgAYHx9neHi4twVlqqPWdevWzcl2li9f7t9thfqh3vHx8ZiYmHjYAcXCXhQzG0NDQ4yNjfW6jMaRNGfbWrduHf32H4nNrZGRkbZ/UD4lmQfmMizMOnFg9LmqwkKSg8gexoHRx/yBtro5MMwsmwPDzLI5MPqUT0esFxwYZpbNgWFm2RwYZpbNgWFm2RwYZpbNgdGn3NfDesGBYWbZHBhmls2B0ceqPi3xaY9N5sDoc1V9qB0W1o4DYx6Y6w+3w8Km4sCYJ+bqQ+6wsE4cGPPIbD/sDgubjgNjnpnJhz4iHBaWpe9uAmzTK3/4p+sG76Cwbjgw5jkHgs0ln5KYWTYHhpllc2CYWTYHhpllqz0wJK2Q9ICkU0rTRiWtl3R23fWYWb5aA0PSAuB9wDWTZo0CK4AjJO1TZ01mlq/uI4y/BC4H7po0vXWxQJSGzaxhagsMSY8FTgIuajN7LTAGjEXE1rpqMrPu1Hnh1geBt0bEA5OvPoyIS4BLplpR0vbW8MDAQFX1mVkyPj6+y+cuIpZAxYEhaRVwRhrdH/h0CoulwPGS7o+IK6bbTqtYgJGREV+6aFax4eFhJiYmlkyeXmlgRMQaYM3k6ZIuBq7KCQszaw5fh2Fm2XrS+SwiTuvFfs1sdnyEYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZsgJD0kGSnpuG95a0b7VlmVkTTRsYks4APgP8a5r0OOCKCmsys4bKOcJYBRwF3AMQEd8HHlllUWbWTDmB8duIuK81ImkhEN3uSNJKSb+WtCG9/rY0b1TSeklnd7tdM6vPwoxlrpf0dmBvSc8D/gL4/Az3d2NEvLjN9FFgBXCppH0iYtsMt29mFco5wngb8AtgE/B64AvAO+e4DqWfURo2s4aZ9ggjInYC/55es/VMSRuBO4FzIuK2NH0tMAZ8IiK2zsF+zKwCimjfHCFpEx3aKiLikK52JO0H7IyIbZKOB/45Ip6cue721vDAwMDiiYmJbnZtZl0aHBxky5YtO1rjEbEEOh9htNoaVqWfH08/XwnsePjiDydpFXBGGj0+Iu5MO/+CpA9LWhoRv5xuO61iAUZGRrpucDWz7gwPDzMxMbFk8vQpAyMi7gCQdFREHFWa9TZJ/w28e7qdRsQaYE3azqMlKSJC0jMo2k98qGDWR3K+JVki6VkR8XUASUcCD0ueDKcAZ0m6H7gXGI2pzofMrJFyAuO1wEcl7Z/G7wZO73ZHEbEaWN3tembWHDnfkqwDDk2NloqIX1dflpk10bSBUb4iM40DEBHTtmGY2fySc0qyvTS8iOLbk+9UU46ZNVnOKcmF5XFJ7weurKwiM2usmdxAZzHwxLkuxMyaL6cNo3zF5wLgAOA9VRZlZs2U04ZR7l16P/DziLi/onrMrMFyTkneGxF3pNdPI+J+SR+ffjUzm29yAuPp5ZF0A53l1ZRjZk02ZWBIOlfSVuAQSfek11bg58DnaqvQzBpjysCIiH+IiH2BCyJiv/TaNyIGI+LcGms0s4aYstFT0sER8V3gMkmHTZ4fEesrrczMGqfTtyRvAs4ELmwzL4BjKqnIzBqr0/0wzkyDx0XEb8rzJC2qtCoza6Scb0m+kTnNzOa5Tm0YjwYeS/F4gT/mobt570dxebiZ7WY6tWG8ADiN4tGI/1SavhV4e4U12TzQug1CO77RWv/q1IZxCXCJpJMj4vIaa7I+1Skk2i3n4Og/Od3bL5f0IoorPheVpvsGOgbkB8VU6zk4+kdOb9WLKNos/gT4CMXNfG+tuC7rAzMNiqm24+BovpxvSY6MiNcAv4qIdwHPBB5fbVlm1kQ5gXFv+rlD0oHA74Dh6kqyfjBXRxdVb9PmVs79MK6S9PvABcB6iqs85+I5q2bWZ3IaPVt317pc0lUUDZ8HV1qVNVqVRwKS3JbRYF3d0zMifpueS3JZRfWYWYPN5CbA8NBVn2Zzzm0ZzTXTwPAx427KH+bdW6e+JJ+nfTAIGKysIjNrrE6Nnu+f4Twzm6c69SW5vs5CzKz5ZtqGYWa7IQeGmWVzYFhX6rioyhduNde0gSHpy+nS8Nb4IyRdU2lVZtZIOUcYSyPi7tZIRPwKeGRlFZlZY+UExk5JT2iNSDqIGV64JWmlpA2SbpN0fWn6qKT1ks6eyXatXlWeMvh0pNlyequ+A/h66QP+bIrnlXQlndZ8GHhhRPxIUvkoZRRYAVwqaZ+I2Nbt9q1eETHnV306LJovp7fq1enJZ0dQXOX5xoj45Qz29QpgbUT8KG33rtK81l9e4H4qZo3V6WHMB6efhwFPAO4Efgo8od2jEzM8BXiEpOskrZP0mtK8tcAYMBYRW2ewbTOrgaY6DJT0bxFxpqSvtZkdEdHVoxIlrQZGgGOBvYGbgBdFxPcy1t3eGh4YGFg8MTHRza6tBjM9PfFpSDMNDg6yZcuWHa3xiFgCFT8qUdIq4Iw0+p/A1RGxHdgu6QbgUGDawGgVCzAyMuK/sAZqffBzg8NB0WzDw8NMTEwsmTy90kclRsSaiFgWEcuAzwJHS1ooaTFwOPCdnO1Y/4iIjmEw3XxrttoelRgR35F0NfAtYCfwkYjY3H3J1g8cCvNT7qMSL+ShwLiHGT4qMSIuoLiZsJn1IT8q0cyy5bRhLG/Tl+S91ZVkZk2VExjHtelLcnxlFZlZY+UExgJJe7VGJO0N7NVheTObp3L6knwC+Kqkj1Fcun06cEmlVZlZI+X0JTlf0iaKKzQFvCcifD8Ms91QzhEGEfFF4IsV12JmDZdzx60jJH1T0jZJ90l6QNI9dRRnZs2S0+i5GjgV+D5Fp7HXAR+qsigza6bcU5IfSFoQEQ8AH5OU1ZfEzOaXnMDYIWlPYIOk84GfAQ/rxWZm81/OKcmrgQXAG4DtwOOBk6ssysyaKedr1TvS4L3Au6otx8yarFP39k10uDt4RBxSSUVm1lidjjBeXFsVZtYXOnVvv2OqeWa2e/KFW2aWzRdumVk2X7hlZtl84ZaZZcu9cGsPfOGW2W6vmwu3foMv3DLbrXV6tupL0pPLWuO3SLo9vU6ppzwza5JOpyRvAa4sje8FrABWAmdVWJPZvCZp2ldTdTol2TMiflwa/3pETAATktzoadalboJAUiOfHtfpCOMR5ZGIeENp9IBqyjGbn2Zy1NDEo41OgXGLpDMmT5T0euDW6koym19m+6FvUmh0OiV5I3CFpFcA69O05RRtGSdWXJfZvDBXH/amnKJ06nx2F3CkpGOAp6fJ/xUR19ZSmVmfm+sjgyaERs51GNcCDgmzLlR1GtHr0Mi50tPMGqSXbRoODLM51qRGyrnmwDCzbA4MM8tWW2BIerOkDem1Od25ayDNG5W0XtLZddVjVoX5fDoCNQZGRFwQEcsiYhlwLnB9RGxJs0cp+qkcIWmfumoys+706pTkVOBTpfFWLEdp2MwapvbAkLQYeCFweWnyWmAMGIuIrXXXZGZ5VPdFIJJeDrwqIk7oYp3treGBgYHFExMTldRmNlt1tWFU/bkdHBxky5YtO0r7WwIVH2FIWlVq6DwwTR5l19ORaUXEktZreHh47gs1myO9vnR7rgwPD+/yuWtNz7pr+ExFxBpgTWtc0v7Ac4BXVblfM6tG3W0YJwFfiojt0y5pZo1Ta2BExMURMVrnPs3qNl9OS9rxlZ5mFagqNCLCvVXN5qP5eKThwDCr0FyGRhMCyIFhVrG5+KA3ISyg4q9VzWx2mhIULT7CMKvBTBormxYW4CMMs1qVQ6DdZeRNDIkyB4ZZjzQ9HNrxKYmZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZastMCTtL+nzkjZKuk3Sn5XmjUpaL+nsuuoxs+7VeYSxCvh2RBwKrAQulLRnmjcKrACOkLRPjTWZWRfqDIwA9lXxjPt9gC3A/WmeSsuozbpm1gB1BsZq4KnAncAm4K8jYmeatxYYA8YiYmuNNZlZFxQR9exIOgU4CngT8CTgy8ChEXFPxrrbW8MDAwOLJyYmKqvTzGBwcJAtW7bsaI1HxBKo+AhD0ipJGyRtoGjDWBuFHwDjwME524mIJa3X8PBwhRWbGcDw8PAun7vW9EoDIyLWRMSyiFgGfBc4FkDSo4A/AG6vcv9mNrcW1riv9wAXS9pE0bD51oj4ZY37N7NZqi0wIuJO4Pl17c/M5p6v9DSzbA4MM8vmwDCzbA4MM8vmwDCzbA4MM8vmwDCzbLX1JZkrS5cujaGhoQfHx8fH6ZfLxV1rNVzr3BsfH4+JiYmHHVD0XWBMJml7+Vr3JnOt1XCt9fEpiZllc2CYWba+PyUxs/r4CMPMsjkwzCybA8PMsjUyMCR9VNJdkjaXpi2TdHO65d+YpGek6a9s3QYwvXZKWpbmrUzLnt+EWtO8QyTdlJ7NsknSoibWKmlI0r2l3+tFpXUaVWtp/hMkbZN0TlNrlfSM0u90o6ST6qx11iKicS/g2cBhwObStC8Bx6Xh44Hr2qz3R8DtpfH/APYGLgQO7nWtFDcs+hbFzY8BBoEFDa11qLzcpO00qtbS/MuBy4BzmlorsBhYmIYfA9xVGq+81tm+GnmEERE3UDy3ZJfJwH5peH+KxxVMdirwqdL4Hmm9nVT0vJMua30+8K2I2JjWnYiIBxpaayeNq1XSiRT3iL2tybVGxI6IaD2PZ1FarrZaZ63XidUhtYfYNbGfCvwI+DHwU+CgNuv8H/CHpfEXAOuBC5tQK3A28HHgmlTXWxpc6xCwHfgf4Hrg6AbXugS4ieIBWeex6xFGo2pN8w6nCLZtwEl11zqr99nrArr4B/gX4OQ0/DLgK5OWPxzY1ORagXMoHq+wlOLQ9Cbg2IbWuhcwmIaXpz/8/Rpa6/uBl6XhXQKjabVOWuepwK3AorrrnfH77HUBXfwD/JqHLjQTcM+k5T8AvL3JtVI8Q/bi0nJ/A7y5ibW2We86YKSJtQI3Aj9Mr7spTg/e0MRa26z3tbp/r7N5NbINYwp3As9Jw8cA32/NkLQH8FLg0z2oq52par0GOETSYkkL0zLf7kF9ZW1rlXSApAVp+InAk+n9c2Ta1hoRR0fEUEQMAR8E/j4iVvekwodM9XsdTv/2SDqI4vk8P+xFgTPS68SaInU/BfwM+B3wE+C1wLOAdcBG4BZgeWn5lcDNfVLrqyjOXzcD5ze1VuDkVOdGivPqE5pa66T1zqPmU5Iuf6+vTr/XDen3emIv/m5n+nJfEjPL1k+nJGbWYw4MM8vmwDCzbA4MM8vmwDCzbA4MM8vmwOgjkh4l6ZOSbpe0LnWTP2madYbK3a673N9pkg4sjX9E0tMy110p6aqZ7Ddz+wdK+kwaXibp+Bls47xyV3ibngOjT0gScAVwQ0Q8MSKWU1xq/rgKd3sa8GBgRMTrIqLXV6YCEBF3RsQpaXQZRRdyq5gDo38cA9wXEQ/eyCYi7oiID8GDRxI3SlqfXkdO3kCnZSS9Jd3QZ6Okf5R0CjACXJpu9rK3pOskjaTlX5i2sVHSV3PfhKRT0342S3pfafo2SX+XtnezpEel6U9K49+U9G5J20rvZbOkPYF3Ay9Pdb588pFDWm4oDb9D0v9K+grFZdmU9nN1OnK7UdLBue9pt9LrS039ynsBfwV8oMP8xaRejxT9PsbS8BCpU1SHZY4DvgEsTuMD6ed1lDpGtcaBAyh6rw6Xl59Uz0rgqknTDqTo8n0Axc2EriVdGk1xH4gT0vD5wDvT8FXAqWn4z4Ftbd7XacDq0n7OY9cu7pvT8suBTen3sB/wg9ZywFeBJ6fhw4Fre/1v3sTXwuxksUaRtIaiv8J9EbEC+D1gtYrbEz4APKXNalMt81zgYxGxAyAiJt8MZrIjKE6NxjOXb1lBceepX6T3cCnF3aquAO6jCAco+mA8Lw0/EzgxDX+Soiv7TB0NfLb1PiVdmX7uAxwJXFac+QFF936bxIHRP26j6BAGQESskrQUGEuT3gj8HDiU4lTzN222MdUyYtc7P02n2+XL603ld5H+e6cIs9n8bd7Prqfbi0rD7ereA7g7IpbNYp+7Bbdh9I9rgUWSzipNW1wa3h/4WUTspOgRuaDNNqZa5kvA6ZIWA0gaSNO3Avu22c5NwHMkDU9afjq3pPWWpq7zp1LczauTm3koKEenWGZynT+kuMcmkg4DWk8/vgE4KbXH7AucABAR9wDjkl6a1pGkQzPf027FgdEn0v++J1J84MYl3QpcArw1LfJh4E8l3UxxqrG9zWbaLhMRVwNXAmOSNlDcGQzgYuCiVqNnqZZfAGcCayVtpLh5bTvHSvpJ60XRjnAuxU1jNgLrI+Jz07z1s4E3pff7GIob00z2NeBprUZPipsBD6T3chbwvVT3+lTrhrTMjaVtvBJ4bXo/twEvmaau3ZK7t1ujpaOeeyMiJI1SNID6w9wjbsOwpltO0VAritvvnd7bcnZvPsIws2xuwzCzbA4MM8vmwDCzbA4MM8vmwDCzbP8P2sgRNynK+PQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask_map_catalog.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the mask from statistically significant pixels in a dataset\n", "\n", "Here we want to determine an exclusion from the data directly. We will estimate the significance of the data using the `ExcessMapEstimator`, and exclude all pixels above a given threshold.\n", "\n", "Here we use the `MapDataset` taken from the Fermi data used above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We apply a significance estimation. We integrate the counts using a correlation radius of 0.4 degree and apply regular significance estimate. " ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.751370Z", "iopub.status.busy": "2021-11-22T21:09:29.751073Z", "iopub.status.idle": "2021-11-22T21:09:29.819630Z", "shell.execute_reply": "2021-11-22T21:09:29.819804Z" } }, "outputs": [], "source": [ "estimator = ExcessMapEstimator(\"0.4 deg\", selection_optional=[])\n", "result = estimator.run(dataset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we create the mask map by applying a threshold of 5 sigma to remove pixels." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.822069Z", "iopub.status.busy": "2021-11-22T21:09:29.821729Z", "iopub.status.idle": "2021-11-22T21:09:29.823259Z", "shell.execute_reply": "2021-11-22T21:09:29.822998Z" } }, "outputs": [], "source": [ "significance_mask = result[\"sqrt_ts\"] < 5.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because the `ExcessMapEstimator` returns NaN for masked pixels, we need to put the NaN values to `True` to avoid incorrectly excluding them. " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.825007Z", "iopub.status.busy": "2021-11-22T21:09:29.824691Z", "iopub.status.idle": "2021-11-22T21:09:29.826101Z", "shell.execute_reply": "2021-11-22T21:09:29.826286Z" } }, "outputs": [], "source": [ "invalid_pixels = np.isnan(result[\"sqrt_ts\"].data)\n", "significance_mask.data[invalid_pixels] = True" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.828136Z", "iopub.status.busy": "2021-11-22T21:09:29.827856Z", "iopub.status.idle": "2021-11-22T21:09:29.885089Z", "shell.execute_reply": "2021-11-22T21:09:29.885369Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEHCAYAAAA0+iR9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATZUlEQVR4nO3de5BkZX2H8ecLKMsul7gsxqiRGVMaNAms7OIFbxRqFKIVKDAuWiaUisZgUsQyXnMx5q4SNYJFSkvYGEIUIWhQwSiiqKCuGxCIiSZOMIaUmEGEXSAI/PJHn4beoaenZ3b79PTs86nqmnM/v5dhv3Oub6eqkKTd3R7jLkCSlgPDUJKY8DBMsn3cNYzKSm2b7Zo8K7ltvSY6DCVpVzEMJQnIpN1NXrduXU1NTQEwMzPD9PT0eAsakZXaNts1eVZS22ZmZmp2drbvQeBebRezs6amptiyZcu4y5A0gTZu3Jj55nmaLEkYhpIEGIaSBBiGkgQYhpIEGIaSBIwhDJMckeSeJCf2TNuUZGuS09quR5Kg5TBMsifwF8Clc2ZtAo4Anpxk3zZrkiRo/8jwN4ELgJvmTO8+CFk9w5LUmtbCMMkjgOOBs/rMvhDYAmypqtvaqkmSutp8He/dwBuq6p5kx4O/qtoMbG6xFknawUjDMMmpwCnN6AHA3zdBuA44NsndVXXRENu5rz+1tWvXjqBSSbuDmZmZHfKkqtZ0h0cahlV1JnDm3OlJzgEuHiYIm+3cV/DGjRsnq5sdScvG9PQ0s7Oza/rN8zlDSWJMXXhV1cnj2K8kzccjQ0nCMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCQA9hp3AVr5krS2r6pqbV9aWTwylCQMQ0kCWgzDJEcl+VGSq5vP7/fM25Rka5LT2qpHknq1fc3wiqp6fp/pm4AjgHOT7FtV21quS9JubrncQOleYa+eYWnRdvZmjTdgdl9tXzN8SpJrknwqyc/1TL8Q2AJsqarbWq5JkkhbfwmT7A/cW1XbkhwLvKeqHjPkutu7w2vXrl09Ozs7qjI1Am0+WrOzPDJc2Q488EBuvvnm27vjVbWmOzzUkWGSg5M8uxneJ8l+Q653aveGCXDftcCq+iTwoCTrhtlOVa3pfqanp4dZRZIeYHp6eoc86Z23YBgmOQX4KPDXzaRHAhcNs+OqOrOq1lfVeuDeNIcISZ7Y7NtDPEnLwjA3UE4Fngh8BaCqvp3koUvY14nAq5PcDdwBbCrPSSQtE8OE4f9V1V3d6z5J9qJz13dRquoM4IzFrqfJMUnXBuczXxv8u73yDXPN8PNJ3gzsk+Q5wPnAP462LElq1zBh+EbgB8C1wKuATwK/O8qiJKltC54mV9W9wPubjyStSPOGYZJrGXBtsKoOHUlFkjQGg44Mu+8Qn9r8/FDz8yXA7Q9cXLuTlXCzZDEW015vtkymecOwqm4ASPLUqnpqz6w3JvkS8LZRFydJbRnmBsqaJE/rjiQ5ElgzYHlJmjjDPGf4cuCDSQ5oxm8BXjayiiRpDIa5m/x14LCmo4VU1Y9GX5YktWvBMOztkboZB6CqvGYoacUY5jR5e8/wKjp3mb85mnKkydfvzrN3mJe/YU6TT+8dT/JO4OMjq0iSxmApPV2vBh69qwuRpHEa5pph75soewIHAX80yqIkqW3DXDPs/Ta7u4HvV9XdI6pHksZimNPkP66qG5rPf1fV3Uk+tPBqkjQ5hgnD3m+x63buumE05UjSeMwbhknelOQ24NAktzaf24DvAx9rrUJJasG8YVhVf1ZV+wHvqKr9m89+VXVgVb2pxRolaeQG9Wd4SFX9K3B+ksPnzq+qrSOtTJJaNOhu8muBVwKn95lXwNEjqUiSxmBQf4avbAaPqao7e+clWTXSqrRs7G6duI6K37q3/A1zN/nLQ06TpIk16Jrhw4BH0PmK0CcA3T9t+9N5JU+SVoxB1wyfC5wMPBL4y57ptwFvHmFNktS6QdcMNwObk5xQVRe0WJMktW6YLrwuSPJLdN5EWdUz3c5dVxhvlmh3tuANlCRnAS8CfpPOdcMXAgePuC5JatUwd5OPrKpfBX5YVX8IPAX46dGWJUntGiYM72h+3p7k4cCPgenRlSRJ7RumP8OLk/wE8A5gK523T94/yqIkqW3D3EDp9mp9QZKL6dxEOWSkVUlSyxb1HShV9X/N9yafP6J6JGkslvKFUHD/2yiStCIsNQx9u1zSijLo3eR/pH/oBThwZBVJ0hgMuoHyziXOk6SJM+jd5M+3WYja5at30o6Wes1QklYUw1CSMAwlCRiu15p/al7H644/JMmlI61Kklo2zLvJ66rqlu5IVf0wyUNHV5La0O+LiLypot3ZMKfJ9yZ5VHckycH40LWkFWaYI8O3AF9M0n3U5hl0vk9ZklaMBY8Mq+oS4HDgw8BHgA1VtaRrhkmOSnJ1kut7wpUkm5JsTXLaUrYrSTtr0Ot4h1TVvyY5vJl0Y/PzUUkeVVVbF7Oj5ibM+4DnVdV351x33AQcAZybZN+q2raYbUvSzhp0mvxaOqfDp/eZV8DRi9zXi4ELq+q7AFV1U8+87pX7wh5xJI3BoNfxutcFj6mqO3vnJVnVZ5WFPBZ4UJLLgf2A91TV3zTzLgS2AH9bVbctYdtaJO8cSzsa5gbKl+lcM1xo2jD72gA8C9gHuDLJVVX1re53NM+3YpLt3eG1a9cucreS1DEzM7NDnlTVmu7woGuGDwMeAeyT5Ancf/q6P7B6mB0nORU4pRn9CHBJVW0Htif5AnAY8K2FttNb8MaNG32sR9KSTE9PMzs7u6bfvEFHhs8FTgYeSee6YTcMbwXePMyOq+pM4EyAJI8DzkiyF/Bg4EnAu4bZjiSN2qBrhpuBzUlOqKoLdnZHVfXNJJcA3wDuBT5QVdft7HYlaVcY5g2UDX3eTf7jpeysqt5RVY+vqp+vqncvZRvaNarqAR9pdzZMGB4z991k4NiRVSRJYzBMGO6ZZO/uSJJ9gL0HLC9JE2eYR2v+FvhskrPpPBT9MgY8BiNJk2jBMKyqtye5ls7zgQH+aKnvJkvScjXMkSFV9SngUyOuRZLGZpierp+c5GtJtiW5K8k9SW5tozhJasswN1DOAE4Cvk3nNbpXAO8dZVGS1LZhT5P/PcmeVXUPcHaSL4+4Lklq1TBheHuSBwNXJ3k78D9A33f7JGlSDXOa/FJgT+A1wHbgp4ETRlmUJLVtmEdrbmgG7wD+cLTlaJz8xrzR8XXH5W9QF17XMuBb8Krq0JFUJEljMOjI8PmtVSFJYzaoC68b5psnSSuND11LEj50LUmAD11rAYu5C7q73Xn2DvHK4kPXksTwD13vgQ9dS1rBFvPQ9Z340LWkFWreI8Mkv9x873F3/CtJvtN8TmynPElqx6Ajw9cDm3rG9waOoHO98GzgoyOsSxNovhsKk3RjxZsiu69BYfjgqvqvnvEvVtUsMJvEGyiSVpRBN1Ae0jtSVa/pGT1oNOVI0ngMCsOvJDll7sQkrwK+OrqSJKl9g06Tfxu4KMmLga3NtA10rh0eN+K6JKlVgzpquAk4MsnRwM81kz9RVZe1UpkktWiY5wwvAwxATSzvEGsYw7yBIkkrnmEoSRiGkgQYhpIEDNmfobQzvIGhSeCRoSRhGEoSYBhKEmAYShJgGEoSYBhKEmAYShJgGEoSYBhKEmAYShJgGEoSYBhKEtBiGCb5nSRXN5/rktyTZG0zb1OSrUlOa6seSerVWhhW1Tuqan1VrQfeBHy+qm5uZm+i8wX1T06yb1s1SVLXuE6TTwLO6xlP87N6hiWpNa2HYZLVwPOAC3omXwhsAbZU1W1t1yRJ4+jc9QXAl3pOkamqzcDmMdQiScCIwzDJqcApzeixVXUjneuD582/Vt/tbO8Or127dtcVKGm3MjMzs0OeVNWa7vBIT5Or6szuTZOqujHJAcAzgY8tcjtrup/p6enRFCtpxZuent4hT3rntX3N8Hjg01W1fcElJalFrYZhVZ1TVZva3KckDcM3UCQJw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAEhVjbuGRVm3bl1NTU3dNz4zM8P09PT4ChoR2zVZbNdkmJmZqdnZ2b4HgRMXhnMl2V5Va8Zdx65muyaL7Zp8niZLEoahJAEr4DRZknYFjwwlCcNQkgDDUJKAZRqGST6Y5KYk1/VMW5/kqiRXJ9mS5InN9Jc007qfe5Osb+Yd1Sz79jE1ZQeLaVcz79AkVya5Psm1SVY10ye2XUmmktzR8/s6q2ediW1Xz/xHJdmW5HU90ya2XUme2PO7uibJ8T3rLKt27bSqWnYf4BnA4cB1PdM+DRzTDB8LXN5nvV8AvtMz/mFgH+B04JBJahewF/AN4LBm/EBgzxXQrqne5eZsZ2Lb1TP/AuB84HUroV3AamCvZvingJt6xpdVu3b2syyPDKvqC8DNcycD+zfDBwA39ln1JOC8nvE9mvXuBbKLy1y0RbbrF4FvVNU1zbqzVXVPM2+S2zXIRLcryXHAd4Dr56wzse2qqtur6u5m+qpmua5l1a6dNu40HvDXa4od/3I9Dvgu8F/AfwMH91nnP4Cf7xl/LrAVOH3c7Vlsu4DTgA8BlzZteP0KadcUsB34Z+DzwNNXSLvWAFcC+wJvZccjw4ltVzPvSXQCfhtw/HJu1079Nxl3AYv4Zf0VcEIz/CvAZ+Ys/yTg2nHXvavaBbwOmAHW0TlVuRJ41rjr3wXt2hs4sBne0Pzj23/c9e+Cdr0T+JVmeIcwXI6fxf77aqY/DvgqsGrc9Y/kv8m4C1jEL+tH3P+QeIBb5yz/LuDN4657V7UL2ASc07Pc7wG/M+76d9Xvq2e5y4GN465/F/y+rgD+s/ncQuc09DXjrn8Ev6/PLeff1858luU1w3ncCDyzGT4a+HZ3RpI9gBcCfz+GunbWfO26FDg0yeokezXL/MsY6luqvu1KclCSPZvhRwOPoXOdbVL0bVdVPb2qpqpqCng38KdVdcZYKlya+X5f083/fyQ5GPhZOoG/4uw17gL6SXIecBSwLsn3gD8ATgHe0/xi7gRe2bPKM4DvVdWy/ke1mHZV1Q+T/CXwNToXqT9ZVZ8YS+ELWOTv6xnA25LcDdwD/HpVzb2Yvyws4f/DibDIdj0NeGOSH9O5UfIbVfW/7Vc9er6bLEks04euJalthqEkYRhKEmAYShJgGEoSYBhKEmAYah5JfjLJ3yX5TpKvN12JHb/AOlO93UItcn8nJ3l4z/gHkjx+yHWPSnLxUvY75PYfnuSjzfD6JMcuYRtv7e3WS8uPYagHSBLgIuALVfXoqtpA5/XAR45wtycD94VhVb2iqpbFGzdVdWNVndiMrqfTxZVWGMNQ/RwN3FVV93W8WlU3VNV74b4jwCuSbG0+R87dwKBlkry+6az2miR/nuREYCNwbtOJ6D5JLk+ysVn+ec02rkny2WEbkeSkZj/XJfmLnunbkvxJs72rkvxkM/1nmvGvJXlbkm09bbkuyYOBtwEvaup80dwjvma5qWb4LUn+Lcln6LzGRs9+LmmOuK9IcsiwbdIIjfvlaD/L7wP8FvCuAfNX0/RcQufd4i3N8BTNy/8DljkG+DKwuhlf2/y8nJ4OALrjwEF0eraZ7l1+Tj1HARfPmfZwOl1SHUTntdPLgOOaeQW8oBl+O/C7zfDFwEnN8K8D2/q062TgjJ79vJUdu+u6rll+A3Bt899hf+Dfu8sBnwUe0ww/Cbhs3L9zP7U8303W8pLkTDrvqN5VVUcADwLOSOfrFe4BHttntfmWeTZwdlXdDlALv5f8ZDqn6zNDLt91BJ3emn/QtOFcOu9FXwTcRSf4AL4OPKcZfgpwXDP8d3S65VqqpwP/0G1nko83P/cFjgTO71yNADrdmmnMDEP1cz1wQnekqk5Nsg7Y0kz6beD7wGF0LrXc2Wcb8y0TduwteSGLXb53vfn8uJrDMjpBvTP/Du5mx8tNq3qG+9W9B3BLVa3fiX1qBLxmqH4uA1YleXXPtNU9wwcA/1NV9wIvBfbss435lvk08LIkqwGSrG2m3wbs12c7VwLPTDI9Z/mFfKVZb13TZdhJdHrWHuQq7v8jsGmeZebW+Z90vk+EJIcD0830LwDHN9c/9wNeAFBVtwIzSV7YrJMkhw3ZJo2QYagHaI6ajqMTJjNJvgpsBt7QLPI+4NeSXEXn9Hd7n830XaaqLgE+DmxJcjWdHr0BzgHO6t5A6anlB3S6k7owyTV0voSon2cl+V73Q+e63ZvodEZ6DbC1qj62QNNPA17btPen6HR4OtfngMd3b6DQ+QKotU1bXg18q6l7a1Pr1c0yV/Rs4yXAy5v2XA/88gJ1qQV24SU1mqPVO6qqkmyiczPFoNpNeM1Qut8GOjd9Qqfr/peNtxy1ySNDScJrhpIEGIaSBBiGkgQYhpIEGIaSBMD/A0rLMUhBzbU+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "significance_mask.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This method frequently yields isolated pixels or weakly significant features if one places the threshold too low. \n", "\n", "To overcome this issue, one can use `~skimage.filters.apply_hysteresis_threshold` . This filter allows to define two thresholds and mask only the pixels between the low and high thresholds if they are not continuously connected to a pixel above the high threshold. This allows to better preserve the structure of the excesses. \n", "\n", "Note that scikit-image is not a required dependency of gammapy, you might need to install it.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Masks operations\n", "\n", "If two masks share the same geometry it is easy to combine them with `Map` arithmetic.\n", "\n", "OR condition is represented by `|` operator :" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.887632Z", "iopub.status.busy": "2021-11-22T21:09:29.887315Z", "iopub.status.idle": "2021-11-22T21:09:29.942095Z", "shell.execute_reply": "2021-11-22T21:09:29.942287Z" }, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEHCAYAAABbS3hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATvUlEQVR4nO3de5RdZX3G8e9DUghJgDoTvKDCHF1a1BZSMhEE0RS8gbKEBergrRQFS2NbdOEFtV146UWQamvioq1LYSlqi0REquAFuVi5OKSJCWrVMuIlLNEzYpIJiiG//rHfAzvDmTPvmWSf2SfzfNbaa/Z9/85k5snee/b7bkUEZmY59prtAsysfzgwzCxbXweGpInZriFXP9UK/VVvP9UK/VdvWV8Hhpn1lgPDzLKp3/5KsmTJkhgaGgJgbGyMRqMxuwVl6qdaob/q7adaoT/qHRsbi2az+YgTivmzUcyuGBoaYnR0dLbLMNujDQ8Pq918X5KYWTYHhpllc2CYWTYHhpllc2CYWTYHhpllc2CYWTYHhpllc2CYWTYHhpllc2CYWTYHhpllc2CYWTYHhk1Jattg0eYwB4a1JYl+6yvFqtd3/WFYtVpnFQ4La8dnGPYQh4VNx4FhgMPC8jgwzGFh2XwPwxwUls1nGGaWzYFhZtkcGGaWzYFhZtkcGGaWzYFhZtl6HhiSlkt6UNJppXkjktZKOrfX9ZhZvp4GhqR5wPuB6yYtGgGWA0dJWtzLmswsX6/PMP4SuBK4d9L8VjvqKI2bWc30LDAkPR44BbikzeI1wCgwGhFbelWTmXWnl4+Gfwh4W0Q8OLljloi4DLhsqg0lTbTGBwYGqqrPzJKxsbGdfu8iYhFUHBiSVgJnpckDgM+ksFgCnChpe0RcNd1+WsUCDA8Pu+GDWcUajQbNZnPR5PmVBkZErAZWT54v6VLgmpywMLP68HMYZpZtVpq3R8QZs3FcM9s1PsMws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2xZgSHpEEnPS+P7Stqv2rLMrI6mDQxJZwGfBf41zXoCcFWFNZlZTeWcYawEjgE2A0TED4BHV1mUmdVTTmD8NiIeaE1Img9EtweStELSryWtS8PflpaNSFor6dxu92tmvTM/Y50bJb0D2FfS84G/AL4ww+PdHBEvaTN/BFgOXC5pcURsneH+zaxCOWcYbwd+AWwA3gB8EXjXbq5D6WuUxs2sZqY9w4iIHcC/p2FXPUvSemATcF5E3JnmrwFGgU9GxJbdcBwzq4Ai2t+OkLSBDvcqIuKwrg4k7Q/siIitkk4E/jkinpK57URrfGBgYGGz2ezm0GbWpcHBQcbHx7e1piNiEXQ+w2jda1iZvn4ifX0VsO2Rqz+SpJXAWWnyxIjYlA7+RUkfkbQkIn453X5axQIMDw93fcPVzLrTaDRoNpuLJs+fMjAi4m4AScdExDGlRW+X9N/Ae6Y7aESsBlan/TxWkiIiJD2T4v6JTxXM+kjOX0kWSXp2RHwDQNLRwCOSJ8NpwDmStgP3AyMx1fWQmdVSTmC8DviYpAPS9H3Amd0eKCJWAau63c7M6iPnryR3AIenm5aKiF9XX5aZ1dG0gVF+IjNNAxAR097DMLM9S84lyURpfAHFX0++W005ZlZnOZckF5enJX0AuLqyisystmbSgc5C4Em7uxAzq7+cexjlJz7nAQcC762yKDOrp5x7GOXWpduBn0fE9orqMbMay7kkeV9E3J2Gn0XEdkmfmH4zM9vT5ATGM8oTqQOdZdWUY2Z1NmVgSDpf0hbgMEmb07AF+Dnw+Z5VaGa1MWVgRMQ/RMR+wEURsX8a9ouIwYg4v4c1mllNTHnTU9KhEfE94ApJR0xeHhFrK63MzGqn019J3gycDVzcZlkAx1VSkZnVVqf+MM5OoydExG/KyyQtqLQqM6ulnL+SfDNznpnt4Trdw3gs8HiK1wv8MQ/35r0/xePhZjbHdLqH8ULgDIpXI/5Taf4W4B0V1mR7gFY3CJO5k7X+1ukexmXAZZJOjYgre1iT9ampQmKqdRwe/SenefuVkl5M8cTngtJ8d6BjD8kJi3bbODT6S05r1Uso7ln8CfBRis58b6+4LusDMwmJTvtweNRfzl9Jjo6I1wK/ioh3A88CnlhtWWZWRzmBcX/6uk3SQcDvgEZ1JZlZXeX0h3GNpN8HLgLWUjzluTves2p9bHdcjrTbpy9L6i3npmerd60rJV1DcePz0EqrslqrIizK+3Zo1FdXfXpGxG/Te0muqKgeM6uxmXQCDA8/9Wlmc8hMA8PnjHNUlZcjvTyGzUyntiRfoH0wCBisrCIzq61ONz0/MMNlZraH6tSW5MZeFmJm9TfTexhmNgc5MMwsmwPDutKLh6r84FZ9TRsYkr6SHg1vTT9K0nWVVmVmtZRzhrEkIu5rTUTEr4BHV1aRmdVWTmDskHRwa0LSIczwwS1JKyStk3SnpBtL80ckrZV07kz2a71V5SWDL0fqLae16juBb5R+wZ9D8b6SrqTLmo8AL4qIH0sqn6WMAMuByyUtjoit3e7feisidvsTmQ6L+stprXptevPZURRPeb4pIn45g2O9ElgTET9O+723tKz1kxe4nYpZbXV6GfOh6esRwMHAJuBnwMHtXp2Y4anAoyTdIOkOSa8tLVsDjAKjEbFlBvs2sx7QVKeBkv4tIs6W9PU2iyMiunpVoqRVwDBwPLAvcAvw4oj4fsa2E63xgYGBhc1ms5tDWw/M9PLElyH1NDg4yPj4+LbWdEQsgopflShpJXBWmvxP4NqImAAmJN0EHA5MGxitYgGGh4f9E1ZDrV/83OBwUNRbo9Gg2Wwumjy/0lclRsTqiFgaEUuBzwHHSpovaSFwJPDdnP1Y/4iIjmEw3XKrt569KjEivivpWuDbwA7goxGxsfuSrR84FPZMua9KvJiHA2MzM3xVYkRcRNGZsJn1Ib8q0cyy5dzDWNamLcn7qivJzOoqJzBOaNOW5MTKKjKz2soJjHmS9mlNSNoX2KfD+ma2h8ppS/JJ4GuSPk7x6PaZwGWVVmVmtZTTluRCSRsontAU8N6IcH8YZnNQzhkGEfEl4EsV12JmNZfT49ZRkr4laaukByQ9KGlzL4ozs3rJuem5Cjgd+AFFo7HXAx+usigzq6fcS5IfSpoXEQ8CH5eU1ZbEzPYsOYGxTdLewDpJFwL3AI9oxWZme76cS5LXAPOANwITwBOBU6ssyszqKefPqnen0fuBd1dbjpnVWafm7Rvo0Dt4RBxWSUVmVludzjBe0rMqzKwvdGrefvdUy8xsbvKDW2aWzQ9umVk2P7hlZtn84JaZZct9cGsv/OCW2ZzXzYNbv8EPbpnNaZ3erfrS9Oay1vRtku5Kw2m9Kc/M6qTTJclbgatL0/sAy4EVwDkV1mRmNdXpkmTviPhJafobEdEEmpJ809NsDup0hvGo8kREvLE0eWA15ZhZnXUKjNsknTV5pqQ3ALdXV5KZ1VWnS5I3AVdJeiWwNs1bRnEv4+SK6zKzGurU+Oxe4GhJxwHPSLP/KyKu70llZlY7Oc9hXA84JMws60lPMzPAgWFmXXBgmFk2B4aZZetZYEh6i6R1adiYeu4aSMtGJK2VdG6v6jGz7vUsMCLioohYGhFLgfOBGyNiPC0eoWincpSkxb2qycy6M1uXJKcDny5NK32N0riZ1UzPA0PSQuBFwJWl2WuAUWA0Irb0uiYzy6OIKd9VVM0BpVcAr46Ik7rYZqI1PjAwsLDZbFZSm5kVBgcHGR8f39aajohFUPEZhqSVpRudB6XZI+x8OTKtiFjUGhqNxu4v1Mx20mg0dvq9a83P6jV8piJiNbC6NS3pAOC5wKurPK6ZVaPX9zBOAb4cERPTrmlmtdPTwIiISyNipJfHNLPdx096mlk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFm2ngWGpAMkfUHSekl3Svqz0rIRSWslnduresyse708w1gJfCciDgdWABdL2jstGwGWA0dJWtzDmsysC70MjAD2kyRgMTAObE/LVFpHbbY1sxroZWCsAp4GbAI2AH8dETvSsjXAKDAaEVt6WJOZdUER0ZsDSacBxwBvBp4MfAU4PCI2Z2w70RofGBhY2Gw2K6vTzGBwcJDx8fFtremIWAQVn2FIWilpnaR1FPcw1kThh8AYcGjOfiJiUWtoNBoVVmxmAI1GY6ffu9b8SgMjIlZHxNKIWAp8DzgeQNJjgD8A7qry+Ga2e83v4bHeC1wqaQPFjc23RcQve3h8M9tFPQuMiNgEvKBXxzOz3c9PeppZNgeGmWVzYJhZNgeGmWVzYJhZNgeGmWVzYJhZtp61JdldlixZEkNDQw9Nj42N0S+Pi7vWarjW3W9sbCyazeYjTij6LjAmkzRRfta9zlxrNVxr7/iSxMyyOTDMLFvfX5KYWe/4DMPMsjkwzCybA8PMstUyMCR9TNK9kjaW5i2VdGvq8m9U0jPT/Fe1ugFMww5JS9OyFWndC+tQa1p2mKRb0rtZNkhaUMdaJQ1Jur/0fb2ktE2tai0tP1jSVknn1bVWSc8sfU/XSzqll7Xusoio3QA8BzgC2Fia92XghDR+InBDm+3+CLirNP0fwL7AxcChs10rRYdF36bo/BhgEJhX01qHyutN2k+tai0tvxK4AjivrrUCC4H5afxxwL2l6cpr3dWhlmcYEXETxXtLdpoN7J/GD6B4XcFkpwOfLk3vlbbbQUXvO+my1hcA346I9WnbZkQ8WNNaO6ldrZJOpugj9s461xoR2yKi9T6eBWm9ntW6y2Y7sTqk9hA7J/bTgB8DPwF+BhzSZpv/A/6wNP1CYC1wcR1qBc4FPgFcl+p6a41rHQImgP8BbgSOrXGti4BbKF6QdQE7n2HUqta07EiKYNsKnNLrWnfpc852AV38A/wLcGoafznw1UnrHwlsqHOtwHkUr1dYQnFqegtwfE1r3QcYTOPL0g/+/jWt9QPAy9P4ToFRt1onbfM04HZgQa/rnfHnnO0CuvgH+DUPP2gmYPOk9T8IvKPOtVK8Q/bS0np/A7yljrW22e4GYLiOtQI3Az9Kw30UlwdvrGOtbbb7eq+/r7sy1PIexhQ2Ac9N48cBP2gtkLQX8DLgM7NQVztT1XodcJikhZLmp3W+Mwv1lbWtVdKBkual8ScBT2H23yPTttaIODYihiJiCPgQ8PcRsWpWKnzYVN/XRvq3R9IhFO/n+dFsFDgjs51YU6Tup4F7gN8BPwVeBzwbuANYD9wGLCutvwK4tU9qfTXF9etG4MK61gqcmupcT3FdfVJda5203QX0+JKky+/ra9L3dV36vp48Gz+3Mx3clsTMsvXTJYmZzTIHhpllc2CYWTYHhpllc2CYWTYHhpllc2D0EUmPkfQpSXdJuiM1kz9lmm2Gys2uuzzeGZIOKk1/VNLTM7ddIemamRw3c/8HSfpsGl8q6cQZ7OOCclN4m54Do09IEnAVcFNEPCkillE8av6ECg97BvBQYETE6yNitp9MBSAiNkXEaWlyKUUTcquYA6N/HAc8EBEPdWQTEXdHxIfhoTOJmyWtTcPRk3fQaR1Jb00d+qyX9I+STgOGgctTZy/7SrpB0nBa/0VpH+slfS33Q0g6PR1no6T3l+ZvlfR3aX+3SnpMmv/kNP0tSe+RtLX0WTZK2ht4D/CKVOcrJp85pPWG0vg7Jf2vpK9SPJZN6TjXpjO3myUdmvuZ5pTZftTUQ94A/BXwwQ7LF5JaPVK0+xhN40OkRlEd1jkB+CawME0PpK83UGoY1ZoGDqRovdoorz+pnhXANZPmHUTR5PtAis6Eric9Gk3RD8RJafxC4F1p/Brg9DT+58DWNp/rDGBV6TgXsHMT941p/WXAhvR92B/4YWs94GvAU9L4kcD1s/1vXsdhfnayWK1IWk3RXuGBiFgO/B6wSkX3hA8CT22z2VTrPA/4eERsA4iIyZ3BTHYUxaXRWOb6Lcspep76RfoMl1P0VnUV8ABFOEDRBuP5afxZwMlp/FMUTdln6ljgc63PKenq9HUxcDRwRXHlBxTN+20SB0b/uJOiQRgAEbFS0hJgNM16E/Bz4HCKS83ftNnHVOuInXt+mk6365e3m8rvIv33ThFmu/KzuZ2dL7cXlMbb1b0XcF9ELN2FY84JvofRP64HFkg6pzRvYWn8AOCeiNhB0SJyXpt9TLXOl4EzJS0EkDSQ5m8B9muzn1uA50pqTFp/Orel7ZakpvOnU/Tm1cmtPByUI1OsM7nOH1H0sYmkI4DW249vAk5J92P2A04CiIjNwJikl6VtJOnwzM80pzgw+kT63/dkil+4MUm3A5cBb0urfAT4U0m3UlxqTLTZTdt1IuJa4GpgVNI6ip7BAC4FLmnd9CzV8gvgbGCNpPUUnde2c7ykn7YGivsI51N0GrMeWBsRn5/mo58LvDl93sdRdEwz2deBp7duelJ0BjyQPss5wPdT3WtTrevSOjeX9vEq4HXp89wJvHSauuYkN2+3WktnPfdHREgaobgB6l/mWeJ7GFZ3yyhu1Iqi+70zZ7ecuc1nGGaWzfcwzCybA8PMsjkwzCybA8PMsjkwzCzb/wOHoajH6BmJegAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask = mask_map | mask_map_catalog\n", "mask.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "AND condition is represented by `&` or `*` operators :" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:29.952042Z", "iopub.status.busy": "2021-11-22T21:09:29.951703Z", "iopub.status.idle": "2021-11-22T21:09:30.001816Z", "shell.execute_reply": "2021-11-22T21:09:30.002106Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEHCAYAAABbS3hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVCElEQVR4nO3de5gddX3H8feHUAgJl7obvKDCrj5a1BZSshEE0RS8gfIID6jBWykKlsa26IMX1PbBSy+CVFsTH9r6KDyK2iIRkSp4QS5WLm7SxAS1alnxAo/oiUguKEK+/WPmwGQ5e/Z3ds/MmXPO5/U88+zc53s2Zz6ZmZ35jSICM7MUu/W6ADPrHw4MM0vW14EhaXuva0jVT7VCf9XbT7VC/9Vb1NeBYWbVcmCYWTL1219JlixZEmNjYwBMTU0xPj7e24ISVVHrunXrurKeZcuW+Xdbon6od2pqKhqNxiMOKHbvRTHzMTY2xuTkZK/LqB1JXVvXunXr6Lf/SKy7JiYmWn6hfEoyALoZFmbtODD6XFlhIclBZI/gwOhj3qGtag4MM0vmwDCzZA6MPuXTEesFB4aZJXNgWMf8F5Th5cCwZNODwqExfBwYlmSmcHBoDBcHhs1qtlBwaAwPB0afqupZj9QwcGgMBweGdY1DY/A5MKyrHBqDzYHRx+r6CLpDY3A5MPqcQ8Oq5MAYAA4Nq4oDY0A4NKwKDowB4tCwsjkwBoxDw8rkwBhAEfFQVycOjf7Xd62GW2fahYZ3YOuUA2MIOShsrnxKMmR6GRZ1O0WyzjkwhojDwubLgTEkHBbWDZUHhqTlkh6UdEph3EpJ6yWdXXU9Vi6HxWCpNDAkLQDeD1wzbdJKYDlwhKS9q6xpGPTq6MJhMXiqPsL4S+By4O5p45vf6Cj0Wxc4LKybKgsMSY8HTgIuajF5LTAJTEbE1qpqGmS9bNnbYTG4qrwP40PA2yLiwelf5Ii4BLhkpgUlbW/2j4yMlFXfwPAFTpuvqampXfa7iFgMJQeGpFXAGfngfsBn8i/zEuB4SQ9ExBWzradZLMDExIS/kW04LKwbxsfHaTQai6ePLzUwImINsGb6eEkXA1elhIX1B4fFcPB9GAPE1yysbD15liQiTuvFdgeZw8Kq4IfP+pyvWViVfErSxxwWVjUHhnXMYTG8HBh9rBc7rsNiuDkw+lyVO7DDwhwYA6CKHdlhYeDAGBhl7tAOC2tyYAyQMnZsh4UVOTAGTDd3cIeFTefAGEDd2NEdFtaKA2NAzWeHd1jYTBwYA2wuO77DwtpxYAy4TgLAYWGzcWAMgZQgcFhYCgfGkGgXCA4LS+XAGCKtgsFhYZ1wYAyZYkA4LKxTbkBnCDkobK58hGFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyRwYZpbMgWFmyZICQ9JBkp6X9+8laZ9yyzKzOpo1MCSdAXwW+Nd81BOAK0qsycxqKuUIYxVwFHAvQET8AHh0mUWZWT2lBMZvI+L+5oCk3YGOW2CRtELSryVtyLu/LUxbKWm9pLM7Xa+ZVSelxa3rJb0D2EvS84G/AL4wx+3dGBEvaTF+JbAcuFTS3hGxbY7rN7MSpRxhvB34BbAJeAPwReBdXa5D+c8o9JtZzcx6hBERO4F/z7v5epakjcCdwDkRcVs+fi0wCXwyIrZ2YTtmVgLN1CCspE20uVYREYd0tCFpX2BnRGyTdDzwzxHxlMRltzf7R0ZGFjUajU42bWYdGh0dZcuWLTuawxGxGNofYTSvNazKf34i//kqYMcjZ38kSauAM/LB4yPiznzjX5T0EUlLIuKXs62nWSzAxMSEm7w2K9n4+DiNRmPx9PEzBkZE3AEg6aiIOKow6e2S/ht4z2wbjYg1wJp8PY+VpIgISc8ku37iQwWzPpLyV5LFkp4dEd8AkHQk8IjkSXAKcJakB4D7gJXhF2SY9ZWUwHgd8DFJ++XD9wCnd7qhiFgNrO50OTOrj5S/kqwDDs0vWioifl1+WWZWR7MGRvGOzHwYgIiY9RqGmQ2WlFOS7YX+hWR/PfluOeWYWZ2lnJJcWByW9AHgytIqMrPamksDOouAJ3W7EDOrv5RrGMU7PhcA+wPvLbMoM6unlGsYxadLHwB+HhEPlFSPmdVYyinJ+yLijrz7WUQ8IOkTsy9mZoMmJTCeURzIG9BZVk45ZlZnMwaGpHMlbQUOkXRv3m0Ffg58vrIKzaw2ZgyMiPiHiNgHuCAi9s27fSJiNCLOrbBGM6uJGS96Sjo4Ir4HXCbpsOnTI2J9qZWZWe20+yvJm4EzgQtbTAvgmFIqMrPaatcexpl573ER8ZviNEkLS63KzGop5a8k30wcZ2YDrt01jMcCjyd7vcAf83Br3vuS3R5uZkOm3TWMFwKnkb0a8Z8K47cC7yixJhsAzWYQWnFDa/2r3TWMS4BLJJ0cEZdXWJP1qXYh0Wo+B0f/SXm8/XJJLya743NhYbwb0DEgPShmWs7B0T9Snla9iOyaxZ8AHyVrzPfWkuuyPjDXoJhpPQ6O+kv5K8mREfFa4FcR8W7gWcATyy3LzOooJTDuy3/ukHQA8DtgvLySrB906+ii7HVad6W0h3GVpN8HLgDWk93l2Y33rJpZn0m56NlsXetySVeRXfg8uNSqrNbKPBKQ5GsZNdZRm54R8dv8vSSXlVSPmdXYXBoBhofv+jTrOl/LqK+5BoaPGYeUd+bh1u5Zki/QOhgEjJZWkZnVVruLnh+Y4zQzG1DtniW5vspCzKz+5noNw8yGkAPDzJI5MKwjVdxU5Ru36mvWwJD0lfzW8ObwoyRdU2pVZlZLKUcYSyLinuZARPwKeHRpFZlZbaUExk5JBzYHJB3EHG/ckrRC0gZJt0m6vjB+paT1ks6ey3qtWmWeMvh0pN5SnlZ9J/CNwg7+HLL3lXQkP635CPCiiPixpOJRykpgOXCppL0jYlun67dqRUTX7/p0WNRfytOqV+dvPjuC7C7PN0XEL+ewrVcCayPix/l67y5Ma37zAj+nYlZb7V7GfHD+8zDgQOBO4GfAga1enZjgqcCjJF0naZ2k1xamrQUmgcmI2DqHdZtZBTTTYaCkf4uIMyV9vcXkiIiOXpUoaTUwARwL7AXcBLw4Ir6fsOz2Zv/IyMiiRqPRyaatAnM9PfFpSD2Njo6yZcuWHc3hiFgMJb8qUdIq4Ix88D+BqyNiO7Bd0g3AocCsgdEsFmBiYsLfsBqayzUNh0V9jY+P02g0Fk8fX+qrEiNiTUQsjYilwOeAoyXtLmkRcDjw3ZT1WH+IiIe6+cxj9VXZqxIj4ruSrga+DewEPhoRmzsv2fqBA2Ewpb4q8UIeDox7meOrEiPiArLGhM2sD/lViWaWLOUaxrIWz5K8r7ySzKyuUgLjuBbPkhxfWkVmVlspgbFA0p7NAUl7AXu2md/MBlTKsySfBL4m6eNkt26fDlxSalVmVkspz5KcL2kT2R2aAt4bEW4Pw2wIpRxhEBFfAr5Uci1mVnMpLW4dIelbkrZJul/Sg5LuraI4M6uXlIueq4FTgR+QPTT2euDDZRZlZvWUekryQ0kLIuJB4OOSkp4lMbPBkhIYOyTtAWyQdD5wF/CIp9jMbPClnJK8BlgAvBHYDjwROLnMosysnlL+rHpH3nsf8O5yyzGzOmv3ePsm2rQOHhGHlFKRmdVWuyOMl1RWhZn1hXaPt98x0zQzG06+ccvMkvnGLTNL5hu3zCyZb9wys2SpN27thm/cMht6ndy49Rt845bZUGv3btWX5m8uaw7fIun2vDulmvLMrE7anZK8FbiyMLwnsBxYAZxVYk1mA03SrF1dtTsl2SMiflIY/kZENICGJF/0NOtQJ0EgqZZvj2t3hPGo4kBEvLEwuH855ZgNprkcNdTxaKNdYNwi6YzpIyW9Abi1vJLMBst8d/o6hUa7U5I3AVdIeiWwPh+3jOxaxokl12U2ELq1s9flFKXdw2d3A0dKOgZ4Rj76vyLi2koqM+tz3T4yqENopNyHcS3gkDDrQFmnEb0OjZQ7Pc2sRnp5TcOBYdZldbpI2W0ODDNL5sAws2SVBYakt0jakHeb85a7RvJpKyWtl3R2VfWYlWGQT0egwsCIiAsiYmlELAXOBa6PiC355JVkz6kcIWnvqmoys8706pTkVODTheFmLEeh38xqpvLAkLQIeBFweWH0WmASmIyIrVXXZGZpVPVNIJJeAbw6Ik7oYJntzf6RkZFFjUajlNrM5quqaxhl77ejo6Ns2bJlR2F7i6HkIwxJqwoXOg/IR69k19ORWUXE4mY3Pj7e/ULNuqTXt253y/j4+C77XXN8UqvhcxURa4A1zWFJ+wHPBV5d5nbNrBxVX8M4CfhyRGyfdU4zq51KAyMiLo6IlVVu06xqg3Ja0orv9DQrQVmhERF+WtVsEA3ikYYDw6xE3QyNOgSQA8OsZN3Y0esQFlDyn1XNbH7qEhRNPsIwq8BcLlbWLSzARxhmlSqGQKvbyOsYEkUODLMeqXs4tOJTEjNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNL5sAws2QODDNLVllgSNpP0hckbZR0m6Q/K0xbKWm9pLOrqsfMOlflEcYq4DsRcSiwArhQ0h75tJXAcuAISXtXWJOZdaDKwAhgH2XvuN8b2AI8kE9TYR61WNbMaqDKwFgNPA24E9gE/HVE7MynrQUmgcmI2FphTWbWAUVENRuSTgGOAt4MPBn4CnBoRNybsOz2Zv/IyMiiRqNRWp1mBqOjo2zZsmVHczgiFkPJRxiSVknaIGkD2TWMtZH5ITAFHJyynohY3OzGx8dLrNjMAMbHx3fZ75rjSw2MiFgTEUsjYinwPeBYAEmPAf4AuL3M7ZtZd+1e4bbeC1wsaRPZhc23RcQvK9y+mc1TZYEREXcCL6hqe2bWfb7T08ySOTDMLJkDw8ySOTDMLJkDw8ySOTDMLJkDw8ySVfYsSbcsWbIkxsbGHhqempqiX24Xd63lcK3dNzU1FY1G4xEHFH0XGNNJ2l68173OXGs5XGt1fEpiZskcGGaWrO9PScysOj7CMLNkDgwzS+bAMLNktQwMSR+TdLekzYVxSyXdnDf5Nynpmfn4VzWbAcy7nZKW5tNW5POeX4da82mHSLopfzfLJkkL61irpDFJ9xV+rxcVlqlVrYXpB0raJumcutYq6ZmF3+lGSSdVWeu8RUTtOuA5wGHA5sK4LwPH5f3HA9e1WO6PgNsLw/8B7AVcCBzc61rJGiz6NlnjxwCjwIKa1jpWnG/aempVa2H65cBlwDl1rRVYBOye9z8OuLswXHqt8+1qeYQRETeQvbdkl9HAvnn/fmSvK5juVODTheHd8uV2UtL7Tjqs9QXAtyNiY75sIyIerGmt7dSuVkknkrURe1uda42IHRHRfB/Pwny+ymqdt14nVpvUHmPXxH4a8GPgJ8DPgINaLPN/wB8Whl8IrAcurEOtwNnAJ4Br8rreWuNax4DtwP8A1wNH17jWxcBNZC/IOo9djzBqVWs+7XCyYNsGnFR1rfP6nL0uoIN/gH8BTs77Xw58ddr8hwOb6lwrcA7Z6xWWkB2a3gQcW9Na9wRG8/5l+Rd/35rW+gHg5Xn/LoFRt1qnLfM04FZgYdX1zvlz9rqADv4Bfs3DN5oJuHfa/B8E3lHnWsneIXtxYb6/Ad5Sx1pbLHcdMFHHWoEbgR/l3T1kpwdvrGOtLZb7etW/1/l0tbyGMYM7gefm/ccAP2hOkLQb8DLgMz2oq5WZar0GOETSIkm75/N8pwf1FbWsVdL+khbk/U8CnkLv3yPTstaIODoixiJiDPgQ8PcRsbonFT5spt/reP5vj6SDyN7P86NeFDgnvU6sGVL308BdwO+AnwKvA54NrAM2ArcAywrzrwBu7pNaX012/roZOL+utQIn53VuJDuvPqGutU5b7jwqPiXp8Pf6mvz3uiH/vZ7Yi+/tXDs/S2JmyfrplMTMesyBYWbJHBhmlsyBYWbJHBhmlsyBYWbJHBh9RNJjJH1K0u2S1uWPyZ80yzJjxceuO9zeaZIOKAx/VNLTE5ddIemquWw3cf0HSPps3r9U0vFzWMd5xUfhbXYOjD4hScAVwA0R8aSIWEZ2q/kTStzsacBDgRERr4+IXt+ZCkBE3BkRp+SDS8keIbeSOTD6xzHA/RHxUEM2EXFHRHwYHjqSuFHS+rw7cvoK2s0j6a15gz4bJf2jpFOACeDSvLGXvSRdJ2kin/9F+To2Svpa6oeQdGq+nc2S3l8Yv03S3+Xru1nSY/LxT86HvyXpPZK2FT7LZkl7AO8BXpHX+YrpRw75fGN5/zsl/a+kr5Ldlk1hO1fnR243Sjo49TMNlV7fauourQP+Cvhgm+mLyJ96JHvuYzLvHyN/KKrNPMcB3wQW5cMj+c/rKDwY1RwG9id7enW8OP+0elYAV00bdwDZI9/7kzUmdC35rdFk7UCckPefD7wr778KODXv/3NgW4vPdRqwurCd89j1EffN+fzLgE3572Ff4IfN+YCvAU/J+w8Hru31v3kdu92Tk8VqRdIasucV7o+I5cDvAauVNU/4IPDUFovNNM/zgI9HxA6AiJjeGMx0R5CdGk0lzt+0nKzlqV/kn+FSstaqrgDuJwsHyJ7BeH7e/yzgxLz/U2SPss/V0cDnmp9T0pX5z72BI4HLsjM/IHu836ZxYPSP28geCAMgIlZJWgJM5qPeBPwcOJTsVPM3LdYx0zxi15afZtPp/MXlZvK7yP97Jwuz+Xw3H2DX0+2Fhf5Wde8G3BMRS+exzaHgaxj941pgoaSzCuMWFfr3A+6KiJ1kT0QuaLGOmeb5MnC6pEUAkkby8VuBfVqs5ybguZLGp80/m1vy5Zbkj86fStaaVzs383BQrpxhnul1/oisjU0kHQY03358A3BSfj1mH+AEgIi4F5iS9LJ8GUk6NPEzDRUHRp/I//c9kWyHm5J0K3AJ8LZ8lo8AfyrpZrJTje0tVtNynoi4GrgSmJS0gaxlMICLgYuaFz0LtfwCOBNYK2kjWeO1rRwr6afNjuw6wrlkjcZsBNZHxOdn+ehnA2/OP+/jyBqmme7rwNObFz3JGgMeyT/LWcD387rX57VuyOe5sbCOVwGvyz/PbcBLZ6lrKPnxdqu1/KjnvogISSvJLoB6Z+4RX8OwultGdqFWZM3vnd7bcoabjzDMLJmvYZhZMgeGmSVzYJhZMgeGmSVzYJhZsv8HjBiJKVWFiloAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask_map &= mask_map_catalog\n", "mask_map.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The NOT operator is represented by `~` symbol:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.008007Z", "iopub.status.busy": "2021-11-22T21:09:30.005933Z", "iopub.status.idle": "2021-11-22T21:09:30.061697Z", "shell.execute_reply": "2021-11-22T21:09:30.061956Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEHCAYAAAA0+iR9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATGElEQVR4nO3deZCkdX3H8fcHUJflShZIDB7smNLgEdiw4oFHKNQoHiUUqIuWCWWCiUFTxFKj5jLmVomaYMpULHXjFUUIGqJoFFES8Rg3IJAYTdzgWWJ2VdhdCALf/NFPw7NjT0/P7PbT073vV1XXPGc/39/2zmee89epKiRpX7ffpAuQpNXAMJQkpjwMk+ycdA3jMqtts13TZ5bb1jbVYShJe4thKElApu1qcpLpKljSqrFu3bratm3bwJ1A9wwl7TPm5uay2DzDUJIwDCUJMAwlCTAMJQkwDCUJMAwlCZhAGCY5IcntSc5oTduUZEuSc7uuR5Kg4zBMsj/w58BHFszaBJwAPCLJwV3WJEnQ/Z7hi4ALgRsWTO/fCFmtYUnqTGdhmORewGnAmwfMvgiYB+ar6qauapKkvgM63NYbgN+qqtuT3Xf+qmozsLnDWiRpN2PtqCHJOcDZzehh3HUIfASwC3h+VV08wvu0+1NbuzdrlLTvWLduHdu3b9/VH6+qg/rDE+m1JsnbgUuq6v0rWNdeayStyMaNG5mfnx94XcL7DCWJbs8Z3qmqzprEdiVpMe4ZShKGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAFwwKQL0Oyrqs62laSzbWm2uGcoSRiGkgR0GIZJTkrygyRXNa/fa83blGRLknO7qkeS2ro+Z3hFVT11wPRNwAnAu5IcXFU7Oq5L0j5utVxA6Z/1rtawtGx7erHGCzD7rq7PGT4yydVJPpzkwa3pFwHzwHxV3dRxTZJEurrtIcmhwB1VtSPJk4E3VtX9R1x3Z2t07VgK1Nh0eWvNnnLPcLatW7eO7du37+qPV9VB/eGR9gyTHJ3k8c3wgUkOGXG9c/oXTIA7zwVW1YeAuyU5YpT3qaqD+q9RlpekQebm5hbNkyXDMMnZwPuBv2km3Ru4eJQNV9WbqmpDVW0A7kjzZzfJw5ptbxu5FZI0RqNcQDkHeBjwWYCq+kqSn1jBts4AXpDkNuBmYFNN0/GTpJk2Shj+X1Xd2j+XkuQAeld9l6WqzgfOX+56mh6z8LdtsTZ4LnH2jXLO8JNJXgkcmOQJwAXAP463LEnq1ihh+HLgu8A1wK8CHwJ+Z5xFSVLXljxMrqo7gL9tXpI0kxYNwyTXMOTcYFUdO5aKJGkChu0Z9p8hPqf5+Y7m53OAXT+6uPYls3CxZDmW014vtkynRcOwqq4HSPKoqnpUa9bLk/wr8OpxFydJXRnlAspBSR7dH0lyIuCTIJJmyij3Gf4y8NYkhzXj3weeN7aKJGkCRrma/AXguKajhVTVD8ZfliR1a8kwbPdI3YwDUFWeM5Q0M0Y5TG53n7WG3lXm/xhPOdL0G3Tl2SvMq98oh8nntceTvA744NgqkqQJWElP12uB++3tQiRpkkY5Z9h+EmV/4EjgD8dZlCR1bZRzhu1vs7sN+E5V3TameiRpIkY5TP6jqrq+eX2zqm5L8o6lV5Ok6TFKGLa/xa7fuevG8ZQjSZOxaBgmeUWSm4Bjk9zYvG4CvgN8oLMKJakDi4ZhVf1pVR0CvLaqDm1eh1TV4VX1ig5rlKSxG9af4TFV9SXggiTHL5xfVVvGWpkkdWjY1eQXA88Hzhswr4CTx1KRJE3AsP4Mn98MnlJVt7TnJVkz1qq0auxrnbiOi9+6t/qNcjX50yNOk6SpNeyc4T2Be9H7itCfA/p/wg6l90ieJM2MYecMnwicBdwb+IvW9JuAV46xJknq3LBzhpuBzUlOr6oLO6xJkjo3ShdeFyZ5Cr0nUda0ptu564zxYon2ZUteQEnyZuBZwIvonTd8BnD0mOuSpE6NcjX5xKr6ReB7VfUHwCOB+4y3LEnq1ihheHPzc1eSo4AfAnPjK0mSujdKf4aXJPkx4LXAFnpPn/ztOIuSpK6NcgGl36v1hUkuoXcR5ZixViVJHVvWd6BU1f8135t8wZjqkaSJWMkXQsFdT6NI0kxYaRh6Q5qkmTLs2eR/ZHDoBTh8bBVJ0gQMu4DyuhXOk6SpM+zZ5E92WYi65aN30u5Wes5QkmaKYShJGIaSBIzWa80/N4/j9cd/PMlHxlqVJHVslGeTj6iq7/dHqup7SX5ifCWpC4O+iMiLKtqXjXKYfEeS+/ZHkhyNN11LmjGj7Bn+NvAvSfq32jyW3vcpS9LMWHLPsKouBY4H3gu8D9hYVSs6Z5jkpCRXJbmuFa4k2ZRkS5JzV/K+krSnhj2Od0xVfSnJ8c2kbzU/75vkvlW1ZTkbai7C/DXwpKr62oLzjpuAE4B3JTm4qnYs570laU8NO0x+Mb3D4fMGzCvg5GVu69nARVX1NYCquqE1r382v7BHHEkTMOxxvP55wVOq6pb2vCRrBqyylAcAd0tyOXAI8Maq+rtm3kXAPPDOqrppBe+tZfLKsbS7LPVLkWRLVR2/1LQlN5ScDzwUeBxwIHAl8JSq+vII6+5sja5dznY1mGG4Ogy6xUnjs27dOrZv376rP15VB/WHh50zvCdwL+DAJD/HXYevhzJiICU5Bzi7GX0fcGlV7QR2JvkUcBywZBi2C07ib7GkFZmbm2Pbtm0HDZo37JzhE4GzgHvTO2/YD8MbgVeOsuGqehPwJoAkDwTOT3IAcHfg4cDrR3kfSRq3YecMNwObk5xeVRfu6Yaq6j+SXAp8EbgDeEtVXbun7ytJe8MoT6BsHPBs8h+tZGNV9dqqelBVPaSq3rCS99DekeRHXtK+bJQwPGXhs8nAk8dWkSRNwChhuH+Se/RHkhwI3GPI8pI0dUZ5NvmdwMeTvI3eTdHPAzaPtSpJ6tiSYVhVr0lyDb37AwP84UqfTZak1WqUPUOq6sPAh8dciyRNzCg9XT8iyeeT7Ehya5Lbk9zYRXGS1JVRLqCcD5wJfIXeY3S/AvzVOIuSpK6Nepj8X0n2r6rbgbcl+fSY65KkTo0ShruS3B24KslrgG8DA5/tk6RpNcph8nOB/YEXAjuB+wCnj7MoSeraKLfWXN8M3gz8wXjL0ST5jXnj4+OOq9+wLryuYci34FXVsWOpSJImYNie4VM7q0KSJmxYF17XLzZPkmaNN11LEt50LUmAN11rCcu5CrqvXXn2CvFs8aZrSWL0m673w5uuJc2w5dx0fQvedC1pRi26Z5jk6c33HvfHP5vkq83rjG7Kk6RuDNszfBmwqTV+D+AEeucL3wa8f4x1aQotdkFhmi6seFFk3zUsDO9eVV9vjf9LVW0DtiXxAoqkmTLsAsqPt0eq6oWt0SPHU44kTcawMPxskrMXTkzyq8DnxleSJHVv2GHybwIXJ3k2sKWZtpHeucNTx1yXJHVqWEcNNwAnJjkZeHAz+Z+q6rJOKpOkDo1yn+FlgAGoqeUVYo1ilCdQJGnmGYaShGEoSYBhKEnAiP0ZSnvCCxiaBu4ZShKGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAR2GYZKXJrmqeV2b5PYk65p5m5JsSXJuV/VIUltnYVhVr62qDVW1AXgF8Mmq2t7M3kTvC+ofkeTgrmqSpL5JHSafCbynNd7v46law5LUmc7DMMla4EnAha3JFwHzwHxV3dR1TZI0ic5dnwb8a+sQmaraDGyeQC2SBIx5zzDJOa2LJkc1kzex+yHyKO+zs//a+1VK2lds3bp10TxJVXVWSJLDgK3AfapqRcGWpLuCJc2UjRs3Mj8/P/C6RNfnDE8DPrrSIJSkcek0DKvq7VW1qcttStIofAJFkjAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQmAAyZdwHIdfvjhrF+//s7xrVu3Mjc3N7mCxsR2TRfbNR22bt1ai81L1aLzpkKSnVV10KTr2Nts13SxXdPPw2RJwjCUJGAGDpMlaW9wz1CSMAwlCTAMJQlYpWGY5K1JbkhybWvahiSfSXJVkvkkD2umP6eZ1n/dkWRDM++kZtnXTKgpu1lOu5p5xya5Msl1Sa5JsqaZPrXtSrI+yc2tz+vNrXWmtl2t+fdNsiPJS1rTprZdSR7W+qyuTnJaa51V1a49VlWr7gU8FjgeuLY17aPAKc3wk4HLB6z3s8BXW+PvBQ4EzgOOmaZ20bsh/ovAcc344cD+M9Cu9e3lFrzP1LarNf9C4ALgJbPQLmAtcEAz/FPADa3xVdWuPX2tyj3DqvoUsH3hZODQZvgw4FsDVj0TeE9rfL9mvTuA7OUyl22Z7foF4ItVdXWz7raqur2ZN83tGmaq25XkVOCrwHUL1pnadlXVrqq6rZm+plmub1W1a49NOo2H/PVaz+5/uR4IfA34OvBN4OgB6/w38JDW+BOBLcB5k27PctsFnAu8A/hI04aXzUi71gM7gX8DPgk8ZkbadRBwJXAw8Cp23zOc2nY18x5OL+B3AKet5nbt0b/JpAtYxof1l8DpzfAzgY8tWP7hwDWTrntvtQt4CbAVOILeocqVwOMmXf9eaNc9gMOb4Y3NL9+hk65/L7TrdcAzm+HdwnA1vpb7+9VMfyDwOWDNpOsfy7/JpAtYxof1A+66STzAjQuWfz3wyknXvbfaBWwC3t5a7neBl066/r31ebWWuxx46KTr3wuf1xXA/zSv79M7DH3hpOsfw+f1idX8ee3Ja1WeM1zEt4Cfb4ZPBr7Sn5FkP+AZwN9PoK49tVi7PgIcm2RtkgOaZf59AvWt1MB2JTkyyf7N8P2A+9M7zzYtBrarqh5TVeuraj3wBuBPqur8iVS4Mot9XnPN/z+SHA38DL3AnzmrsguvJO8BTgKOSPIN4PeBs4E3Nh/MLcDzW6s8FvhGVa3qX6rltKuqvpfkL4DP0ztJ/aGq+qeJFL6EZX5ejwVeneQ24Hbg16pq4cn8VWEF/w+nwjLb9Wjg5Ul+SO9Cya9X1f92X/X4+WyyJLFKb7qWpK4ZhpKEYShJgGEoSYBhKEmAYShJgGGoRST5ySTvTvLVJF9ouhI7bYl11re7hVrm9s5KclRr/C1JHjTiuicluWQl2x3x/Y9K8v5meEOSJ6/gPV7V7tZLq49hqB+RJMDFwKeq6n5VtZHe44H3HuNmzwLuDMOq+pWqWhVP3FTVt6rqjGZ0A70urjRjDEMNcjJwa1Xd2fFqVV1fVX8Fd+4BXpFkS/M6ceEbDFsmycuazmqvTvJnSc4AHgq8q+lE9MAklyd5aLP8k5r3uDrJx0dtRJIzm+1cm+TPW9N3JPnj5v0+k+Qnm+k/3Yx/Psmrk+xoteXaJHcHXg08q6nzWQv3+Jrl1jfDv53kP5N8jN5jbLS2c2mzx31FkmNGbZPGaNIPR/tafS/gN4DXD5m/lqbnEnrPFs83w+tpHv4fsswpwKeBtc34uubn5bQ6AOiPA0fS69lmrr38gnpOAi5ZMO0oel1SHUnvsdPLgFObeQU8rRl+DfA7zfAlwJnN8K8BOwa06yzg/NZ2XsXu3XVd2yy/Ebim+Xc4FPiv/nLAx4H7N8MPBy6b9Gfuq1bns8laXZK8id4zqrdW1QnA3YDz0/t6hduBBwxYbbFlHg+8rap2AdTSzyU/gt7h+tYRl+87gV5vzd9t2vAues9FXwzcSi/4AL4APKEZfiRwajP8bnrdcq3UY4B/6LczyQebnwcDJwIX9M5GAL1uzTRhhqEGuQ44vT9SVeckOQKYbyb9JvAd4Dh6p1puGfAeiy0Tdu8teSnLXb693mJ+WM1uGb2g3pPfg9vY/XTTmtbwoLr3A75fVRv2YJsaA88ZapDLgDVJXtCatrY1fBjw7aq6A3gusP+A91hsmY8Cz0uyFiDJumb6TcAhA97nSuDnk8wtWH4pn23WO6LpMuxMej1rD/MZ7vojsGmRZRbW+T/0vk+EJMcDc830TwGnNec/DwGeBlBVNwJbkzyjWSdJjhuxTRojw1A/otlrOpVemGxN8jlgM/BbzSJ/DfxSks/QO/zdOeBtBi5TVZcCHwTmk1xFr0dvgLcDb+5fQGnV8l163UldlORqel9CNMjjknyj/6J33u4V9DojvRrYUlUfWKLp5wIvbtr7U/Q6PF3oE8CD+hdQ6H0B1LqmLS8AvtzUvaWp9apmmSta7/Ec4Jeb9lwHPH2JutQBu/CSGs3e6s1VVUk20buYYlDtIzxnKN1lI72LPqHXdf/zJluOuuSeoSThOUNJAgxDSQIMQ0kCDENJAgxDSQLg/wFTUfM6WEUWJwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "significance_mask_inv = ~significance_mask\n", "significance_mask_inv.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mask modifications\n", "\n", "### Mask dilation and erosion\n", "\n", "One can reduce or extend a mask using `binary_erode` and `binary_dilate` methods, respectively." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.064969Z", "iopub.status.busy": "2021-11-22T21:09:30.064653Z", "iopub.status.idle": "2021-11-22T21:09:30.121819Z", "shell.execute_reply": "2021-11-22T21:09:30.121999Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEHCAYAAAA0+iR9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS6ElEQVR4nO3de7SldV3H8fcHUIfhVgNoosIcWxpawciAF1RkYaWgLmGBOeSqWCpewlrkMlOzUrSLKF4KW7RyqZO3CiEyUrQkLiVejhMIZKUxYWZLbEaFmQER+PbHfjY8c9hnn33OzH72OZv3a61nnee+v7/Zcz7nueznt1NVSNID3R6TLkCSlgPDUJJY4WGYZPukaxiXaW2b7Vp5prltbSs6DCVpdzEMJQnISrubnGRlFSxp2VizZk1t2bJl4EGgR4aSHjBmZmYy3zLDUJIwDCUJMAwlCTAMJQkwDCUJMAwlCZhAGCY5JsndSU5rzduQZFOSs7uuR5Kg4zBMsifwNuDTcxZtAI4Bnpxk3y5rkiTo/sjwV4GLgFvmzO9/ELJa45LUmc7CMMkjgFOACwYsvhiYBWar6rauapKkvr06fK13A79ZVXcnOx/8VdVGYGOHtUjSTsbaUUOSs4Azm8kDuO8U+CBgB/CyqrpkhP20+1NbvTtrlPTAsWbNGrZu3bqjP11V+/THJ9JrTZIPApdW1ceXsK291khakvXr1zM7OzvwvoSfM5Qkur1meK+qOmMSrytJ8/HIUJIwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQnoMAyTHJ/k+0mubYbfaS3bkGRTkrO7qkeS2vbq+PWurqrnDpi/ATgG+EiSfatqW8d1SXqAWy6nyWl+VmtckjrTdRg+Jcl1ST6V5Cdb8y8GZoHZqrqt45okiVRVNy+U7A/cU1XbkpwEvKeqHjPitttbk6vHUqCkqbdmzRq2bt26oz9dVfv0x0c6MkxyWJKfacb3TrLfiNud1b9hAtx7LbCqPgk8KMlBo+ynqvbpD6OsL0mDzMzMzJsnC4ZhkjOBjwN/2sx6JHDJKC9cVe+tqnVVtQ64J0mafT6xee0tI7dCksZolLvJZwFPBL4AUFVfS/LQJbzWacArk9wF3A5sqK7O0SVpAaOE4Q+q6s7moI4ke9G767soVXU+cP5it5OkLoxyzfDKJG8A9k7ys8CFwN+OtyxJ6tYoYfg64DvA9cDLgU8CbxxnUZLUtQVPk6vqHuDPmkGSptK8YZjkeoZcG6yqI8ZSkSRNwLAjw/4zxGc1Pz/U/HwRsOP+q0u7blc/YNC/0Sct1rxhWFU3AyR5alU9tbXodUn+GThn3MVJUldGuYGyT5Kn9SeSHAv4JIikqTLK5wxfArw/yQHN9PeAF4+tIkmagFHuJn8ZOLLpaCFV9f3xlyVJ3VowDNs9UjfTAFSV1wwlTY1RTpPb3WetoneX+avjKUcPFON6LH0x+/XOs9pGOU0+rz2d5B3AJ8ZWkSRNwFJ6ul4NPHp3FyJJkzTKNcP2kyh7AgcDbxlnUZLUtVGuGba/ze4u4NtVddeY6pGkiRjlNPmtVXVzM/xPVd2V5EMLbyYtb1U18qDpN0oYtr/Frt+56/rxlCNJkzFvGCZ5fZLbgCOS3NoMtwHfBv6mswolqQPzhmFV/UFV7Qe8var2b4b9qurAqnp9hzVK0tgN68/w8Kr6N+DCJEfNXV5Vm8ZamSR1aNjd5FcDLwPOG7CsgBPGUpEkTcCw/gxf1oyeWFV3tJclWTXWqjRVpuFu7KA2+DjfdBnlbvLnRpwnSSvWsGuGPwY8gt5XhD4B6P8Z3J/eI3mSNDWGXTN8FnAG8Ejgna35twFvGGNNktS5YdcMNwIbk5xaVRd1WJMkdW6ULrwuSvIcek+irGrNt3NXSVNjwRsoSS4AXgj8Kr3rhi8ADhtzXZLUqVHuJh9bVb8EfLeq3gw8BXjUeMuSpG6NEoa3Nz93JDkE+CEwM76SJKl7o/RneGmSHwHeDmyi9/TJn42zKEnq2ig3UPq9Wl+U5FJ6N1EOH2tVktSxRX0HSlX9oPne5AvHVI9WMDtG1Uq2lC+EgvueRpGkqbDUMPTPvaSpMuzZ5L9lcOgFOHBsFUnSBAy7gfKOJS6TpBVn2LPJV3ZZiFYWb4xo2iz1mqEkTRXDUJIwDCUJGK3Xmr9vHsfrT/9okk+PtSpJ6tgozyYfVFXf609U1XeTPHR8JWklGPRlSN5U0Uo2ymnyPUkO7U8kOQw/dC1pyoxyZPhbwD8l6X/U5jh636csSVNjwSPDqroMOAr4S+CvgPVVtaRrhkmOT3Jtkhtb4UqSDUk2JTl7KfuVpF017HG8w6vq35Ic1cz6VvPz0CSHVtWmxbxQcxPmT4BnV9U35lx33AAcA3wkyb5VtW0x+5akXTXsNPnV9E6HzxuwrIATFvlavwBcXFXfAKiqW1rL+lfjC3vEkTQBwx7H618XPLGq7mgvS7JqwCYLeSzwoCRXAPsB76mqP2+WXQzMAh+uqtuWsG91zDvHmjZZ6D91kk1VddRC8xZ8oeR84GjgmcDewDXAc6rqP0bYdntrcvViXlfjYRgO/niRlrc1a9awdevWHf3pqtqnPz7smuGPAY8A9k7yBO47fd2fEQMpyVnAmc3kXwGXVdV2YHuSq4AjgQXDsF1wEn8LJS3JzMwMW7Zs2WfQsmHXDJ8FnAE8kt51w34Y3gq8YZQXrqr3Au8FSPI44PwkewEPBp4EvGuU/UjSuA27ZrgR2Jjk1Kq6aFdfqKq+muQy4CvAPcD7quqGXd2vJO0OozyBsn7As8lvXcqLVdXbq+rxVfVTVfXupexDy0OS+w3SSjZKGJ4499lk4KSxVSRJEzBKGO6Z5CH9iSR7Aw8Zsr4krTijPJv8YeCzST5A70PRLwY2jrUqSerYgmFYVecmuZ7e5wMDvGWpzyZL0nI1ypEhVfUp4FNjrkWSJmaUnq6fnORLSbYluTPJ3Ulu7aI4rSyD7jDPNyxX09AGLc0oN1DOB04HvkbvMbqXAn88zqIkqWujniZ/PcmeVXU38IEknxtzXZLUqVHCcEeSBwPXJjkX+F9g4LN9krRSjXKa/IvAnsCrgO3Ao4BTx1mUJHVtlI/W3NyM3g68ebzl6IFivhsQ4+gazJsdGsWwLryuZ8i34FXVEWOpSJImYNiR4XM7q0KSJmxYF143z7dMkqaNH7qWJPzQtSQBfuhay4x3fjUpfuhakhj9Q9d74IeuJU2xxXzo+g780LWkKTXvkWGS5zffe9yf/kKSm5rhtG7Kk6RuDDtNfi3widb0Q4BjgOOBV46xJknq3LDT5AdX1X+3pv+pqrYAW5J4A0XSVBl2ZPij7YmqelVr8uDxlCNJkzEsDL+Q5My5M5O8HPji+EqSpO4NO03+deCSJL8AbGrmrad37fDkMdclSZ0a1lHDLcCxSU4AfrKZ/XdVdXknlUlSh0b5nOHlgAEoaaqN8gSKJE09w1CSMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJ6DAMk/xGkmub4YYkdydZ0yzbkGRTkrO7qkeS2joLw6p6e1Wtq6p1wOuBK6tqa7N4A70vqH9ykn27qkmS+iZ1mnw68LHWdJqf1RqXpM50HoZJVgPPBi5qzb4YmAVmq+q2rmuSpAW/HW8Mngf8c+sUmaraCGycQC2SBIz5yDDJWa2bJoc0szew8ynyKPvZ3h92f5WSHig2b948b56kqjorJMkBwGbgUVW1pGBL0l3BkqbK+vXrmZ2dHXhfoutrhqcAn1lqEErSuHQahlX1wara0OVrStIofAJFkjAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQmAvSZdwGIdeOCBrF279t7pzZs3MzMzM7mCxsR2rSy2a2XYvHlzzbcsVfMuWxGSbK+qfSZdx+5mu1YW27XyeZosSRiGkgRMwWmyJO0OHhlKEoahJAGGoSQByzQMk7w/yS1JbmjNW5fk80muTTKb5InN/Bc18/rDPUnWNcuOb9Y9d0JN2cli2tUsOyLJNUluTHJ9klXN/BXbriRrk9zeer8uaG2zYtvVWn5okm1JXtOat2LbleSJrffquiSntLZZVu3aZVW17AbgOOAo4IbWvM8AJzbjJwFXDNjup4GbWtN/CewNnAccvpLaRe8D8V8BjmymDwT2nIJ2rW2vN2c/K7ZdreUXARcCr5mGdgGrgb2a8YcDt7Sml1W7dnVYlkeGVXUVsHXubGD/ZvwA4FsDNj0d+Fhreo9mu3uA7OYyF22R7fo54CtVdV2z7ZaqurtZtpLbNcyKbleSk4GbgBvnbLNi21VVO6rqrmb+qma9vmXVrl026TQe8tdrLTv/5Xoc8A3gv4H/AQ4bsM1/Aj/Vmn4WsAk4b9LtWWy7gLOBDwGfbtrw2ilp11pgO/AvwJXA06ekXfsA1wD7Am9i5yPDFduuZtmT6AX8NuCU5dyuXfo3mXQBi3iz/gg4tRn/eeAf5qz/JOD6Sde9u9oFvAbYDBxE71TlGuCZk65/N7TrIcCBzfj65pdv/0nXvxva9Q7g55vxncJwOQ6L/f1q5j8O+CKwatL1j+XfZNIFLOLN+j73fUg8wK1z1n8X8IZJ17272gVsAD7YWu+3gd+YdP276/1qrXcFcPSk698N79fVwH81w/fonYa+atL1j+H9+sfl/H7tyrAsrxnO41vAM5rxE4Cv9Rck2QN4AfAXE6hrV83Xrk8DRyRZnWSvZp1/nUB9SzWwXUkOTrJnM/5o4DH0rrOtFAPbVVVPr6q1VbUWeDfw+1V1/kQqXJr53q+Z5v8fSQ4DfoJe4E+dZdmFV5KPAccDByX5JvC7wJnAe5o35g7gZa1NjgO+WVXL+pdqMe2qqu8meSfwJXoXqT9ZVX83kcIXsMj36zjgnCR3AXcDr6iquRfzl4Ul/D9cERbZrqcBr0vyQ3o3Sn6lqv6v+6rHz2eTJYll+qFrSeqaYShJGIaSBBiGkgQYhpIEGIaSBBiGmkeShyX5aJKbkny56UrslAW2WdvuFmqRr3dGkkNa0+9L8vgRtz0+yaVLed0R939Iko834+uSnLSEfbyp3a2Xlh/DUPeTJMAlwFVV9eiqWk/v8cBHjvFlzwDuDcOqemlVLYsnbqrqW1V1WjO5jl4XV5oyhqEGOQG4s6ru7Xi1qm6uqj+Ge48Ar06yqRmOnbuDYeskeW3TWe11Sf4wyWnA0cBHmk5E905yRZKjm/Wf3ezjuiSfHbURSU5vXueGJG9rzd+W5Pea/X0+ycOa+T/eTH8pyTlJtrXackOSBwPnAC9s6nzh3CO+Zr21zfhvJfn3JP9A7zE2Wq9zWXPEfXWSw0dtk8Zo0g9HOyy/Afg14F1Dlq+m6bmE3rPFs834WpqH/4escyLwOWB1M72m+XkFrQ4A+tPAwfR6tplprz+nnuOBS+fMO4Rel1QH03vs9HLg5GZZAc9rxs8F3tiMXwqc3oy/Atg2oF1nAOe3XudN7Nxd1w3N+uuB65t/h/2Br/fXAz4LPKYZfxJw+aTfc4dans8ma3lJ8l56z6jeWVXHAA8Czk/v6xXuBh47YLP51vkZ4ANVtQOgFn4u+cn0Ttc3j7h+3zH0emv+TtOGj9B7LvoS4E56wQfwZeBnm/GnACc34x+l1y3XUj0d+Ot+O5N8ovm5L3AscGHvagTQ69ZME2YYapAbgVP7E1V1VpKDgNlm1q8D3waOpHep5Y4B+5hvnbBzb8kLWez67e3m88NqDsvoBfWu/B7cxc6Xm1a1xgfVvQfwvapatwuvqTHwmqEGuRxYleSVrXmrW+MHAP9bVfcAvwjsOWAf863zGeDFSVYDJFnTzL8N2G/Afq4BnpFkZs76C/lCs91BTZdhp9PrWXuYz3PfH4EN86wzt87/ovd9IiQ5Cphp5l8FnNJc/9wPeB5AVd0KbE7ygmabJDlyxDZpjAxD3U9z1HQyvTDZnOSLwEbgN5tV/gT45SSfp3f6u33AbgauU1WXAZ8AZpNcS69Hb4APAhf0b6C0avkOve6kLk5yHb0vIRrkmUm+2R/oXbd7Pb3OSK8DNlXV3yzQ9LOBVzftfTi9Dk/n+kfg8f0bKPS+AGpN05ZXAv/R1L2pqfXaZp2rW/t4EfCSpj03As9foC51wC68pEZztHp7VVWSDfRuphhUDxBeM5Tus57eTZ/Q67r/xZMtR13yyFCS8JqhJAGGoSQBhqEkAYahJAGGoSQB8P9L/w36fmvs7wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask = significance_mask_inv.binary_erode(width=0.2 * u.deg, kernel=\"disk\")\n", "mask.plot();" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.124734Z", "iopub.status.busy": "2021-11-22T21:09:30.124422Z", "iopub.status.idle": "2021-11-22T21:09:30.180740Z", "shell.execute_reply": "2021-11-22T21:09:30.180918Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEHCAYAAAA0+iR9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATI0lEQVR4nO3deZBlZXnH8e8PUIdhSwZwwYVpUxrUBCaMuJCoFGoUlxIK1EErCaWiMWgKLTVqNmPMIopLxJQpLXXiFkUIKiqaSFBUXMYJCESjxhHXEjOjwsxAEHjyxz0XzrS3b9/u6Xvu7Z7vp+pUn/0879zpp8953/O+N1WFJO3p9pp0AJI0DUyGksQyT4ZJdkw6hnFZqWWzXMvPSi5b27JOhpK0VEyGkgRkubUmJ1leAUuaGmvWrKmtW7cOvAn0zlDSHmNmZiZzbTMZShImQ0kCTIaSBJgMJQkwGUoSYDKUJGACyTDJMUluSXJKa92GJJuTnNl1PJIEHSfDJHsDrwY+MWvTBuAY4CFJ9u8yJkmC7u8Mnw+cB1w7a33/RchqzUtSZzpLhknuDpwEvGXA5vOBTcCmqrq+q5gkqW+fDq/1BuBPquqWZNebv6raCGzsMBZJ2sVYB2pIcgZwerN4ELc/Ah8C7ASeXVUXjHCe9nhqq5cyRkl7jjVr1rBt27ad/eWq2q8/P5FRa5K8E7iwqj64iGMdtUbSoqxfv55NmzYNbJfwPUNJots6w9tU1WmTuK4kzcU7Q0nCZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgDYZ9IBaOWrqs6ulaSza2ll8c5QkjAZShJgMpQkwGQoSYANKFqkLhtFFmIhcdnYojbvDCUJk6EkASZDSQJMhpIEmAwlCbA1WXswW57V5p2hJGEylCTAZChJgMlQkgAbUDSPae121zUbW1Y+7wwlCZOhJAEdJsMkxyX5eZLLm+kvWts2JNmc5Myu4pGktq7rDC+tqicMWL8BOAZ4T5L9q2p7x3FJ2sNNSwNKv8a5WvPqmI0lS2PQv6ONKtOv6zrDhya5IsnHkzygtf58YBOwqaqu7zgmSSJd3Q0kORC4taq2J3kc8Maqus+Ix+5oLa4eS4DyznCMvDOcDmvWrGHbtm07+8tVtV9/fqQ7wySHJ3lUM79vkgNGPO6MfoMJcFtdYFV9DLhDkkNGOU9V7defRtlfkgaZmZmZM5/MmwyTnA58EPinZtU9gAtGuXBVvbmq1lXVOuDWNH8ekzyoufbWkUshSWM0SgPKGcCDgC8CVNU3k9x5Edc6BXhukpuBG4AN5XOZpCkxSjL8v6q6qV/nkWQfeq2+C1JV5wDnLPQ4aSWY6+++dYnTY5Q6w08neTmwb5JHA+cCHxlvWJLUrVGS4UuBnwBXAs8BPgb82TiDkqSuzfuYXFW3Am9tJklakeZMhkmuZEjdYFUdOZaIJGkCht0Z9vsQn9H8fFfz8+nAzl/eXcuJDfnTwa5702POZFhV1wAk+e2q+u3Wppcm+RzwynEHJ0ldGaUBZb8kv9NfSHIsYE8QSSvKKO8ZPhN4e5KDmuWfAc8YW0SSNAGjtCZ/BTiqGWghVfXz8YclSd2aNxm2R6RulgGoKusMJa0Yozwmt4fPWkWvlflr4wlHki3MkzHKY/LZ7eUkrwU+PLaIJGkCFjPS9Wrg3ksdiCRN0ih1hu2eKHsDhwJ/Pc6gJKlro9QZtr/N7mbgx1V185jikaSJGOUx+VVVdU0z/aCqbk7yrvkPk6TlY5Rk2P4Wu/7gruvHE44kTcacyTDJy5JcDxyZ5Lpmuh74MfChziKUpA7MmQyr6u+q6gDgNVV1YDMdUFUHV9XLOoxRksZu2HiGR1TV14Fzkxw9e3tVbR5rZJLUoWGtyS8Eng2cPWBbAcePJSJJmoBh4xk+u5k9oapubG9LsmqsUWnJOIirNJpRWpM/P+I6SVq2htUZ3hW4O72vCP0toN9T/EB6XfIkacUYVmf4GOA04B7A61rrrwdePsaYJKlzw+oMNwIbk5xcVed1GJMkdW6UIbzOS/J4ej1RVrXWO7jrlLGxRFq8eRtQkrwFeCrwfHr1hk8GDh9zXJLUqVFak4+tqt8HflpVfwU8FLjneMOSpG6NkgxvaH7uTHIY8AtgZnwhSVL3RhnP8MIkvwK8BthMr/fJW8cZlCR1bZQGlP6o1ucluZBeI8oRY41Kkjq2oO9Aqar/a743+dwxxSNJE7GYL4SC23ujSNKKsNhk6AttklaUYX2TP8LgpBfg4LFFJEkTMKwB5bWL3CZJy86wvsmf7jIQSZqkxdYZStKKYjKUJEyGkgSMNmrNvzXd8frLv5rkE2ONSpI6Nsqd4SFV9bP+QlX9FLjz2CKSpAkYJRnemuRe/YUkh+NL15JWmFFGrflT4LNJ+q/aPJze9ylL0oox751hVV0EHA28H/gAsL6qFlVnmOS4JJcnubqVXEmyIcnmJGcu5ryStLuGdcc7oqq+nuToZtUPm5/3SnKvqtq8kAs1jTD/CDy2qr6bpF3vuAE4BnhPkv2ravtCzi1Ju2vYY/IL6T0Onz1gWwHHL/BaTwPOr6rvAlTVta1t/VFwCkfEkTQBw7rj9esFT6iqG9vbkqwacMh87gvcIcklwAHAG6vqn5tt5wObgHdX1fWLOLck7ZbM9/WSSTZX1dHzrZv3Qsk5wAOBRwL7ApcBj6+qb4xw7I7W4uqFXHdP4leFrlyJD0xLYc2aNWzbtm1nf7mq9uvPD6szvCtwd2DfJL/F7Y+vBzJiQkpyBnB6s/gB4KKq2gHsSPIZ4Chg3mTYDjiJv/GSFmVmZoatW7fuN2jbsDrDxwCnAfegV2/YT4bXAS8f5cJV9WbgzQBJ7geck2Qf4I7Ag4HXj3IeSRq3YXWGG4GNSU6uqvN290JV9bUkFwFfBW4F3lZVV+3ueSVpKYzSA2X9gL7Jr1rMxarqNVV1/6r6jap6w2LOIUnjMEoyPGFA3+THjS0iSZqAUZLh3knu1F9Isi9wpyH7S9KyM0rf5HcDn0ryDnovRT8D2DjWqCSpY/Mmw6o6K8mV9N4PDPDXi+2bLEnTapQ7Q6rq48DHxxyLJE3MKCNdPyTJl5NsT3JTkluSXNdFcJLUlVEaUM4BTgW+Sa8b3bOAN40zKEnq2qiPyd9KsndV3QK8I8nnxxyXJHVqlGS4M8kdgcuTnAX8CBjYt0+SlqtRHpN/D9gbeB6wA7gncPI4g5Kkro3yas01zewNwF+NNxztjkHDPDmslzSaYUN4XcmQb8GrqiPHEpEkTcCwO8MndBaFJE3YsCG8rplrmyStNL50LUn40rUkAb50veLN9UVCtjJLu/Kla0li9Jeu98KXriWtYAt56fpGfOla0go1551hkic133vcX/5ikm830yndhCdJ3Rj2mPwS4MOt5TsBxwDHAc8dY0zqQJJfmjQd/GwmY9hj8h2r6nut5c9W1VZgaxIbUCStKMPuDH+1vVBVz2stHjqecCRpMoYlwy8mOX32yiTPAb40vpAkqXvDHpNfAFyQ5GnA5mbdenp1hyeOOS5J6tSwgRquBY5NcjzwgGb1R6vq4k4ik6QOjfKe4cWACXAP4OCw42OL8PQbpQeKJK14JkNJwmQoSYDJUJKAEccz1J7L8RDnZqPIyuKdoSRhMpQkwGQoSYDJUJIAk6EkAbYma5FWatc9W4j3XN4ZShImQ0kCTIaSBJgMJQmwAUVLaFq77tkoolF4ZyhJdJgMk7w4yeXNdFWSW5KsabZtSLI5yZldxSNJbZ0lw6p6TVWtq6p1wMuAT1fVtmbzBnpfUP+QJPt3FZMk9U3qMflU4H2t5X6lTrXmJakznTegJFkNPBZofyn9+cAm4N1VdX3XMWm8bMDQcjCJ1uQnAp9rPSJTVRuBjROIRZKAMT8mJzmj1WhyWLN6A7s+Io9ynh39aemjlLSn2LJly5z5JF2+A5bkIGALcM+qWlRiS7L8RwOQNBHr169n06ZNA+ttum5AOQn45GIToSSNS6fJsKreWVUburymJI3CHiiShMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJMBlKEmAylCTAZChJgMlQkgCToSQBJkNJAkyGkgSYDCUJgH0mHcBCHXzwwaxdu/a25S1btjAzMzO5gMbEci0vlmt52LJlS821LVVzblsWkuyoqv0mHcdSs1zLi+Va/nxMliRMhpIErIDHZElaCt4ZShImQ0kCTIaSBExpMkzy9iTXJrmqtW5dki8kuTzJpiQPatY/vVnXn25Nsq7Zdlyz71kTKsouFlKuZtuRSS5LcnWSK5OsatYv23IlWZvkhtbn9ZbWMcu2XK3t90qyPcmLWuuWbbmSPKj1WV2R5KTWMVNVrt1WVVM3AQ8Hjgauaq37JHBCM/844JIBx/0m8O3W8vuBfYGzgSOWU7novRD/VeCoZvlgYO8VUK617f1mnWfZlqu1/TzgXOBFK6FcwGpgn2b+bsC1reWpKtfuTlN5Z1hVnwG2zV4NHNjMHwT8cMChpwLvay3v1Rx3K5AlDnPBFliu3wW+WlVXNMdurapbmm3LuVzDLOtyJTkR+DZw9axjlm25qmpnVd3crF/V7Nc3VeXabZPOxkP+eq1l179c9wO+C3wP+AFw+IBj/gf4jdbyY4DNwNmTLs9CywWcCbwL+ERThpeskHKtBXYA/wl8GnjYCinXfsBlwP7AK9j1znDZlqvZ9mB6CX47cNI0l2u3/k0mHcACPqx/AE5u5p8C/Pus/R8MXDnpuJeqXMCLgC3AIfQeVS4DHjnp+JegXHcCDm7m1ze/fAdOOv4lKNdrgac087skw2mcFvr71ay/H/AlYNWk4x/Lv8mkA1jAh/Vzbn9JPMB1s/Z/PfDySce9VOUCNgDvbO3358CLJx3/Un1erf0uAR446fiX4PO6FPhOM/2M3mPo8yYd/xg+r/+Y5s9rd6aprDOcww+BRzTzxwPf7G9IshfwZOBfJhDX7pqrXJ8AjkyyOsk+zT7/NYH4FmtguZIcmmTvZv7ewH3o1bMtFwPLVVUPq6q1VbUWeAPwt1V1zkQiXJy5Pq+Z5v8fSQ4Hfp1ewl9xpnIIryTvA44DDknyfeAvgdOBNzYfzI3As1uHPBz4flVN9S/VQspVVT9N8jrgy/QqqT9WVR+dSODzWODn9XDglUluBm4B/rCqZlfmT4VF/D9cFhZYrt8BXprkF/QaSv6oqv63+6jHz77JksSUvnQtSV0zGUoSJkNJAkyGkgSYDCUJMBlKEmAy1ByS3CXJe5N8O8lXmqHETprnmLXtYaEWeL3TkhzWWn5bkvuPeOxxSS5czHVHPP9hST7YzK9L8rhFnOMV7WG9NH1MhvolSQJcAHymqu5dVevpdQ+8xxgvexpwWzKsqmdV1VT0uKmqH1bVKc3iOnpDXGmFMRlqkOOBm6rqtoFXq+qaqnoT3HYHeGmSzc107OwTDNsnyUuawWqvSPL3SU4BHgi8pxlEdN8klyR5YLP/Y5tzXJHkU6MWIsmpzXWuSvLq1vrtSf6mOd8XktylWf9rzfKXk7wyyfZWWa5KckfglcBTmzifOvuOr9lvbTP/p0n+O8m/0+vGRus6FzV33JcmOWLUMmmMJt052mn6JuCPgdcP2b6aZuQSen2LNzXza2k6/w/Z5wTg88DqZnlN8/MSWgMA9JeBQ+mNbDPT3n9WPMcBF85adxi9IakOpdft9GLgxGZbAU9s5s8C/qyZvxA4tZn/Q2D7gHKdBpzTus4r2HW4rqua/dcDVzb/DgcC3+rvB3wKuE8z/2Dg4kl/5k41nX2TNV2SvJleH9WbquoY4A7AOel9vcItwH0HHDbXPo8C3lFVOwFq/n7JD6H3uL5lxP37jqE3WvNPmjK8h16/6AuAm+glPoCvAI9u5h8KnNjMv5fesFyL9TDgX/vlTPLh5uf+wLHAub3aCKA3rJkmzGSoQa4GTu4vVNUZSQ4BNjWrXgD8GDiKXlXLjQPOMdc+YdfRkuez0P3bx83lF9XcltFL1Lvze3Azu1Y3rWrND4p7L+BnVbVuN66pMbDOUINcDKxK8tzWutWt+YOAH1XVrcDvAXsPOMdc+3wSeEaS1QBJ1jTrrwcOGHCey4BHJJmZtf98vtgcd0gzZNip9EbWHuYL3P5HYMMc+8yO8zv0vk+EJEcDM836zwAnNfWfBwBPBKiq64AtSZ7cHJMkR41YJo2RyVC/pLlrOpFeMtmS5EvARuBPml3+EfiDJF+g9/i7Y8BpBu5TVRcBHwY2Jbmc3ojeAO8E3tJvQGnF8hN6w0mdn+QKel9CNMgjk3y/P9Grt3sZvcFIrwA2V9WH5in6mcALm/Lejd6Ap7P9B3D/fgMKvS+AWtOU5bnAN5q4NzexXt7sc2nrHE8HntmU52rgSfPEpQ44hJfUaO5Wb6iqSrKBXmOKiWoPYZ2hdLv19Bp9Qm/o/mdMNhx1yTtDScI6Q0kCTIaSBJgMJQkwGUoSYDKUJAD+H03D/pJDoBeyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mask = significance_mask_inv.binary_dilate(width=0.2 * u.deg)\n", "mask.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Boundary mask\n", "\n", "In the following example we use the Fermi dataset previously loaded and add its `mask_fit` taking into account a margin based on the psf width. The margin width is determined using the `containment_radius` method of the psf object and the mask is created using the `boundary_mask` method available on the geometry object." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.184475Z", "iopub.status.busy": "2021-11-22T21:09:30.184155Z", "iopub.status.idle": "2021-11-22T21:09:30.251467Z", "shell.execute_reply": "2021-11-22T21:09:30.251654Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEHCAYAAAA0+iR9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAASu0lEQVR4nO3deZBlZX2H8ecLqMOwxQFccGHalAaXwMiIIm4UahTUEgrUQcuEUnEJmkLLDTWJWxZRXCJYpLSUiaJRhKCiolFEURFtJyAQjRpHXLDEzIgwMyACv/xxz4Uzbfft2z1zz+1unk/VrT7rPb+37/R3zvreVBWSdEe3w7gLkKSFwDCUJBZ5GCbZPO4aRmWpts12LT5LuW1tizoMJWl7MQwlCchiu5qcZHEVLGnBWLFiRW3YsGHanUD3DCXdYUxMTGSmeYahJGEYShJgGEoSYBhKEmAYShJgGEoSMIYwTHJQkluSHNOatibJuiQndl2PJEHHYZhkR+BtwBemzFoDHAQcnGTXLmuSJOh+z/BlwNnANVOm92+ErNawJHWmszBMci/gKOD0aWafA0wCk1V1fVc1SVLfTh1u693Aa6rqlmTrnb+qWgus7bAWSdrKSDtqSHICcHwzuge3HwLvBWwBXlhV5w7xPu3+1JZvzxol3XGsWLGCjRs3bumPV9Uu/eGx9FqT5AzgvKr65DzWtdcaSfOyevVqJicnp70u4X2GkkS35wxvU1XHjWO7kjQT9wwlCcNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSAMNQkgDDUJIAw1CSgDF17rrYjeOrEqTFbuoXwS007hlKEoahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAEdhmGSQ5P8LsmlzevvWvPWJFmX5MSu6pGktp063t5FVfXUaaavAQ4Czkyya1Vt6rguSXdwC+UwOc3Pag1LUme6DsNHJrksyeeTPLg1/RxgEpisqus7rkmSSFV1s6Fkd+DWqtqU5AjgPVV1/yHX3dwaXT6SAuegq9+ZtJQk4z/oW7FiBRs3btzSH6+qXfrDQ+0ZJtk3yROa4Z2T7Dbkeif0L5gAt50LrKrPAXdKstcw71NVu/RfwywvSdOZmJiYMU9mDcMkxwOfBP61mXRv4NxhNlxVp1XVqqpaBdya5r+GJA9vtr1h6FZI0ggNczX5BODhwCUAVfWjJHebx7aOAV6S5GbgBmBNebwpaYEYJgx/X1U39Y/3k+xE76rvnFTVqcCpc11PkrowzDnDryZ5HbBzkicCZwGfGW1ZktStYcLwtcBvgMuBFwGfA94wyqIkqWuzHiZX1a3A+5uXJC1JM4ZhkssZcG6wqvYfSUWSNAaD9gz7zxCf0Pz8cPPzOcCWP15ckhavGcOwqq4CSPKoqnpUa9Zrk3wDePOoi5OkrgxzAWWXJI/ujyQ5BPBJEElLyjD3GT4f+GCSPZrxa4HnjawiSRqDYa4mfxc4oOloIVX1u9GXJUndmjUM2z1SN+MAVJXnDCUtGcMcJre7z1pG7yrz90dTjiSNxzCHyae0x5O8A/j0yCqSpDGYT0/Xy4H7be9CJGmchjln2H4SZUdgb+AtoyxKkro2zDnD9rfZ3Qz8uqpuHlE9kjQWwxwmv7Wqrmpev6yqm5N8ePbVJGnxGCYM299i1+/cdfVoypGk8ZgxDJOclOR6YP8k1zWv64FfA5/qrEJJ6sCMYVhV/1RVuwFvr6rdm9duVbVnVZ3UYY2SNHKD+jPcr6p+AJyV5MCp86tq3Ugrk6QODbqa/ArghcAp08wr4LCRVCRJYzCoP8MXNoOHV9WN7XlJlo20Kknq2DBXk7855DRJWrQGnTO8B3Avel8R+lAgzazd6T2SJ0lLxqBzhk8CjgPuDbyzNf164HUjrEmSOjfonOFaYG2So6vq7A5rkqTODdOF19lJnkLvSZRlrel27ippyZj1AkqS04FnAS+jd97wGcC+I65Lkjo1zNXkQ6rqL4HfVtWbgEcC9xltWZLUrWHC8Ibm55Yk+wB/ACZGV5IkdW+Y/gzPS/InwNuBdfSePnn/KIuSpK4NcwGl36v12UnOo3cRZb+RViVJHZvTd6BU1e+b700+a0T1SNJYzOcLoeD2p1EkaUmYbxjW7ItI0uIx6NnkzzB96AXYc2QVSdIYDLqA8o55zpOkRWfQs8lf7bIQSRqn+Z4zlKQlxTCUJAxDSQKG67XmP5vH8frjd03yhZFWJUkdG2bPcK+qurY/UlW/Be42sookaQyGCcNbk9y3P5JkX7zpWtISM0yvNa8Hvp6kf6vNY+l9n7IkLRmz7hlW1fnAgcDHgU8Aq6tqXucMkxya5NIkV7bClSRrkqxLcuJ83leSttWgx/H2q6ofJDmwmXR18/O+Se5bVevmsqHmIsz7gCdX1c+StM87rgEOAs5MsmtVbZrLe0vSthp0mPwKeofDp0wzr4DD5ritZwPnVNXPAKrqmta8fi84hT3iSBqDQY/j9c8LHl5VN7bnJVk2zSqzeQBwpyQXArsB76mqf2vmnQNMAh+pquvn8d6StE1SNfjCcJJ1VXXgbNNm3VByKvAw4PHAzsDFwFOq6odDrLu5Nbp8Ltsdhdl+Z5L+WDL+g74VK1awcePGLf3xqtqlPzzonOE9gHsBOyd5KLcfvu7OkIGU5ATg+Gb0E8D5VbUZ2Jzka8ABwKxh2C44iUkkaV4mJibYsGHDLtPNG3TO8EnAccC96Z037IfhdcDrhtlwVZ0GnAaQ5IHAqUl2Au4MPAJ41zDvI0mjNuic4VpgbZKjq+rsbd1QVX0/yfnA94BbgQ9U1RXb+r6StD0M8wTK6mmeTX7rfDZWVW+vqgdV1UOq6t3zeQ9JGoVhwvDwaZ5NPmJkFUnSGAwThjsmuUt/JMnOwF0GLC9Ji84wzyZ/BPhykg/Ruyn6ecDakVYlSR2bNQyr6uQkl9O7PzDAW+b7bLIkLVTD7BlSVZ8HPj/iWiRpbIbp6frgJN9JsinJTUluSXJdF8VJUleGuYByKnAs8CN6j9G9AHjvKIuSpK4Ne5j84yQ7VtUtwIeSfHPEdUlSp4YJwy1J7gxcmuRk4FfAtM/2SdJiNcxh8nOBHYGXApuB+wBHj7IoSeraMLfWXNUM3gC8abTlSNJ4DOrC63IGfAteVe0/kookaQwG7Rk+tbMqJGnMBnXhddVM8yRpqfGma0nCm64lCfCma0kCvOlakoDhb7reAW+6lrSEzeWm6xvxpmtJS9SMe4ZJnt5873F//JIkP2lex3RTniR1Y9Bh8quBT7fG7wIcBBwKvGSENUlS5wYdJt+5qn7eGv96VW0ANiTxAoqkJWXQnuFd2yNV9dLW6N6jKUeSxmNQGF6S5PipE5O8CPj26EqSpO4NOkx+OXBukmcD65ppq+mdOzxyxHVJUqcGddRwDXBIksOABzeTP1tVF3RSmSR1aJj7DC8ADEBJS9owT6BI0pJnGEoShqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBhqEkAYahJAGGoSQBHYZhklclubR5XZHkliQrmnlrkqxLcmJX9UhSW2dhWFVvr6pVVbUKOAn4alVtbGavofcF9Qcn2bWrmiSpb1yHyccCH2uNp/lZrWFJ6kznYZhkOfBk4OzW5HOASWCyqq7vuiZJmvXb8UbgacA3WofIVNVaYO0YapEkYMR7hklOaF002aeZvIatD5GHeZ/N/df2r1LSHcX69etnzJNUVWeFJNkDWA/cp6rmFWxJuit4Bl3+zqSlIhn/5YDVq1czOTk5bSFdnzM8CvjifINQkkal0zCsqjOqak2X25SkYfgEiiRhGEoSYBhKEmAYShJgGEoSMJ4nUBa9hXC/lKTtyz1DScIwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCTAMJQkwDCUJMAwlCQAdhp3AXO15557snLlytvG169fz8TExPgKGhHbtbjYrsVh/fr1NdO8VM04b1FIsrmqdhl3Hdub7VpcbNfi52GyJGEYShKwBA6TJWl7cM9QkjAMJQkwDCUJWKBhmOSDSa5JckVr2qok30pyaZLJJA9vpj+nmdZ/3ZpkVTPv0GbZk8fUlK3MpV3NvP2TXJzkyiSXJ1nWTF+07UqyMskNrc/r9NY6i7Zdrfn3TbIpyStb0xZtu5I8vPVZXZbkqNY6C6pd26yqFtwLeCxwIHBFa9oXgcOb4SOAC6dZ78+Bn7TGPw7sDJwC7LeY2kXvhvjvAQc043sCOy6Bdq1sLzflfRZtu1rzzwbOAl65FNoFLAd2aobvCVzTGl9Q7drW14LcM6yqrwEbp04Gdm+G9wCunmbVY4GPtcZ3aNa7Fch2LnPO5tiuvwC+V1WXNetuqKpbmnmLuV2DLOp2JTkS+Alw5ZR1Fm27qmpLVd3cTF/WLNe3oNq1zcadxgP+91rJ1v9zPRD4GfBz4JfAvtOs87/AQ1rjTwLWAaeMuz1zbRdwIvBh4AtNG169RNq1EtgM/BfwVeAxS6RduwAXA7sCb2TrPcNF265m3iPoBfwm4KiF3K5t+p2Mu4A5fFj/AhzdDD8T+NKU5R8BXD7uurdXu4BXAuuBvegdqlwMPH7c9W+Hdt0F2LMZXt388e0+7vq3Q7veATyzGd4qDBfia65/X830BwLfBpaNu/6R/E7GXcAcPqzfcftN4gGum7L8u4DXjbvu7dUuYA1wRmu5vwVeNe76t9fn1VruQuBh465/O3xeFwE/bV7X0jsMfem46x/B5/WVhfx5bctrQZ4znMHVwOOa4cOAH/VnJNkBeAbw72Ooa1vN1K4vAPsnWZ5kp2aZ/x5DffM1bbuS7J1kx2b4fsD96Z1nWyymbVdVPaaqVlbVSuDdwD9W1aljqXB+Zvq8Jpp/fyTZF/gzeoG/5CzILrySfAw4FNgryS+AvweOB97TfDA3Ai9srfJY4BdVtaD/qObSrqr6bZJ3At+hd5L6c1X12bEUPos5fl6PBd6c5GbgFuDFVTX1ZP6CMI9/h4vCHNv1aOC1Sf5A70LJX1fV/3Vf9ej5bLIksUBvupakrhmGkoRhKEmAYShJgGEoSYBhKEmAYagZJLl7ko8m+UmS7zZdiR01yzor291CzXF7xyXZpzX+gSQPGnLdQ5OcN5/tDvn++yT5ZDO8KskR83iPN7a79dLCYxjqjyQJcC7wtaq6X1Wtpvd44L1HuNnjgNvCsKpeUFUL4ombqrq6qo5pRlfR6+JKS4xhqOkcBtxUVbd1vFpVV1XVe+G2PcCLkqxrXodMfYNByyR5ddNZ7WVJ/jnJMcDDgDObTkR3TnJhkoc1yz+5eY/Lknx52EYkObbZzhVJ3taavinJPzTv960kd2+m/2kz/p0kb06yqdWWK5LcGXgz8KymzmdN3eNrllvZDL8+yf8k+RK9x9hobef8Zo/7oiT7DdsmjdC4H472tfBewN8A7xowfzlNzyX0ni2ebIZX0jz8P2CZw4FvAsub8RXNzwtpdQDQHwf2ptezzUR7+Sn1HAqcN2XaPvS6pNqb3mOnFwBHNvMKeFozfDLwhmb4PODYZvjFwKZp2nUccGprO29k6+66rmiWXw1c3vwedgd+3F8O+DJw/2b4EcAF4/7MfdXCfDZZC0uS0+g9o3pTVR0E3Ak4Nb2vV7gFeMA0q820zBOAD1XVFoCa/bnkg+kdrq8fcvm+g+j11vybpg1n0nsu+lzgJnrBB/Bd4InN8COBI5vhj9Lrlmu+HgP8R7+dST7d/NwVOAQ4q3c2Auh1a6YxMww1nSuBo/sjVXVCkr2AyWbSy4FfAwfQO9Vy4zTvMdMyYevekmcz1+Xb683kD9XsltEL6m35O7iZrU83LWsNT1f3DsC1VbVqG7apEfCcoaZzAbAsyUta05a3hvcAflVVtwLPBXac5j1mWuaLwPOSLAdIsqKZfj2w2zTvczHwuCQTU5afzSXNens1XYYdS69n7UG+xe3/CayZYZmpdf6U3veJkORAYKKZ/jXgqOb8527A0wCq6jpgfZJnNOskyQFDtkkjZBjqjzR7TUfSC5P1Sb4NrAVe0yzyPuCvknyL3uHv5mneZtplqup84NPAZJJL6fXoDXAGcHr/Akqrlt/Q607qnCSX0fsSouk8Pskv+i965+1OotcZ6WXAuqr61CxNPxF4RdPee9Lr8HSqrwAP6l9AofcFUCuatrwE+GFT97qm1kubZS5qvcdzgOc37bkSePosdakDduElNZq91RuqqpKsoXcxxaC6g/CcoXS71fQu+oRe1/3PG2856pJ7hpKE5wwlCTAMJQkwDCUJMAwlCTAMJQmA/wc6HMjVACNE5gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# get PSF 95% containment radius\n", "energy_true = dataset.exposure.geom.axes[0].center\n", "psf_r95 = dataset.psf.containment_radius(\n", " fraction=0.95, energy_true=energy_true\n", ")\n", "# create mask_fit with margin based on PSF\n", "mask_fit = dataset.counts.geom.boundary_mask(psf_r95.max())\n", "dataset.mask_fit = mask_fit\n", "dataset.mask_fit.sum_over_axes().plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reading and writing masks\n", "\n", "`gammapy.maps` can directly read/write maps with boolean content as follows:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.254960Z", "iopub.status.busy": "2021-11-22T21:09:30.254601Z", "iopub.status.idle": "2021-11-22T21:09:30.277836Z", "shell.execute_reply": "2021-11-22T21:09:30.278021Z" } }, "outputs": [], "source": [ "# To save masks to disk\n", "mask_map.write(\"exclusion_mask.fits\", overwrite=\"True\")" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:09:30.279754Z", "iopub.status.busy": "2021-11-22T21:09:30.279436Z", "iopub.status.idle": "2021-11-22T21:09:30.282577Z", "shell.execute_reply": "2021-11-22T21:09:30.282795Z" } }, "outputs": [], "source": [ "# To read maps from disk\n", "mask_map = Map.read(\"exclusion_mask.fits\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.0" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 4 }