{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing spatial data with Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import geopandas\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = geopandas.read_file(\"data/ne_110m_admin_0_countries.zip\")\n",
    "cities = geopandas.read_file(\"data/ne_110m_populated_places.zip\")\n",
    "rivers = geopandas.read_file(\"data/ne_50m_rivers_lake_centerlines.zip\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GeoPandas visualization functionality\n",
    "\n",
    "GeoPandas itself provides some visualization functionality, and together with matplotlib for further customization, you can already get decent results for visualizing vector data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Basic plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adjusting the figure size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries.plot(figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Removing the box / x and y coordinate labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = countries.plot(figsize=(15, 15))\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Coloring based on column values\n",
    "\n",
    "Let's first create a new column with the GDP per capita:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = countries[(countries['pop_est'] >0 ) & (countries['name'] != \"Antarctica\")].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries['gdp_per_cap'] = countries['gdp_md_est'] / countries['pop_est'] * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and now we can use this column to color the polygons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = countries.plot(figsize=(15, 15), column='gdp_per_cap')\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a classification scheme to bin the values (using [`mapclassify`](https://pysal.org/mapclassify/)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = countries.plot(figsize=(15, 15), column='gdp_per_cap', scheme='quantiles', legend=True)\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Combining different dataframes on a single plot\n",
    "\n",
    "The `.plot` method returns a matplotlib Axes object, which can then be re-used to add additional layers to that plot with the `ax=` keyword:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = countries.plot(figsize=(15, 15))\n",
    "cities.plot(ax=ax, color='red', markersize=10)\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = countries.plot(edgecolor='k', facecolor='none', figsize=(15, 10))\n",
    "rivers.plot(ax=ax)\n",
    "cities.plot(ax=ax, color='C1')\n",
    "ax.set(xlim=(-20, 60), ylim=(-40, 40))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a background map with contextily"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The contextily package allow to easily add a web-tile based backgroubd (basemap) to your GeoPandas plots.\n",
    "\n",
    "Currently, the only requirement is that your data is already in the WebMercator projection (EPSG:3857)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# selecting the cities in Europe\n",
    "cities_europe = cities[cities.within(countries[countries['continent'] == 'Europe'].unary_union)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# converting to WebMercator\n",
    "cities_europe2 = cities_europe.to_crs(epsg=3857)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = cities_europe2.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import contextily"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = cities_europe2.plot(figsize=(10, 6))\n",
    "contextily.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = cities_europe2.plot(figsize=(10, 6))\n",
    "contextily.add_basemap(ax, url=contextily.providers.Stamen.TonerLite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection-aware maps with Cartopy\n",
    "\n",
    "Cartopy is the base matplotlib cartographic library, and it is used by `geoplot` under the hood to provide projection-awareness (http://scitools.org.uk/cartopy/docs/latest/index.html).\n",
    "\n",
    "The following example is taken from the docs: http://geopandas.readthedocs.io/en/latest/gallery/cartopy_convert.html#sphx-glr-gallery-cartopy-convert-py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy import crs as ccrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the CartoPy CRS object.\n",
    "crs = ccrs.AlbersEqualArea()\n",
    "\n",
    "# This can be converted into a `proj4` string/dict compatible with GeoPandas\n",
    "crs_proj4 = crs.proj4_init\n",
    "countries_ae = countries.to_crs(crs_proj4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's what the plot looks like in GeoPandas\n",
    "countries_ae.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's what the plot looks like when plotting with cartopy\n",
    "fig, ax = plt.subplots(subplot_kw={'projection': crs})\n",
    "ax.add_geometries(countries_ae['geometry'], crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's what the plot looks like when plotting with cartopy and geopandas combined\n",
    "fig, ax = plt.subplots(subplot_kw={'projection': crs})\n",
    "countries_ae['geometry'].plot(ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using `geoplot`\n",
    "\n",
    "The `geoplot` packages provides some additional functionality compared to the basic `.plot()` method on GeoDataFrames:\n",
    "\n",
    "- High-level plotting API (with more plot types as geopandas)\n",
    "- Native projection support through cartopy\n",
    "\n",
    "https://residentmario.github.io/geoplot/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geoplot\n",
    "import geoplot.crs as gcrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={\n",
    "    'projection': gcrs.Orthographic(central_latitude=40.7128, central_longitude=-74.0059)\n",
    "})\n",
    "geoplot.choropleth(countries, hue='gdp_per_cap', projection=gcrs.Orthographic(), ax=ax,\n",
    "                   cmap='magma', linewidth=0.5, edgecolor='white')\n",
    "ax.set_global()\n",
    "ax.spines['geo'].set_visible(True)\n",
    "#ax.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive web-based visualizations\n",
    "\n",
    "There are nowadays many libraries that target interactive web-based visualizations and that can handle geospatial data. Some packages with an example for each:\n",
    "\n",
    "- Bokeh: https://bokeh.pydata.org/en/latest/docs/gallery/texas.html\n",
    "- GeoViews (other interface to Bokeh/matplotlib): http://geo.holoviews.org\n",
    "- Altair: https://altair-viz.github.io/gallery/choropleth.html\n",
    "- Plotly: https://plot.ly/python/#maps\n",
    "- ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another popular javascript library for online maps is [Leaflet.js](https://leafletjs.com/), and this has python bindings in the [folium](https://github.com/python-visualization/folium) and [ipyleaflet](https://github.com/jupyter-widgets/ipyleaflet) packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example with ipyleaflet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipyleaflet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = ipyleaflet.Map(center=[48.8566, 2.3429], zoom=6)\n",
    "\n",
    "layer = ipyleaflet.GeoJSON(data=cities.__geo_interface__)\n",
    "m.add_layer(layer)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = ipyleaflet.Map(center=[48.8566, 2.3429], zoom=3)\n",
    "geo_data = ipyleaflet.GeoData(\n",
    "    geo_dataframe = countries,\n",
    "    style={'color': 'black', 'fillColor': '#3366cc', 'opacity':0.05, 'weight':1.9, 'dashArray':'2', 'fillOpacity':0.6},\n",
    "    hover_style={'fillColor': 'red' , 'fillOpacity': 0.2},\n",
    "    name = 'Countries')\n",
    "m.add_layer(geo_data)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More: https://ipyleaflet.readthedocs.io/en/latest/api_reference/geodata.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example with folium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map([48.8566, 2.3429], zoom_start=6, tiles=\"OpenStreetMap\")\n",
    "folium.GeoJson(countries).add_to(m)\n",
    "folium.GeoJson(cities).add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map([0, 0], zoom_start=1)\n",
    "folium.Choropleth(geo_data=countries, data=countries, columns=['iso_a3', 'gdp_per_cap'],\n",
    "             key_on='feature.properties.iso_a3', fill_color='BuGn', highlight=True).add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\" style=\"font-size:120%\">\n",
    "\n",
    "**NOTE**: <br>\n",
    "\n",
    "Making a quick plot using Folium is now also available as the `.explore()` method on a GeoDataFrame or GeoSeries.\n",
    "\n",
    "See https://geopandas.org/en/stable/docs/user_guide/interactive_mapping.html for more examples.\n",
    "</div>"
   ]
  }
 ],
 "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.10.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}