{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# :suppress:\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": [
    "# Edge behavior and interiors\n",
    "\n",
    "This notebook illustrates the edge behavior (when a grid point falls on the edge of a polygon) and how polygon interiors are treated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "Import regionmask and check the version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import regionmask\n",
    "\n",
    "regionmask.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib import colors as mplc\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some colors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color1 = \"#9ecae1\"\n",
    "color2 = \"#fc9272\"\n",
    "\n",
    "cmap1 = mplc.ListedColormap([color1])\n",
    "cmap2 = mplc.ListedColormap([color2])\n",
    "\n",
    "cmap12 = mplc.ListedColormap([color1, color2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methods\n",
    "\n",
    "Regionmask offers two backends (internally called \"methods\"*) to rasterize regions\n",
    "\n",
    "1. `rasterize`: fastest but only for equally-spaced grid, uses `rasterio.features.rasterize` internally.\n",
    "2. `shapely`: for irregular grids, uses `shapely.STRtree` internally.\n",
    "\n",
    "All methods use the `lon` and `lat` coordinates to determine if a grid cell is in a region. `lon` and `lat` are assumed to indicate the *center* of the grid cell. All methods have the same edge behavior and consider 'holes' in the regions. `regionmask` automatically determines which `method` to use.\n",
    "\n",
    "(2) and (3) subtract a tiny offset from `lon` and `lat` to achieve a edge behaviour consistent with (1). Due to [mapbox/rasterio#1844](https://github.com/mapbox/rasterio/issues/1844) this is unfortunately also necessary for (1).\n",
    "\n",
    "\n",
    "\\*Note that all \"methods\" yield the same results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Edge behavior\n",
    "\n",
    "The edge behavior determines how points that fall on the outline of a region are treated. It's easiest to see the edge behaviour in an example.\n",
    "\n",
    "### Example\n",
    "\n",
    "Define a region and a lon/ lat grid, such that some gridpoints lie exactly on the border:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outline = np.array([[-80.0, 44.0], [-80.0, 28.0], [-100.0, 28.0], [-100.0, 44.0]])\n",
    "\n",
    "region = regionmask.Regions([outline])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_US = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(\n",
    "    *(-161, -29, 2), *(75, 13, -2)\n",
    ")\n",
    "\n",
    "print(ds_US)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Let's create a mask with each of these methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_rasterize = region.mask(ds_US, method=\"rasterize\")\n",
    "mask_shapely = region.mask(ds_US, method=\"shapely\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "regionmask automatically detects which method to use, so there is no need to specify the ``method`` keyword.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the masked regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(1, 2, subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "opt = dict(add_colorbar=False, ec=\"0.5\", lw=0.5, transform=ccrs.PlateCarree())\n",
    "\n",
    "mask_rasterize.plot(ax=axes[0], cmap=cmap1, **opt)\n",
    "mask_shapely.plot(ax=axes[1], cmap=cmap2, **opt)\n",
    "\n",
    "\n",
    "for ax in axes:\n",
    "    ax = region.plot_regions(ax=ax, add_label=False)\n",
    "    ax.set_extent([-105, -75, 25, 49], ccrs.PlateCarree())\n",
    "    ax.coastlines(lw=0.5)\n",
    "\n",
    "    ax.plot(\n",
    "        ds_US.LON, ds_US.lat, \"*\", color=\"0.5\", ms=0.5, transform=ccrs.PlateCarree()\n",
    "    )\n",
    "\n",
    "axes[0].set_title(\"backend: rasterize\")\n",
    "axes[1].set_title(\"backend: shapely\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points indicate the grid cell centers (`lon` and `lat`), lines the grid cell borders, colored grid cells are selected to be part of the region. The top and right grid cells now belong to the region while the left and bottom grid cells do not. This choice is arbitrary but follows what `rasterio.features.rasterize` does. This avoids spurious columns of unassigned grid points as the following example shows.\n",
    "\n",
    "### SREX regions\n",
    "\n",
    "Create a global dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_GLOB = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(\n",
    "    *(-180, 181, 2), *(90, -91, -2)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srex = regionmask.defined_regions.srex\n",
    "\n",
    "srex_new = srex.mask(ds_GLOB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, 1, subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "opt = dict(add_colorbar=False, cmap=\"viridis_r\")\n",
    "\n",
    "srex_new.plot(ax=ax, ec=\"0.7\", lw=0.25, **opt)\n",
    "\n",
    "srex.plot_regions(ax=ax, add_label=False, line_kws=dict(lw=1))\n",
    "ax.set_extent([-135, -50, 24, 51], ccrs.PlateCarree())\n",
    "ax.coastlines(resolution=\"50m\", lw=0.25)\n",
    "\n",
    "ax.plot(\n",
    "    ds_GLOB.LON, ds_GLOB.lat, \"*\", color=\"0.5\", ms=0.5, transform=ccrs.PlateCarree()\n",
    ")\n",
    "\n",
    "\n",
    "sel = ((ds_GLOB.LON == -105) | (ds_GLOB.LON == -85)) & (ds_GLOB.LAT > 28)\n",
    "\n",
    "\n",
    "ax.plot(\n",
    "    ds_GLOB.LON.values[sel],\n",
    "    ds_GLOB.LAT.values[sel],\n",
    "    \"*\",\n",
    "    color=\"r\",\n",
    "    ms=0.5,\n",
    "    transform=ccrs.PlateCarree(),\n",
    ")\n",
    "\n",
    "\n",
    "ax.set_title(\"edge points are assigned to the left polygon\", fontsize=9);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not assigning the grid cells falling exactly on the border of a region (red points) would leave vertical stripes of unassigned cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Points at -180°E (0°E) and -90°N\n",
    "\n",
    "The described edge behaviour leads to a consistent treatment of points on the border. However, gridpoints at -180°E (or 0°E) and -90°N would *never* fall in any region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "From version 0.8 this applies only if ``wrap_lon`` is *not* set to ``False``. If wrap_lon is set to False `regionmask` assumes the coordinates are not lat and lon coordinates.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We exemplify this with a region spanning the whole globe and a coarse longitude/ latidude grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# almost 360 to avoid wrap-around for the plot\n",
    "lon_max = 360.0 - 1e-10\n",
    "outline_global = np.array([[0, 90], [0, -90], [lon_max, -90], [lon_max, 90]])\n",
    "\n",
    "region_global = regionmask.Regions([outline_global])\n",
    "\n",
    "lon = np.arange(0, 360, 30)\n",
    "lat = np.arange(90, -91, -30)\n",
    "\n",
    "LON, LAT = np.meshgrid(lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the masks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting `wrap_lon=False` turns this feature off\n",
    "mask_global_nontreat = region_global.mask(LON, LAT, wrap_lon=False)\n",
    "\n",
    "mask_global = region_global.mask(LON, LAT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And illustrate the issue:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.PlateCarree(central_longitude=180)\n",
    "f, axes = plt.subplots(1, 2, subplot_kw=dict(projection=proj))\n",
    "f.subplots_adjust(wspace=0.05)\n",
    "\n",
    "opt = dict(add_colorbar=False, ec=\"0.2\", lw=0.25, transform=ccrs.PlateCarree())\n",
    "\n",
    "ax = axes[0]\n",
    "mask_global_nontreat.plot(ax=ax, cmap=cmap1, x=\"lon\", y=\"lat\", **opt)\n",
    "\n",
    "ax.set_title(\"Not treating points at 0°E and -90°N\", size=6)\n",
    "ax.set_title(\"(a)\", loc=\"left\", size=6)\n",
    "\n",
    "ax = axes[1]\n",
    "mask_global.plot(ax=ax, cmap=cmap1, x=\"lon\", y=\"lat\", **opt)\n",
    "\n",
    "ax.set_title(\"Treating points at 0°E and -90°N\", size=6)\n",
    "ax.set_title(\"(b)\", loc=\"left\", size=6)\n",
    "\n",
    "\n",
    "for ax in axes:\n",
    "\n",
    "    ax = region_global.plot(\n",
    "        ax=ax,\n",
    "        line_kws=dict(lw=2, color=\"#b15928\"),\n",
    "        add_label=False,\n",
    "    )\n",
    "    ax.plot(LON, LAT, \"o\", color=\"0.3\", ms=1, transform=ccrs.PlateCarree(), zorder=5)\n",
    "    ax.spines[\"geo\"].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example the region spans the whole globe and there are gridpoints at 0°E and -90°N. Just applying the approach above leads to gridpoints that are not assigned to any region even though the region is global (as shown in a). Therefore, points at -180°E (or 0°E) and -90°N are treated specially (b):\n",
    "\n",
    "Points at -180°E (0°E) are mapped to 180°E (360°E). Points at -90°N are slightly shifted northwards (by 1 * 10 ** -10). Then it is tested if the shifted points belong to any region\n",
    "\n",
    "This means that (i) a point at -180°E is part of the region that is present at 180°E and not the one at -180°E (this is consistent with assigning points to the polygon *left* from it) and (ii) only the points at -90°N get assigned to the region above.\n",
    "\n",
    "This is illustrated in the figure below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outline_global1 = np.array([[-180.0, 60.0], [-180.0, -60.0], [0.0, -60.0], [0.0, 60.0]])\n",
    "outline_global2 = np.array([[0.0, 60.0], [0.0, -60.0], [180.0, -60.0], [180.0, 60.0]])\n",
    "\n",
    "region_global_2 = regionmask.Regions([outline_global1, outline_global2])\n",
    "\n",
    "mask_global_2regions = region_global_2.mask(lon, lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = region_global_2.plot(\n",
    "    line_kws=dict(color=\"#b15928\", zorder=3, lw=1.5),\n",
    "    add_label=False,\n",
    ")\n",
    "\n",
    "ax.plot(\n",
    "    LON, LAT, \"o\", color=\"0.3\", lw=0.25, ms=2, transform=ccrs.PlateCarree(), zorder=5\n",
    ")\n",
    "mask_global_2regions.plot(ax=ax, cmap=cmap12, **opt)\n",
    "\n",
    "ax.set_title(\"Points at -180°E are mapped to 180°E\", size=6)\n",
    "\n",
    "ax.spines[\"geo\"].set_lw(0.25)\n",
    "ax.spines[\"geo\"].set_zorder(1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    ":::{note}\n",
    "\n",
    "This only applies if the border of the region falls *exactly* on the point. One way to avoid the problem is to calculate the fractional overlap of each gridpoint with the regions, see [mask_3D_frac_approx](mask_3D_frac_approx.ipynb).\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polygon interiors\n",
    "\n",
    "`Polygons` can have interior boundaries ('holes'). regionmask unmasks these regions.\n",
    "\n",
    "### Example\n",
    "\n",
    "Let's test this on an example and define a `region_with_hole`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interior = np.array(\n",
    "    [\n",
    "        [-86.0, 40.0],\n",
    "        [-86.0, 32.0],\n",
    "        [-94.0, 32.0],\n",
    "        [-94.0, 40.0],\n",
    "    ]\n",
    ")\n",
    "\n",
    "poly = Polygon(outline, holes=[interior])\n",
    "\n",
    "region_with_hole = regionmask.Regions([poly])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_hole_rasterize = region_with_hole.mask(ds_US, method=\"rasterize\")\n",
    "mask_hole_shapely = region_with_hole.mask(ds_US, method=\"shapely\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(1, 2, subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "opt = dict(add_colorbar=False, ec=\"0.5\", lw=0.5)\n",
    "\n",
    "mask_hole_rasterize.plot(ax=axes[0], cmap=cmap1, **opt)\n",
    "mask_hole_shapely.plot(ax=axes[1], cmap=cmap2, **opt)\n",
    "\n",
    "for ax in axes:\n",
    "    region_with_hole.plot_regions(ax=ax, add_label=False, line_kws=dict(lw=1))\n",
    "\n",
    "    ax.set_extent([-105, -75, 25, 49], ccrs.PlateCarree())\n",
    "    ax.coastlines(lw=0.25)\n",
    "\n",
    "    ax.plot(\n",
    "        ds_US.LON, ds_US.lat, \"o\", color=\"0.5\", ms=0.5, transform=ccrs.PlateCarree()\n",
    "    )\n",
    "\n",
    "axes[0].set_title(\"backend: rasterize\")\n",
    "axes[1].set_title(\"backend: shapely\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the edge behavior of the interior is inverse to the edge behavior of the outerior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Caspian Sea\n",
    "\n",
    "The Caspian Sea is defined as polygon interior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "land110 = regionmask.defined_regions.natural_earth_v5_0_0.land_110\n",
    "\n",
    "mask_land110 = land110.mask(ds_GLOB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, 1, subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "mask_land110.plot(ax=ax, cmap=cmap2, add_colorbar=False)\n",
    "\n",
    "\n",
    "ax.set_extent([15, 75, 25, 50], ccrs.PlateCarree())\n",
    "ax.coastlines(resolution=\"50m\", lw=0.5)\n",
    "\n",
    "ax.plot(\n",
    "    ds_GLOB.LON, ds_GLOB.lat, \".\", color=\"0.5\", ms=0.5, transform=ccrs.PlateCarree()\n",
    ")\n",
    "\n",
    "\n",
    "ax.text(52, 43.5, \"Caspian Sea\", transform=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_title(\"Polygon interiors are unmasked\");"
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}