{
 "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\")\n",
    "\n",
    "# turn off pandas html repr:\n",
    "# does not gracefully survive the ipynb -> rst -> html conversion\n",
    "import pandas as pd\n",
    "\n",
    "pd.set_option(\"display.notebook_repr_html\", False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "{{ prolog }}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with geopandas (shapefiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "regionmask includes support for regions defined as geopandas GeoDataFrame. These are often shapefiles, which can be opened in the formats `.zip`, `.shp`, `.geojson` etc. with `geopandas.read_file(url_or_path).`\n",
    "\n",
    "There are two possibilities:\n",
    "\n",
    "1. Directly create a mask from a geopandas GeoDataFrame or GeoSeries using `mask_geopandas` or `mask_3D_geopandas`. \n",
    "2. Convert a GeoDataFrame to a `Regions` object (regionmask's internal data container) using `from_geopandas`.\n",
    "\n",
    "As always, start with the imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import geopandas as gp\n",
    "import matplotlib.patheffects as pe\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pooch\n",
    "\n",
    "import regionmask\n",
    "\n",
    "regionmask.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Opening an example shapefile\n",
    "\n",
    "The U.S. Geological Survey (USGS) offers a shapefile containing the outlines of continens [1]. We use the library pooch to locally cache the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = pooch.retrieve(\n",
    "    \"https://pubs.usgs.gov/of/2006/1187/basemaps/continents/continents.zip\", None\n",
    ")\n",
    "\n",
    "continents = gp.read_file(\"zip://\" + file)\n",
    "\n",
    "display(continents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a mask from a GeoDataFrame\n",
    "\n",
    "`mask_geopandas` and `mask_3D_geopandas` allow to directly create a mask from a GeoDataFrame or GeoSeries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.arange(-180, 180)\n",
    "lat = np.arange(-90, 90)\n",
    "\n",
    "mask = regionmask.mask_geopandas(continents, lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the new mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "mask.plot(\n",
    "    ax=ax,\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    add_colorbar=False,\n",
    ")\n",
    "\n",
    "ax.coastlines(color=\"0.1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly a 3D boolean mask can be created from a GeoDataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_3D = regionmask.mask_3D_geopandas(continents, lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and plotted:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import colors as mplc\n",
    "\n",
    "cmap1 = mplc.ListedColormap([\"none\", \"#9ecae1\"])\n",
    "\n",
    "f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "\n",
    "mask_3D.sel(region=0).plot(\n",
    "    ax=ax,\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    add_colorbar=False,\n",
    "    cmap=cmap1,\n",
    ")\n",
    "\n",
    "ax.coastlines(color=\"0.1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "\n",
    "Set ``regionmask.mask_3D_geopandas(..., overlap=True)`` if some of the regions overlap. See the tutorial on [overlapping regions](overlap) for details.\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Convert GeoDataFrame to a Regions object\n",
    "\n",
    "Creating a `Regions` object with `regionmask.from_geopandas` requires a GeoDataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "continents_regions = regionmask.from_geopandas(continents)\n",
    "continents_regions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This creates default names (`\"Region0\"`, ..., `\"RegionN\"`) and abbreviations (`\"r0\"`, ..., `\"rN\"`).\n",
    "\n",
    "However, it is often advantageous to use columns of the GeoDataFrame as names and abbrevs. If no column with abbreviations is available, you can use `abbrevs='_from_name'`, which creates unique abbreviations using the names column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "continents_regions = regionmask.from_geopandas(\n",
    "    continents, names=\"CONTINENT\", abbrevs=\"_from_name\", name=\"continent\"\n",
    ")\n",
    "continents_regions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual the newly created `Regions` object can be plotted on a world map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "text_kws = dict(\n",
    "    bbox=dict(color=\"none\"),\n",
    "    path_effects=[pe.withStroke(linewidth=2, foreground=\"w\")],\n",
    "    color=\"#67000d\",\n",
    "    fontsize=9,\n",
    ")\n",
    "\n",
    "continents_regions.plot(label=\"name\", add_coastlines=False, text_kws=text_kws);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to create mask a mask for arbitrary latitude/ longitude grids:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.arange(0, 360)\n",
    "lat = np.arange(-90, 90)\n",
    "\n",
    "mask = continents_regions.mask(lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can then be plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))\n",
    "h = mask.plot(\n",
    "    ax=ax,\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    cmap=\"Reds\",\n",
    "    add_colorbar=False,\n",
    "    levels=np.arange(-0.5, 8),\n",
    ")\n",
    "\n",
    "cbar = plt.colorbar(h, shrink=0.625, pad=0.025, aspect=12)\n",
    "cbar.set_ticks(np.arange(8))\n",
    "cbar.set_ticklabels(continents_regions.names)\n",
    "\n",
    "ax.coastlines(color=\"0.2\")\n",
    "\n",
    "continents_regions.plot_regions(add_label=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] Environmental Systems Research , Inc. (ESRI), 20020401, World Continents: ESRI Data & Maps 2002, Environmental Systems Research Institute, Inc. (ESRI), Redlands, California, USA."
   ]
  }
 ],
 "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}