{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Plotting with CartoPy and GeoPandas\n",
    "\n",
    "Converting between GeoPandas and CartoPy for visualizing data.\n",
    "\n",
    "[CartoPy](http://scitools.org.uk/cartopy/) is a Python library\n",
    "that specializes in creating geospatial\n",
    "visualizations. It has a slightly different way of representing\n",
    "Coordinate Reference Systems (CRS) as well as constructing plots.\n",
    "This example steps through a round-trip transfer of data\n",
    "between GeoPandas and CartoPy.\n",
    "\n",
    "First we'll load in the data using GeoPandas.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import geopandas\n",
    "from cartopy import crs as ccrs\n",
    "from geodatasets import get_path\n",
    "\n",
    "path = get_path(\"naturalearth.land\")\n",
    "df = geopandas.read_file(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll visualize the map using GeoPandas\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting with CartoPy\n",
    "=====================\n",
    "\n",
    "Cartopy also handles Shapely objects well, but it uses a different system for\n",
    "CRS. To plot this data with CartoPy, we'll first need to project it into a\n",
    "new CRS. We'll use a CRS defined within CartoPy and use the GeoPandas\n",
    "``to_crs`` method to make the transformation.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the CartoPy CRS object.\n",
    "crs = ccrs.AzimuthalEquidistant()\n",
    "\n",
    "# This can be converted into a `proj4` string/dict compatible with GeoPandas\n",
    "crs_proj4 = crs.proj4_init\n",
    "df_ae = df.to_crs(crs_proj4)\n",
    "\n",
    "# Here's what the plot looks like in GeoPandas\n",
    "df_ae.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our data is in a CRS based off of CartoPy, we can easily\n",
    "plot it.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(subplot_kw={\"projection\": crs})\n",
    "ax.add_geometries(df_ae[\"geometry\"], crs=crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we could have easily done this with an EPSG code like so:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs_epsg = ccrs.epsg(\"3857\")\n",
    "df_epsg = df.to_crs(epsg=\"3857\")\n",
    "\n",
    "# Generate a figure with two axes, one for CartoPy, one for GeoPandas\n",
    "fig, axs = plt.subplots(1, 2, subplot_kw={\"projection\": crs_epsg}, figsize=(10, 5))\n",
    "# Make the CartoPy plot\n",
    "axs[0].add_geometries(\n",
    "    df_epsg[\"geometry\"], crs=crs_epsg, facecolor=\"white\", edgecolor=\"black\"\n",
    ")\n",
    "# Make the GeoPandas plot\n",
    "df_epsg.plot(ax=axs[1], color=\"white\", edgecolor=\"black\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CartoPy to GeoPandas\n",
    "====================\n",
    "\n",
    "Next we'll perform a CRS projection in CartoPy, and then convert it\n",
    "back into a GeoPandas object.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs_new = ccrs.AlbersEqualArea()\n",
    "new_geometries = [\n",
    "    crs_new.project_geometry(ii, src_crs=crs) for ii in df_ae[\"geometry\"].values\n",
    "]\n",
    "\n",
    "fig, ax = plt.subplots(subplot_kw={\"projection\": crs_new})\n",
    "ax.add_geometries(new_geometries, crs=crs_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've created new Shapely objects with the CartoPy CRS,\n",
    "we can use this to create a GeoDataFrame.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_aea = geopandas.GeoDataFrame(\n",
    "    df.drop(columns=\"geometry\"), geometry=new_geometries, crs=crs_new.proj4_init\n",
    ")\n",
    "df_aea.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even combine these into the same figure. Here we'll plot the\n",
    "shapes of the countries with CartoPy. We'll then calculate the centroid\n",
    "of each with GeoPandas and plot it on top.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "# Generate a CartoPy figure and add the countries to it\n",
    "fig, ax = plt.subplots(subplot_kw={\"projection\": crs_new})\n",
    "ax.add_geometries(new_geometries, crs=crs_new)\n",
    "\n",
    "# Calculate centroids and plot\n",
    "df_aea_centroids = df_aea.geometry.centroid\n",
    "# Need to provide \"zorder\" to ensure the points are plotted above the polygons\n",
    "df_aea_centroids.plot(ax=ax, markersize=5, color=\"r\", zorder=10)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}