geo-vis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geo-Spatial Visualisation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn the basics of geospatial visualisation using code. If you're new to geospatial analysis, you should look at the introduction page first.\n",
    "\n",
    "You should be aware when following this chapter that installing geographic analysis packages isn't always the easiest and things can and do go wrong! (Some geospatial analysis courses recommend running everything in a Docker container.)\n",
    "\n",
    "### Imports and packages\n",
    "\n",
    "We'll be using [**geopandas**](https://geopandas.org/index.html), the go-to package for vector spatial analysis in Python. The easiest way to install this package is using `conda install geopandas`; if you want to install it via pip then look at the [install instructions](https://geopandas.org/install.html). \n",
    "\n",
    "Let's import some of the packages we'll be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import os\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "from pathlib import Path\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "# For this page, use data limits and bigger default fig size\n",
    "plt.style.use(\n",
    "    {\n",
    "        \"axes.autolimit_mode\": \"data\",\n",
    "        \"figure.figsize\": (10, 8),\n",
    "        \"figure.dpi\": 125,\n",
    "        \"patch.linewidth\": 0.2,\n",
    "    }\n",
    ")\n",
    "# Turn off warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maps\n",
    "\n",
    "If you've looked at the introductory page on geospatial analysis, you'll know it's easy to make basic maps: you just need to load a shapefile and use the `.plot` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To run this example, you will need to download the files at this url:\n",
    "# https://github.com/aeturrell/coding-for-economists/tree/main/data/geo/uk_lad\n",
    "# and save them in the path 'data/geo/uk_lad' (path relative to where you're running the code)\n",
    "df = gpd.read_file(\n",
    "    Path(\"data/geo/uk_lad/Local_Authority_Districts__May_2020__UK_BUC.shp\")\n",
    ")\n",
    "df.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it goes, this is not very attractive, so let's see some options for customisation that will make it a little better. It's rare that you'll want to include the axes on maps, and these can be turned off by turning everything to do with the axes off. There are two ways to do further manipulations of the figure axis: calling plot returns an axis object or we can create one and then pass `ax=ax_name` to plot as a keyword argument. Colour can be changed using the `color=` keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = df.plot(color=\"#2ca25f\")\n",
    "ax.axis(\"off\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lines that divide up the different local authority districts are faint. They can be controlled with the `edgecolor` and `linewidth` keywords. We can also change the background using the `fig.patch.set_facecolor` method, and add a scale using an extra package, [**matplotlib-scalebar**](https://github.com/ppinard/matplotlib-scalebar)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_scalebar.scalebar import ScaleBar\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "df.plot(color=\"#2ca25f\", edgecolor=\"k\", linewidth=0.2, facecolor=\"blue\", ax=ax)\n",
    "ax.axis(\"off\")\n",
    "fig.patch.set_facecolor(\"#9ecae1\")\n",
    "# Create scale bar\n",
    "scalebar = ScaleBar(\n",
    "    1,\n",
    "    box_alpha=0,\n",
    "    location=\"lower right\",\n",
    "    length_fraction=0.25,\n",
    "    font_properties={\"size\": 12},\n",
    ")\n",
    "ax.add_artist(scalebar)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choropleths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A choropleth map shows different areas in colours according to a statistic that represents an aggregate summary of a geographic characteristic within each area. Population density or per-capita income are good examples of characteristics. The statistic shown might be unique values, equal intervals, quantiles, or the Fisher-Jenks natural breaks.\n",
    "\n",
    "First, though, let's create a basic choropleth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pay = pd.read_csv(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/data/geo/ashe_lad_median_pay_2020.csv\"\n",
    ")\n",
    "pay = pay.rename(columns={\"lad\": \"LAD20CD\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.merge(pay, on=[\"LAD20CD\"], how=\"inner\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col = \"Log median weekly pay (2020 GBP)\"\n",
    "df[col] = np.log(df[\"Median weekly pay (2020 GBP)\"])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(col, loc=\"left\")\n",
    "df.plot(\n",
    "    ax=ax,\n",
    "    column=col,\n",
    "    legend=True,\n",
    "    legend_kwds={\"label\": \"\", \"shrink\": 0.6},\n",
    "    vmin=round(df[col].min()),\n",
    "    vmax=round(df[col].max()),\n",
    ")\n",
    "ax.axis(\"off\")\n",
    "plt.tight_layout()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This used **geopandas**. There's a dedicated plotting tool called [**geoplot**](https://residentmario.github.io/geoplot/index.html) as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geoplot as gplt\n",
    "import geoplot.crs as gcrs\n",
    "\n",
    "gplt.choropleth(\n",
    "    df.to_crs(\"EPSG:4326\"),\n",
    "    hue=col,\n",
    "    projection=gcrs.AlbersEqualArea(),\n",
    "    cmap=\"viridis\",\n",
    "    legend=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to create choropleths is to split the variable into a distinct number of ranges according to a scheme. In the below, we use `scheme='Quantiles'` with `k=4` to produce a choropleth with four distinct groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 10))\n",
    "ax.set_title(col, loc=\"left\")\n",
    "ax.axis(\"off\")\n",
    "df.plot(ax=ax, column=col, legend=True, scheme=\"Quantiles\", k=4, legend_kwds={\"loc\": 2});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A third kind of choropleth has distinct levels based on pre-existing categories. Our data doesn't have any of those, so let's generate some just to show how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"cat_col\"] = df[\"LAD20CD\"].apply(lambda x: x[0])\n",
    "df.iloc[:5, -3:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "df.plot(\n",
    "    column=\"cat_col\",\n",
    "    categorical=True,\n",
    "    ax=ax,\n",
    "    legend=True,\n",
    "    legend_kwds={\"loc\": 1, \"frameon\": True},\n",
    ")\n",
    "ax.set_axis_off()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is useful for, for example, plotting streets of different types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cartogram\n",
    "\n",
    "A cartogram is a thematic map in which the geographic size of regions is altered to be proportional to a variable. The shape of the region is warped or shrunk. They can be especially useful when trying to refectly the fact that regions with the largest area do not always have proportionate variable values. Actually, due to the tendency of the shape of political regions to reflect choices made 100s or 1000s of years previously, and for urban areas to be under separate political arrangements, quite often economic variables are anti-correlated with areas.\n",
    "\n",
    "A cartogram of pay in Wales demonstrates this. Some areas with higer median incomes, such as Monmouthshire and Conwy, almost completely fill their potential region areas. But others, including Blaenau Gwent and Powys, are shown much smaller than their actual areas.\n",
    "\n",
    "The important part of the plot below is the `gplt.cartogram` but, along with other bits to make the plot look better, we're adding `gplt.polyplot` to show what the true size of each region is when it is not proportional to another variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.to_crs(\"EPSG:4326\")\n",
    "# Get a representative point for each region to annotate\n",
    "df[\"coords\"] = df[\"geometry\"].representative_point().apply(lambda x: x.coords[:][0])\n",
    "df_wales = df[df[\"LAD20CD\"].str.contains(\"W\")].fillna(0.0)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 5), dpi=125)\n",
    "gplt.cartogram(df_wales, scale=\"Median weekly pay (2020 GBP)\", ax=ax)\n",
    "gplt.polyplot(df_wales, facecolor=\"lightgray\", edgecolor=\"white\", linewidth=0.5, ax=ax)\n",
    "# Add text annotation to the largest polygons\n",
    "for idx, row in df_wales.iterrows():\n",
    "    if row[\"geometry\"].area > np.quantile(df.area, q=0.7):\n",
    "        ax.annotate(\n",
    "            text=row[\"LAD20NM\"],\n",
    "            xy=row[\"coords\"],\n",
    "            horizontalalignment=\"center\",\n",
    "            weight=\"bold\",\n",
    "            fontsize=8,\n",
    "            color=\"black\",\n",
    "        )\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadtree\n",
    "\n",
    "A quadtree is a tree data structure that splits a space into increasingly small rectangular fractals. This plot takes a sequence of point or polygonal geometries as input and builds a choropleth out of their centroids. Quadtrees are good at illustrating density, and are more flexible than a conventional choropleth: remember that choropleths can be the result of binning point occurrences into geographical regions or of data that are already aggregated to the region level. Quadtree is not a replacement for the latter, because the data are already aggregated. But, if you have point data, quadtree allows you to aggregate them *not* according to a pre-defined geography. Given pre-defined geographies such as Local Authority Districts may not be useful for the question you're thinking about (or worse could be downright misleading), this is a very helpful property.\n",
    "\n",
    "We'll use an example from the **geoplot** documentation to illustrate them. The most basic layer just turns a series of points into a quadtree. We'll use the lats and longs of collisions in New York:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "collisions = gpd.read_file(gplt.datasets.get_path(\"nyc_collision_factors\"))\n",
    "gplt.quadtree(collisions, nmax=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's project this onto a background of NY's boroughs. Because this is computationally expensive, we'll use the `simplify` method to reduce the complexity of the geometries we're using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boroughs = gpd.read_file(gplt.datasets.get_path(\"nyc_boroughs\"))\n",
    "gplt.quadtree(\n",
    "    collisions,\n",
    "    nmax=1,\n",
    "    projection=gcrs.AlbersEqualArea(),\n",
    "    clip=boroughs.simplify(0.001),\n",
    "    facecolor=\"lightgray\",\n",
    "    edgecolor=\"white\",\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can enjoy the best of a choropleth's ability to show us magnitudes *alongside* the ability of quadtree to show us geographic density through smaller rectangles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gplt.quadtree(\n",
    "    collisions,\n",
    "    nmax=1,\n",
    "    agg=np.mean,\n",
    "    projection=gcrs.AlbersEqualArea(),\n",
    "    clip=boroughs,\n",
    "    hue=\"NUMBER OF PEDESTRIANS INJURED\",\n",
    "    cmap=\"plasma\",\n",
    "    edgecolor=\"k\",\n",
    "    legend=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KDE plot\n",
    "\n",
    "You probably know kernel density estimation from 1D distribution functions, but there's no reason not to have the same fun in 2D. Taking the collisions data again, below is an example for New York. The `thresh=0` keyword argument just tells the KDE estimation to leave no empty whitespace where the estimated values are at their lowest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = gplt.polyplot(boroughs, projection=gcrs.AlbersEqualArea(), linewidth=0.5, zorder=1)\n",
    "gplt.kdeplot(\n",
    "    collisions, cmap=\"plasma\", shade=True, thresh=0, clip=boroughs, ax=ax, zorder=0\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatio-Temporal Plots\n",
    "\n",
    "This is really going to combine two things we already have at our fingertips: space and time. There are various ways we could approach this. The first we'll see is to do a series of *small multiples*, also known as a *facet chart*, and advance time by one unit in each plot. The second is just a heat map in which the two dimensions are space and time.\n",
    "\n",
    "The data we'll be using tell a tragic story. They are counts of deaths that occurred within 28 days of (a known case of) coronavirus organised by the death date. Note that there are various issues with these data and they do not tell the whole story of coronavirus by any means. But this is just an illustrative example.\n",
    "\n",
    "We'll just bring in data for London. They don't come with their own geometry, so our first job is to merge them onto our existing UK local authority geodataframe, which does have a geometry. Fortunately, both data sets have the 'LAD20CD' and 'LAD20NM' columns, which makes this easier than it might have been."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df[\"LAD20CD\"].str.contains(\"E09\")]\n",
    "cv_df = pd.read_parquet(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/data/geo/cv_ldn_deaths.parquet\"\n",
    ")\n",
    "df = df.merge(cv_df, on=[\"LAD20CD\", \"LAD20NM\"], how=\"inner\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will create the small multiple chart. There's quite a lot going on in the below so let's talk through the moving pieces. There's a perpectually uniform heatmap from the **colorcet** package that ranges between a min of 0 and max of 250 as set by `vmin` and `vmax`. To plot every one of the 12 time periods, `plt.subplots` is called with `nrows=4` and `ncols=3` for 12 axes in total. We `.flatten()` the axes to make it easier to iterate over them. We turn the legend off for all but one axis. With all of this set up, we iterate through the unique values of the date column, which is at monthly frequency, and plot each in turn on a separate axis oject. For each, we also add a title that is the month and year.\n",
    "\n",
    "The overall effect is quite powerful: you can really see how deaths peaked around April 2020 and January 2021, having begun to pick up in November, but you can also see the long period of relatively few deaths 28 days after diagnosis during the summer months."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "import colorcet as cc\n",
    "\n",
    "col = \"newDeaths28DaysByDeathDate\"\n",
    "\n",
    "fig, axes = plt.subplots(nrows=4, ncols=3)\n",
    "axes = axes.flatten()\n",
    "legend_choice = [False] * len(axes)\n",
    "legend_choice[2] = True\n",
    "for i, date in enumerate(df[\"date\"].unique()):\n",
    "    df_cut = df[df[\"date\"] == date]\n",
    "    df_cut.plot(\n",
    "        ax=axes[i],\n",
    "        column=col,\n",
    "        legend_kwds={\"label\": \"\", \"shrink\": 1.5},\n",
    "        vmin=0,\n",
    "        vmax=250,\n",
    "        legend=legend_choice[i],\n",
    "        cmap=cc.cm.CET_L19,\n",
    "    )\n",
    "    axes[i].axis(\"off\")\n",
    "    axes[i].set_title(pd.to_datetime(str(date)).strftime(\"%B %Y\"), size=10)\n",
    "plt.suptitle(\"Coronavirus - deaths within 28 days of diagnosis\", size=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to show time and space is using a heatmap. The easiest way to plot a heatmap is to put the data into wide format first. To ensure we have nice labels though, we're going to cast the datetime variable in a month and year in the format 'Jan 20' using `strftime` first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"Time\"] = df.date.apply(lambda x: x.strftime(\"%b \\n %y\"))\n",
    "hmap_df = df.pivot(\"Time\", \"LAD20NM\", \"newDeaths28DaysByDeathDate\").T\n",
    "# Puts the datetimes in the correct order\n",
    "hmap_df = hmap_df[list(df[\"Time\"].unique())]\n",
    "hmap_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the data where we want it, we can call the `ax.imshow` method with a colourmap and set the labels to the index and columns in order to show how, for each London Borough, the number of deaths has changed over time. Note that, unlike the geospatial map, this does not capture the linkages in space very well: but it arguably does a better job of capturing linkages in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "im = ax.imshow(hmap_df.values, cmap=cc.cm.CET_L19)\n",
    "cbar = ax.figure.colorbar(im, ax=ax, aspect=50)\n",
    "ax.set_xticks(np.arange(len(hmap_df.columns)))\n",
    "ax.set_yticks(np.arange(len(hmap_df.index)))\n",
    "# Labels\n",
    "ax.set_xticklabels(hmap_df.columns, rotation=0, fontsize=8)\n",
    "ax.set_yticklabels(hmap_df.index, fontsize=8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using basemaps\n",
    "\n",
    "All of the examples we've seen have just been lines around coloured polygons. This is, fortunately, not how most maps look. Instead, they have lots of detail or other features. Sometimes you want to incorporate other features into your map, or just give some context to where things are beyond the boundaries of your polygons. Enter different basemaps, ie different backgrounds for you to draw the geospatial data on top of.\n",
    "\n",
    "There are a few options to do this. Let's begin with **geoplot**, the package from the examples above. It has projection called `WebMercator` that is a 'real' map. It's easiest to illustrate with an example; let's use Wales again. But the whole point is that we can combine this map with other things, so we'll do a few points to show where some places are: LLandudno, Cardiff, Rhyll, Newport, and St Davids. St Davids is the smallest city in Britain, and we'll make the points proportional to population, so you may have to squint to see it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wales_places = {\n",
    "    \"name\": [\"Cardiff\", \"Llandudno\", \"Rhyl\", \"Newport\", \"St Davids\"],\n",
    "    \"lat\": [51.481583, 53.3241, 53.3191, 51.5842, 51.8812],\n",
    "    \"lon\": [-3.179090, -3.8276, -3.4916, -2.9977, -5.2660],\n",
    "    \"pop\": [335145, 20701, 25149, 145700, 1600],\n",
    "}\n",
    "wdf = pd.DataFrame(wales_places)\n",
    "gwdf = gpd.GeoDataFrame(wdf, geometry=gpd.points_from_xy(wdf.lon, wdf.lat))\n",
    "gwdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create basemap\n",
    "ax = gplt.webmap(df_wales, projection=gcrs.WebMercator())\n",
    "# Add points for places\n",
    "gplt.pointplot(gwdf, ax=ax, hue=\"name\", legend=True, sizes=gwdf[\"pop\"] / 300);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contextily\n",
    "\n",
    "**geoplot** has some basic functionality for basemaps but the [**contextlily**](https://contextily.readthedocs.io/en/latest/intro_guide.html) package provides a whole lot more flexibility and is solely focused on different 'context tiles' for your map. It's also designed to work with **geopandas** (it's built by the same people) and the syntax is quite similar to what you've seen already.\n",
    "\n",
    "Let's see an example of it in action. We'll ask for a bounding box around a place we're interested in, London's Canary Wharf. The option that brings in different basemaps is the `source=` keyword argument. There are a range of sources of basemaps but the real magic is that, given a location, **contextily** *downloads a background map for you*. Pretty cool. There is a full list of providers [here](https://contextily.readthedocs.io/en/latest/providers_deepdive.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import contextily as cx\n",
    "\n",
    "west, south, east, north = (-0.030251, 51.499019, 0.002017, 51.509511)\n",
    "cw_img, cw_ext = cx.bounds2img(\n",
    "    west, south, east, north, ll=True, source=cx.providers.Stamen.Toner\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the map gets downloaded when you make this request, it takes a long time the first time it's run. But the map is locally cached so that, if you call the same function again, it will be much faster the second time. Let's plot the map out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.axis(\"off\")\n",
    "ax.imshow(cw_img, extent=cw_ext);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also get a map through a text search, though be warned you may get another place with the same name in a different region. You can do this type of map-grab using `cx.Place(\"...\", source=...)`.\n",
    "\n",
    "Okay, it's great to download a map but what about combining it with some useful info? Well, we can do that too. Let's use **osmnx** to get some data on coffee shops in an area more tightly focused around Canary Wharf. We'll then pop these onto a map of the area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import osmnx as ox\n",
    "\n",
    "coffee_shops = ox.geometries_from_place(\n",
    "    \"Canary Wharf\", tags={\"amenity\": \"cafe\"}, buffer_dist=300\n",
    ")\n",
    "coffee_shops = coffee_shops.to_crs(\"EPSG:3857\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(dpi=150)\n",
    "coffee_shops.plot(ax=ax, markersize=80, color=\"darkviolet\", edgecolor=\"k\", marker=\"X\")\n",
    "ax.axis(\"off\")\n",
    "cx.add_basemap(\n",
    "    ax, crs=coffee_shops.crs.to_string(), source=cx.providers.OpenStreetMap.Mapnik\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also add boundary and polygon objects to **contextily** basemaps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive maps\n",
    "\n",
    "We'll use [**folium**](https://python-visualization.github.io/folium/index.html), a wrapper for the leaflet javascript library, to create interactive maps that we can layer information on top of. The library has a number of built-in tilesets from OpenStreetMap, Mapbox, and Stamen, and supports custom tilesets too. If it has a disadvantage, it's that it doesn't play that nicely with **geopandas**. It isn't easy to create an interactive choropleth using a **geopandas** dataframe, for example. However, choropleths can be added as a layer from a URL that points to a geoJSON file.\n",
    "\n",
    "In the simple example below, we'll do two things: create a basemap showing the City of London and add some markers to it that show new information on a mouse hover and when clicked. Let's put markers down for the places to get coffee in the serenity of a City church."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "# create data for markers\n",
    "caff_df = pd.DataFrame(\n",
    "    {\n",
    "        \"name\": [\"Host Coffee\", \"Cosy Coffee Corner\", \"The Wren Coffee\"],\n",
    "        \"loc\": [[51.5128, -0.0933], [51.5128, -0.0882], [51.512106, -0.096870]],\n",
    "    }\n",
    ")\n",
    "# create hover over marker msg\n",
    "tooltip = \"Click here.\"\n",
    "\n",
    "# create map\n",
    "m = folium.Map(location=[51.5128, -0.0933], zoom_start=16)\n",
    "# add markers to map\n",
    "for index, row in caff_df.iterrows():\n",
    "    folium.Marker(row[\"loc\"], popup=row[\"name\"], tooltip=tooltip).add_to(m)\n",
    "# show map\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Folium can also be combined with a choropleth; see the [documentation](https://python-visualization.github.io/folium/quickstart.html#) for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Static maps\n",
    "\n",
    "This is an honourable mention to [**py-staticmaps**](https://github.com/flopp/py-staticmaps), which provides the kind of tilesets we've seen already with **contextlily** and **folium** but includes easy methods to add:\n",
    "\n",
    "- markers\n",
    "- image (PNG) markers\n",
    "- geodesic lines, i.e. the shortest path between two points given the geometry. These are *very* important in general relativity, where they appear as solutions to the geodesic equation  ${\\displaystyle {d^{2}x^{\\mu } \\over ds^{2}}+\\Gamma ^{\\mu }{}_{\\alpha \\beta }{dx^{\\alpha } \\over ds}{dx^{\\beta } \\over ds}=0\\ }$\n",
    "- geodesic circles\n",
    "- polygons; and\n",
    "- GPX Tracks, i.e. the paths traced out when you record the route of a run using GPS.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "If you know how to :\n",
    "\n",
    "- ✅  plot geographic data on a map;\n",
    "- ✅  plot choropleths of different kinds on maps;\n",
    "- ✅  create cartographs;\n",
    "- ✅  create quadtrees and when you might use one over a choropleth;\n",
    "- ✅  create geographic kernel density estimate plots;\n",
    "- ✅  use different basemaps in your geospatial visualisations;\n",
    "- ✅  show time and space dimensions on plots; and\n",
    "- ✅  produce interactive geospatial maps\n",
    "\n",
    "then you are well on your way to becoming a geospatial visualisation pro!\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "interpreter": {
   "hash": "671f4d32165728098ed6607f79d86bfe6b725b450a30021a55936f1af379a247"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 64-bit ('codeforecon': conda)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}