{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from matplotlib import rcParams\n",
    "\n",
    "rcParams[\"figure.dpi\"] = 300\n",
    "rcParams[\"font.size\"] = 8\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "{{ prolog }}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overlapping regions\n",
    "\n",
    "Two or more on regions can share the same area - they overlap. One example are region 3 and 4 of the [PRUDENCE regions](prudence-regions). This notebook illustrates how overlapping regions can be treated in regionmask. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In short"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    ":::{important}\n",
    "\n",
    "Starting with v0.11.0 regionmask checks if gridpoints belong to more than one region (unless ``overlap=False``). Thus, it creates the correct mask (or raises an error) per default.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In detail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When creating your own `Regions` you can tell regionmask if they are overlapping:\n",
    "\n",
    "```python\n",
    "region = regionmask.Regions(..., overlap=overlap)\n",
    "region = regionmask.from_geopandas(..., overlap=overlap)\n",
    "```\n",
    "\n",
    "where ``overlap`` must be one of ``None`` (default), ``True`` or ``False``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For `overlap=None` regionmask checks if any grid point belongs to more than one region. If so, 2D masks raise and error (because overlapping regions cannot be represented by a 2D mask). `mask_3D(...)` returns the correct mask.\n",
    "- For `overlap=True` 2D masks raise an error and ``mask_3D(...)`` returns the correct mask. \n",
    "- If you have two overlapping regions and `overlap=False` regionmask will _silently_ assign the gridpoints of the overlapping regions to the one with the higher number.\n",
    "- Setting ``overlap`` to ``True`` or ``False`` is (slightly) faster than the default of ``None``.\n",
    "- `overlap` is correctly set for regions in `regionmask.defined_regions`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "To illustrate the problem we construct two regions in North America that partially overlap. One is horizontal, the other vertical.\n",
    "\n",
    "**Preparation**\n",
    "\n",
    "Imports & regionmask version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.patheffects as pe\n",
    "import numpy as np\n",
    "from matplotlib import colors as mplc\n",
    "from shapely.geometry import box\n",
    "\n",
    "import regionmask\n",
    "\n",
    "regionmask.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_US = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(\n",
    "    *(-160, -29, 2), *(76, 13, -2)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some colors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = mplc.ListedColormap([\"none\", \"#9ecae1\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define helper function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_region_vh(mask, region):\n",
    "\n",
    "    fg = mask.plot(\n",
    "        subplot_kws=dict(projection=ccrs.PlateCarree()),\n",
    "        col=\"region\",\n",
    "        cmap=cmap,\n",
    "        add_colorbar=False,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        ec=\"0.5\",\n",
    "        lw=0.25,\n",
    "    )\n",
    "\n",
    "    for ax in fg.axs.flatten():\n",
    "        region[[0]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#6a3d9a\"))\n",
    "        region[[1]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#ff7f00\"))\n",
    "\n",
    "        ax.set_extent([-105, -75, 25, 55], ccrs.PlateCarree())\n",
    "        ax.plot(\n",
    "            ds_US.LON, ds_US.lat, \"*\", color=\"0.5\", ms=0.5, transform=ccrs.PlateCarree()\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define overlapping polygons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords_v = box(-100, 28, -90, 50)\n",
    "coords_h = box(-100, 40, -80, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Default behavior (`overlap=None`)\n",
    "\n",
    "First test what happens if we keep the default value of `overlap=None`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_vh = regionmask.Regions([coords_v, coords_h])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_vh = region_vh.mask_3D(ds_US)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the masked regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_region_vh(mask_vh, region_vh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The small gray points show the gridpoint center and the vertical and horizontal lines are the gridpoint boundaries. The colored rectangles are the two regions. Here regionmask detects that certain gridpoints belong to two regions and correctly assigns them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting `overlap=False`\n",
    "\n",
    "We test what happens for `overlap=False`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_vh_overlap_false = regionmask.Regions([coords_v, coords_h], overlap=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_vh_overlap_false = region_vh_overlap_false.mask_3D(ds_US)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the masked regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_region_vh(mask_vh_overlap_false, region_vh_overlap_false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that only the gridpoints in the lower part of the vertical (magenta) region were assigned to it. All gridpoints of the overlapping part are now assigned to the horizontal (orange) region - the gridpoints are assigned to the region with the higher number. By switching the order of the regions you could have the common gridpoints assigned to the vertical region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting `overlap=True`\n",
    "\n",
    "Passing `overlap=True` has the same effect as the default:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_overlap = regionmask.Regions([coords_v, coords_h], overlap=True)\n",
    "\n",
    "region_overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it says `overlap:  True` - and we can again create a mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_overlap = region_overlap.mask_3D(ds_US)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_region_vh(mask_overlap, region_vh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the gridpoints in the overlapping part are assigned to both regions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PRUDENCE regions\n",
    "\n",
    "The PRUDENCE regions are a real-world example of overlapping areas. The prudence regions already set `overlap=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prudence = regionmask.defined_regions.prudence\n",
    "prudence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regions 3 and 4 overlap in Western France:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertConformal(central_longitude=10)\n",
    "\n",
    "text_kws = dict(\n",
    "    bbox=dict(color=\"none\"),\n",
    "    path_effects=[pe.withStroke(linewidth=3, foreground=\"w\")],\n",
    "    color=\"#67000d\",\n",
    ")\n",
    "\n",
    "ax = prudence.plot(\n",
    "    projection=proj, text_kws=text_kws, resolution=\"50m\", line_kws=dict(lw=0.75)\n",
    ")\n",
    "\n",
    "\n",
    "ax.set_extent([-10.0, 30.0, 40.0, 55.0], ccrs.PlateCarree())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create mask of PRUDENCE regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.arange(-12, 33, 0.25)\n",
    "lat = np.arange(72, 33, -0.25)\n",
    "\n",
    "mask_prudence = prudence.mask_3D(lon, lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.LambertConformal(central_longitude=10)\n",
    "\n",
    "fg = mask_prudence.sel(region=[3, 4]).plot(\n",
    "    subplot_kws=dict(projection=proj),\n",
    "    col=\"region\",\n",
    "    cmap=cmap,\n",
    "    add_colorbar=False,\n",
    "    transform=ccrs.PlateCarree(),\n",
    ")\n",
    "\n",
    "\n",
    "for ax in fg.axs.flatten():\n",
    "    regionmask.defined_regions.prudence.plot(\n",
    "        ax=ax, add_label=False, resolution=\"50m\", line_kws=dict(lw=0.75)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As above the gridpoints below the overlapping part is now assigned to both regions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}