geo-intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro to Geo-Spatial Analysis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn the basics of geospatial analysis with code.\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.) If you're not familiar with Docker and you're installing everything locally, a strong recommendation is to use conda. There's a brief guide on installing packages in the Chapter on {ref}`code-preliminaries`.\n",
    "\n",
    "There are two types of spatial data geographic information systems (GIS): vector and raster. Vector spatial data are made up of vertices and paths. These are represented computationally by points (zero-dimensional objects, if you like), lines (1D objects), and polygons (2D objects, aka an area). Vector data are analogous to vector image formats, like .eps and .pdf, they are smooth and well-defined regardless of the level of zoom you look at them (on a map). Raster data are composed of pixels, also known as grid cells. Each cell has a different value and the raster map is overlaid on a map. Raster data are better for continuous variables such as temperature or rainfall, while vector data are better suited to political boundaries and features like roads. In this book, we'll only cover vector geospatial data analysis (if you need to work with rasters, check out [rasterio](https://rasterio.readthedocs.io/en/latest/)).\n",
    "\n",
    "Common file formats for vector data include Shapefile (.shp...), GeoJSON/JSON, KML and GeoPackage.\n",
    "\n",
    "### Coordinate Reference Systems\n",
    "\n",
    "A Coordinate Reference System (CRS) associates numerical coordinates with a position on the surface of the Earth. Latitude and longitude are an example (and have units of degrees) but there are many CRSs. They are important because they define how your data will look on a map. Think about it like this: in the case of the usual charts you plot, you usually take it as given that you are working in a space that is 2D with an X- and Y-basis vector. But the most basic object in mapping is a *sphere*, which is fundamentally different to a 2D plane. This means you have to choose whether to show part of a globe or all of a 2D representation of a globe (a *projection*), which inevitably introduces distortions.\n",
    "\n",
    "The type main classes of CRS are geographic or projection. A geographic CRS uses a 3D model of the Earth and 'drapes' the data over it.\n",
    "\n",
    "A projected CRS has the geographic CRS plus a map projection that has co-ordinates on a 2D plane. It is well known that there is no distortion-free projection from a globe to a plane; you cannot preserve areas, shapes, distances, and angles simultaneously. Different map projections introduce different distortions, as lovingly shown in this [xkcd](https://xkcd.com/977/) cartoon.\n",
    "\n",
    "One example of a map projection is the Mercator projection, which is a *conformal mapping*, i.e. a mapping that locally preserves angles, but not necessarily lengths. In fact, Mercator distorts areas, especially the further away an area is from the equator. Some projections are better for some purposes than others.\n",
    "\n",
    "![XKCD Bad Map Projection: South America](https://imgs.xkcd.com/comics/bad_map_projection_south_america.png)\n",
    "\n",
    "XKCD: Bad Map Projection: South America\n",
    "\n",
    "Some analysis tools expect geospatial data to be in a projected CRS-this is true of the main package we'll use, **geopandas**. This is usually not a problem for economic data; it's rare that the curvature of the Earth becomes a factor (though distances might in some rare situations). Most spatial libraries expect that all of the data that are being combined be in the same CRS.\n",
    "\n",
    "CRSs are usually referenced by an [EPSG code](https://epsg.io/). Two widely used CRSs are WGS84 (aka EPSG: 4326) which provides a good representaiton of most places on Earth and NAD83 (aka EPSG: 4269), which provides a good representation of the USA.\n",
    "\n",
    "Why are we banging on about this? Because maps and geometry objects come in different CRS and it's worth being aware of that so that you can ensure they are in the same format, and that you have the right CRS for your purposes.\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"
   ]
  },
  {
   "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 chapter, some bespoke settings:\n",
    "plt.style.use(\n",
    "    {\n",
    "        \"axes.autolimit_mode\": \"data\",\n",
    "        \"patch.linewidth\": 0.2,\n",
    "        \"figure.figsize\": (10, 8),\n",
    "        \"figure.dpi\": 125,\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geopandas dataframes\n",
    "\n",
    "Quite literally, **GeoPandas** is a combination of geo and pandas so the good news is that everything you know about using **pandas** dataframes can be re-used here for geospatial data. The geo part adds functionality for geo-spatial data.\n",
    "\n",
    "### Input and Output\n",
    "\n",
    "So, first, we need some geo-spatial data to analyse. There are several different file formats for geo-spatial data, such as Shapefile (.shp), GeoJSON/JSON, KML, and GeoPackage.\n",
    "\n",
    "We'll use a Shapefile of the countries of the world from [Natural Earth](https://www.naturalearthdata.com/downloads/50m-cultural-vectors/50m-admin-0-countries-2/). It comes as a zip file; unzip it and one of the files ends in .shp, which is the one we load with **geopandas**.\n",
    "\n",
    "Let's load the data and look at the first few rows:"
   ]
  },
  {
   "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/world\n",
    "# and save them in the path 'data/geo/world' (path relative to where you're running the code)\n",
    "df = gpd.read_file(Path(\"data/geo/world/ne_50m_admin_0_countries.shp\"))\n",
    "df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's a lot of info here, but a lot of it is different labelling. The dataset has one country per row.\n",
    "\n",
    "To save the data frame, use `df.to_file(\"output_name.shp\")` for shapefiles. For other output formats use, for example, `df.to_file(\"output_name.geojson\", driver=\"GeoJSON\")`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basics\n",
    "\n",
    "Let's see what we get when we call the humble plot function on the data we already read in!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot(color=\"green\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think it's glorious just how easy this is to use. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because **geopandas** builds on **pandas**, we can do all of the usual pandas-like operations such as filtering based on values. Here's how to filter to an individual country (one with a recognisable shape!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df[\"ADMIN\"] == \"Italy\"].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that it is the last column in the dataframe, the `geometry` column, that makes this sorcery possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with Coordinate Reference Systems\n",
    "\n",
    "We can check what the CRS of this entire **geopandas** dataframe is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can switch between CRS using the `to_crs` function. Let's see the world map we plotted earlier using the WGS84 projection and also using some other projectsion, including the dreaded Meractor projection. Mercator looks completely ridiculous if we don't drop Antarctica first, so let's get shot of it for all the projections. (You can find a list of projections [here](https://proj.org/operations/projections/index.html).)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclude_admins = [\"Antarctica\", \"French Southern and Antarctic Lands\"]\n",
    "proj_names = [\"WGS 84\", \"Mercator\", \"Winkel Tripel\", \"Quartic Authalic\"]\n",
    "crs_names = [\"EPSG:4326\", \"EPSG:3395\", \"+proj=wintri\", \"+proj=qua_aut\"]\n",
    "\n",
    "world_no_antrtc = df[~df[\"ADMIN\"].isin(exclude_admins)]\n",
    "\n",
    "fig, axes = plt.subplots(\n",
    "    nrows=2, ncols=2, figsize=(12, 6), dpi=150, sharex=False, sharey=False\n",
    ")\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    world_no_antrtc.to_crs(crs_names[i]).plot(ax=ax, color=\"#b2f3b2\", linewidth=0.1)\n",
    "    ax.set_title(proj_names[i])\n",
    "    ax.axis(\"off\")\n",
    "plt.tight_layout()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manipulating Space\n",
    "\n",
    "### Basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look more closely at the objects that encode the shapes of the countries on our world map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[df[\"ADMIN\"] == \"Italy\", \"geometry\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object is a multipolygon. Remember that we have points (0D), lines (1D), and polygons (aka areas, 2D) that we can embed in a projection. A line is at least 2 vertices that are connected; a polygon is the area contained within 3 or more vertices. Multipolygons are the union of two or more non-contiguous polygons: in this case, the Italian mainland, Sicily, and Sardinia.\n",
    "\n",
    "The `plot` function works just as happily if our basic objects are points rather than polygons though. In the below example, we'll grab the centroid (the spatial midpoint) of each country as a point and plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"centroid\"] = df.centroid\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.axis(\"off\")\n",
    "df.boundary.plot(ax=ax, lw=0.2, color=\"k\")\n",
    "df[\"centroid\"].plot(ax=ax, color=\"red\", markersize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore those basic building blocks a bit more. A point at position (1, 2) is defined as follows (**shapely** is used by **geopandas**):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "\n",
    "point = Point(1, 2)\n",
    "point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A point doesn't have much other than a single position in 2D space. But lines have length, and polygons have area.\n",
    "\n",
    "There are different kinds of lines but the simplest is the `LineString` which can be constructed from a sequence of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import LineString\n",
    "\n",
    "line = LineString(\n",
    "    [Point(0.0, 1.0), Point(2.0, 2.0), Point(2.5, 5.0), Point(4, 5), Point(4, 0)]\n",
    ")\n",
    "print(f\"Length of line is {line.length:.2f}\")\n",
    "line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already saw Polygons in the shape of Italy. But here's a much simpler one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "\n",
    "poly = Polygon([(0, 0), (1, 1), (1, 0)])\n",
    "print(f\"The area of the poly is {poly.area}\")\n",
    "poly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial Operations\n",
    "\n",
    "We've seen the basic builidng blocks of geometries: points, lines, and polygons. We've also seen how to retireve some basic properties of these such as length, area, and centroid. In this section, we'll see some slightly more advanced spatial operations that you can perform with **geopandas**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Point-in-polygon\n",
    "\n",
    "This does pretty much what you'd expect! It's useful in practice because we might want to know if a particular place falls within one area or another. As a test, let's see if the centroid for Germany actually falls within the Germany polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_row = df.loc[df[\"ADMIN\"] == \"Germany\", :]\n",
    "df_row[\"geometry\"].iloc[0].contains(df_row[\"centroid\"].iloc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But be careful with this! Countries can have complex multi-polygon geometries for which a centroid is not guaranteed to be within any of the polygons. France is a great example as French Guiana is so far away that it pulls the centroid just out of the mainland:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "df_row = df.loc[df[\"ADMIN\"] == \"France\", :]\n",
    "df_row[\"geometry\"].plot(ax=ax)\n",
    "df_row[\"centroid\"].plot(ax=ax, color=\"red\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Buffers\n",
    "\n",
    "Buffers are just an area drawn around a particular geometry, given a specific radious. They have a great practical use in computing catchment areas and so on. To create a buffer around the geometry you're currently using, use the `df.buffer(number)` command to return a column with the buffered version of your geometry in. Be aware that the units of the CRS you're using will determine the units of the buffer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial set operations\n",
    "\n",
    "More complex spatial manipulation can be achieved through spatial set operations. The figure below shows some examples for polygons (but the same principles apply for lines and polygons too):\n",
    "\n",
    "<img src=\"https://geopandas.org/_images/overlay_operations.png\" alt=\"Spatial operations\"></a>\n",
    "\n",
    "Different set operations with two polygons. Source: QGIS documentation.\n",
    "\n",
    "In addition to these, there are 'crosses' (for two lines) and 'touches'. You can find information about all of the set operations that are available [here](https://geopandas.org/set_operations.html).\n",
    "\n",
    "To demonstrate one of these, let's see if we can pick out a few regions that are *intersected* by a river. We'll try out the river Trent, which passes through England, and we'll see which *Local Authority Districts* (LADs) it passes through. First, we load the UK data, which is from the [ONS Open Geography Portal](https://geoportal.statistics.gov.uk/datasets/local-authority-districts-may-2020-boundaries-uk-buc):"
   ]
  },
  {
   "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",
    "dfuk = gpd.read_file(\n",
    "    Path(\"data/geo/uk_lad/Local_Authority_Districts__May_2020__UK_BUC.shp\")\n",
    ")\n",
    "dfuk.plot(color=\"dimgrey\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we bring in data on major rivers and pass it through the same projection as our map of England is using:"
   ]
  },
  {
   "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/rivers\n",
    "# and save them in the path 'data/geo/rivers' (path relative to where you're running the code)\n",
    "river = gpd.read_file(Path(\"data/geo/rivers/rivers.shp\"))\n",
    "river = river.to_crs(dfuk.crs)\n",
    "river.plot(color=\"mediumblue\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now combine these for a view on how these rivers goes through the UK. We will use `buffer` to modify how easy it is to visualise the river by asking for everything within 1.5 km from the river."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "ax.set_frame_on(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "dfuk.plot(color=\"grey\", ax=ax)\n",
    "river.buffer(1.5e3).plot(ax=ax, color=\"lightblue\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's turn our focus to one river: the Trent. We will ask which of the regions are intersected by that river and plot only those parts. We'll also add some annotations for the names of the regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict river to just one\n",
    "river_name = \"Trent\"\n",
    "river = river[river[\"name\"] == river_name]\n",
    "\n",
    "# Find a subset of the regions that is intersected by the river by creating new True/False column\n",
    "dfuk[\"river\"] = dfuk[\"geometry\"].apply(\n",
    "    lambda x: river[\"geometry\"].buffer(1e3).intersects(x)\n",
    ")\n",
    "\n",
    "# The rest of the code is just related to plotting:\n",
    "\n",
    "# Create a cut of dfuk for convenience\n",
    "df_th = dfuk[dfuk[\"river\"]].copy()\n",
    "\n",
    "# Get a representative point for each region to annotate\n",
    "df_th[\"coords\"] = (\n",
    "    df_th[\"geometry\"].representative_point().apply(lambda x: x.coords[:][0])\n",
    ")\n",
    "\n",
    "# Plotting bits\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.set_frame_on(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "df_th.plot(color=\"grey\", ax=ax, edgecolor=\"0.6\")\n",
    "river.plot(color=\"lightblue\", ax=ax)\n",
    "\n",
    "# Add text annotation to the largest polygons\n",
    "for idx, row in df_th.iterrows():\n",
    "    if row[\"geometry\"].area > np.quantile(df_th.area, q=0.4):\n",
    "        ax.annotate(\n",
    "            text=row[\"LAD20NM\"],\n",
    "            xy=row[\"coords\"],\n",
    "            horizontalalignment=\"center\",\n",
    "            weight=\"bold\",\n",
    "            fontsize=10,\n",
    "            color=\"black\",\n",
    "        )\n",
    "ax.set_title(\n",
    "    f\"Local Authority Districts that are intersected by the {river_name}\",\n",
    "    loc=\"left\",\n",
    "    fontsize=20,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the way we set this up to use `river_name` as an input variable means we could oh so easily wrap everything up in a function that did this for other rivers. Oh go on then, let's see how that works *because* it shows how scalable operations are once you do them in *code*:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def river_intersect_plot(river_name):\n",
    "    \"\"\"\n",
    "    Given the name of a river, shows a plot of the LADs that it intersects.\n",
    "    \"\"\"\n",
    "    river = gpd.read_file(os.path.join(\"data\", \"geo\", \"rivers\", \"rivers.shp\"))\n",
    "    river = river[river[\"name\"] == river_name]\n",
    "    river = river.to_crs(dfuk.crs)\n",
    "    dfuk[\"river\"] = dfuk[\"geometry\"].apply(\n",
    "        lambda x: river[\"geometry\"].buffer(2e3).intersects(x)\n",
    "    )\n",
    "    df_th = dfuk[dfuk[\"river\"]].copy()\n",
    "    df_th[\"coords\"] = (\n",
    "        df_th[\"geometry\"].representative_point().apply(lambda x: x.coords[:][0])\n",
    "    )\n",
    "    fig, ax = plt.subplots(figsize=(10, 10))\n",
    "    ax.set_frame_on(False)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    df_th.plot(color=\"grey\", ax=ax, edgecolor=\"0.6\")\n",
    "    river.plot(color=\"lightblue\", ax=ax)\n",
    "    for idx, row in df_th.iterrows():\n",
    "        if row[\"geometry\"].area > np.quantile(df_th.area, q=0.6):\n",
    "            ax.annotate(\n",
    "                text=row[\"LAD20NM\"],\n",
    "                xy=row[\"coords\"],\n",
    "                horizontalalignment=\"center\",\n",
    "                weight=\"bold\",\n",
    "                fontsize=10,\n",
    "                color=\"black\",\n",
    "            )\n",
    "    ax.set_title(\n",
    "        f\"Local Authority Districts that are intersected by the {river_name}\",\n",
    "        loc=\"left\",\n",
    "        fontsize=20,\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our function defined, we can do the whole thing for a completely different river:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "river_intersect_plot(\"Thames\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing distances\n",
    "\n",
    "Computing the distances between two points `a` and `b` (in the same CRS) is achieved using the `a.distance(b)` command. \n",
    "\n",
    "As an example, let's now find out how far it is between two regions of interest in our cut of regions around the Trent. We'll pick the East Riding of Yorkshire and Stafford, and then plot them on the map in different colours. Because `geopandas` is so flexible, and there could be multiple geometries with the names 'East Riding of Yorkshire' and 'Stafford' we have to not only select the names we want but, in case there are multiple rows, specify which row. If you just want the first you can use `iloc[0]`, which is what happens in the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the rows we're interested in out of the dataframe:\n",
    "name_a = \"East Riding of Yorkshire\"\n",
    "name_b = \"Stafford\"\n",
    "# This selects the *all rows* that match these conditions, which is why we have to use .iloc[0] thereafter\n",
    "# to make sure we're only passing a single row.\n",
    "place_a = df_th.loc[df_th[\"LAD20NM\"] == \"East Riding of Yorkshire\", :]\n",
    "place_b = df_th.loc[df_th[\"LAD20NM\"] == \"Stafford\", :]\n",
    "# Compute the distance using representative points\n",
    "dist_a_b = (\n",
    "    place_a[\"geometry\"]\n",
    "    .representative_point()\n",
    "    .iloc[0]\n",
    "    .distance(place_b[\"geometry\"].representative_point().iloc[0])\n",
    ")\n",
    "\n",
    "# Plot the map\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.set_frame_on(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "df_th.plot(color=\"grey\", ax=ax, edgecolor=\"0.6\")\n",
    "river.plot(color=\"lightblue\", ax=ax)\n",
    "place_a.plot(ax=ax, color=\"red\")\n",
    "place_b.plot(ax=ax, color=\"green\")\n",
    "for i, place in enumerate([place_a, place_b]):\n",
    "    ax.annotate(\n",
    "        text=place[\"LAD20NM\"].iloc[0],\n",
    "        xy=place[\"coords\"].iloc[0],\n",
    "        horizontalalignment=\"center\",\n",
    "        weight=\"bold\",\n",
    "        fontsize=15,\n",
    "        color=\"black\",\n",
    "    )\n",
    "# Uncomment below to add a connecting line\n",
    "## Create a line between two rep points in a and b\n",
    "# connector = LineString([place_a['geometry'].representative_point().iloc[0],\n",
    "#                        place_b['geometry'].representative_point().iloc[0]])\n",
    "## Convert it to a geopandas dataframe for easy plotting\n",
    "# gpd.GeoDataFrame([connector],columns=['line'], geometry='line').plot(ax=ax, linestyle='-.', color='black')\n",
    "ax.set_title(\n",
    "    f\"Distance between {name_a} and {name_b} is {dist_a_b/1e3:1.0f} km.\", loc=\"left\"\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neighbours\n",
    "\n",
    "It's really useful to have a way to determine all of the adjacent geometries to a given geometry. The `touch` geospatial operation is the best for this *but* do be careful because the level of the detail in polygons varies for different maps and you may find that they sometimes overlap rather than simply touch. So, in the example of finding all of the neighbours for all of the local authority districts in the UK below, not only do we include those LADs that touch, we also include those that overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Empty list of lists\n",
    "neighbours_all = [[] for _ in range(len(dfuk))]\n",
    "# Iterate over rows\n",
    "for index, row in dfuk.iterrows():\n",
    "    # Find any touches of this row\n",
    "    neighbours = dfuk.loc[dfuk[\"geometry\"].touches(row[\"geometry\"]), \"LAD20NM\"].tolist()\n",
    "    # Find any overlaps of this row\n",
    "    overlap = dfuk.loc[dfuk[\"geometry\"].overlaps(row[\"geometry\"]), \"LAD20NM\"].tolist()\n",
    "    # Take the union of touches and overlaps\n",
    "    neighbours = np.union1d(neighbours, overlap)\n",
    "    # Store the result\n",
    "    neighbours_all[index] = neighbours\n",
    "\n",
    "# Put result back into geopandas dataframe (RHS is a list of lists of strings)\n",
    "dfuk[\"neighbours\"] = neighbours_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we should have an extra column that tells us which regions are adjacent to which by name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfuk.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot an example to be sure we've done it all right!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name_of_region = \"High Peak\"\n",
    "# Set colors now for ease in plotting\n",
    "dfuk[\"color\"] = \"#b2df8a\"\n",
    "dfuk.loc[dfuk[\"LAD20NM\"] == name_of_region, \"color\"] = \"#33a02c\"\n",
    "# Get representative points to use for text annotation. Note that we want xy coords, hence the fussy use of .apply\n",
    "dfuk[\"coords\"] = dfuk[\"geometry\"].representative_point().apply(lambda x: x.coords[:][0])\n",
    "# Get a list of the neighbours--again, this comes out as a dataframe row so we need to call iloc[0] to get first entry\n",
    "names_neighbours = list(\n",
    "    dfuk.loc[dfuk[\"LAD20NM\"] == name_of_region, \"neighbours\"].iloc[0]\n",
    ")\n",
    "# Cut to only those areas of interest\n",
    "df_cut = dfuk.loc[dfuk[\"LAD20NM\"].isin(names_neighbours + [name_of_region])]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "df_cut.plot(ax=ax, color=df_cut[\"color\"], edgecolor=\"k\")\n",
    "ax.set_frame_on(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "# Add annotation\n",
    "for idx, row in df_cut.iterrows():\n",
    "    ax.annotate(\n",
    "        text=row[\"LAD20NM\"],\n",
    "        xy=row[\"coords\"],\n",
    "        horizontalalignment=\"center\",\n",
    "        weight=\"bold\",\n",
    "        fontsize=15,\n",
    "        color=\"black\",\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hooray! It looks like it worked. And notice how the only variable we need to change to get a different map is `name_of_region`, so this could easily be wrapped up in a function to re-use again. \n",
    "\n",
    "For an exercise, why don't you see if you can figure out what's so special about Powys and the Isle of Wight."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aggregation and Dissolve\n",
    "\n",
    "### Dissolve\n",
    "\n",
    "A useful common map operation is to aggregate by spatial region: i.e. to eliminate some of the granularity of a map by combining smaller regions into larger regions. This is called a 'dissolve' because it dissolves the boundaries between regions to form larger ones.\n",
    "\n",
    "To demonstrate it, let's take the many UK local authority districts and `dissolve` them into the UK countries. We can do this using the `LAD20CD` column because the code used for it always begins with a letter that represents the country: E for England, W for Wales, N for Northern Ireland, and S for Scotland. So to do the dissolve, we first create a new column just like we would do for a `groupby` in a regular **pandas** dataframe. But instead of using `groupby`, we `dissolve` by the new column. This will give us a new map that is made up of only 4 polygons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfuk[\"country\"] = dfuk[\"LAD20CD\"].apply(lambda x: str(x)[0])\n",
    "dfuk.iloc[:, [1, 2, -1]].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_uk = dfuk.dissolve(by=\"country\")\n",
    "country_uk.plot(edgecolor=\"k\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default way that columns other than the geometry are combined is by taking the first row of each kind, but dissolve accepts and `agg_func=` argument that you can also pass last, min, max, sum, mean, and median."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregation\n",
    "\n",
    "Because **geopandas** builds on **pandas**, all the usual operations (i.e. `groupby`) that you might do on a **pandas** dataframe also work for non-geometry columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging data\n",
    "\n",
    "### Attribute Joins\n",
    "\n",
    "In an attribute join, a GeoSeries or GeoDataFrame is combined with a regular pandas Series or DataFrame based on a common variable. This is analogous to normal merging or joining in pandas.\n",
    "\n",
    "### Spatial Joins\n",
    "\n",
    "Spatial joins are more complex because we take two different sets of geometries and say how to combine them. The key to this is i) ensuring that the datasets have the same CRS and ii) using `sjoin` to perform the spatial join. We'll show how to do a spatial join that is an inner intersection using `op=intersects`. \n",
    "\n",
    "We'll take some data on, and a map of, the world, which has per country geometries in plus columns for, eg, gdp per capita, and we'll merge it with a dataframe with the names and positions of cities in. Our intersects spatial join will combine each city (a lat/lon point) with the shapes of countries and determine which city goes in which country (even if it's on the boundary).\n",
    "\n",
    "This is just one example of a spatial join: the other possible arguments to `op` are `within` and `contains`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab a world map that has extra data in\n",
    "world = gpd.read_file(gpd.datasets.get_path(\"naturalearth_lowres\"))\n",
    "world.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab data on cities\n",
    "cities = gpd.read_file(gpd.datasets.get_path(\"naturalearth_cities\"))\n",
    "# Force cities and world to share the world CRS:\n",
    "cities = cities.to_crs(world.crs)\n",
    "cities.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cities_with_country = gpd.sjoin(cities, world, how=\"inner\", op=\"intersects\")\n",
    "cities_with_country.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have our cities, and their points, merged into a dataset that also gives the continent and other data that was originally only in our `world` dataframe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geocoding\n",
    "\n",
    "Geocoding is the process of taking an address or the name of a place and retrieving its position in a geographic co-ordinate system. There are many entities you might want to do this for, such as cities, landmarks, geographic features like rivers, and postal code. Reverse geocoding takes a latitute and longitude and works backwards to work out the name of the object that's there.\n",
    "\n",
    "Geocoding as a service is offered by several providers, including Geocodio, Google’s geocode API service, Bing, OpenStreetMap, and more. Some of these charge, especially for volume. \n",
    "\n",
    "There are two parts of geocoding: one is the interface that you use to issue the command (the front-end, the bit you see), and another is that does the actual geocoding (the back-end, an API you don't see). The back-end bit is what you sometimes have to pay for as a service if you're geocoding at scale.\n",
    "\n",
    "We'll look at a geocoding package that comes in two different forms: [**Geopy**](https://geopy.readthedocs.io/en/stable/), which is a Python package that provides a front end for a large number of geocoding APIs, including OpenStreetMap, Bing, Google, ArcGIS, and more. You may also want to explore the [**geocoder**](https://geocoder.readthedocs.io/) Python package as an alternative.\n",
    "\n",
    "### Geopy stand-alone\n",
    "\n",
    "Below is an example of using **Geopy** as a stand-alone package. We’ll use the OpenStreetMap API to do the geocoding, which is called `Nominatim`. It’s important to note that this API has some fair usage conditions including a maximum of 1 request per second, that you provide an informative ‘user agent’ parameter, and that you clearly display attribution (thank you OpenStreetMap!). For bulk geocoding, you may need to pay a fee to a provider.\n",
    "\n",
    "Let's see an example retrieving the full address of a landmark:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you don't have it, install geopy using 'pip install geopy'\n",
    "from geopy.geocoders import Nominatim\n",
    "\n",
    "# Create a geolocator using Open Street Map (aka Nominatim)\n",
    "# Use your own user agent identifier here\n",
    "geolocator = Nominatim(user_agent=\"codeforecon\")\n",
    "\n",
    "# Pass an address to retrieve full location information:\n",
    "location = geolocator.geocode(\"Bank of England\")\n",
    "\n",
    "location.address"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The location parameter also contains a latitude and longitude, which is useful for further geospatial operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(location.latitude, location.longitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's an example of using the same API to do a reverse geocode from a given lat/lon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also reverse geocode from a lat and lon:\n",
    "scnd_location = geolocator.reverse(\"51.529969, -0.127688\")\n",
    "\n",
    "scnd_location.address"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geopy built-in to geopandas\n",
    "\n",
    "If you don't want to break out into a separate package, you need not! **geopandas** has its own interface with **geopy**. We'll use it to get the boroughs of ..\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boros_df = gpd.GeoDataFrame(\n",
    "    {\"boro_name\": [\"Staten Island\", \"Queens\", \"Brooklyn\", \"Manhattan\", \"Bronx\"]},\n",
    "    index=range(5),\n",
    ")\n",
    "boros_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we carry out the geocoding step using `gpd.tools.geocode` and pass an entire column of values that we would like to get `POINT` geometries and address for:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boros_df[[\"geometry\", \"address\"]] = gpd.tools.geocode(boros_df.boro_name, provider=\"photon\")\n",
    "boros_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenStreetMap and Osmnx\n",
    "\n",
    "OpenStreetMap (OSM) is a crowd-sourced project that aims to create a free map of the world containing information about streets, buildings, landmarks, landuse, and more. Its users regularly update the information it holds. Think Wikipedia but for a world map. It excels at highly detailed and freely available renderings of streets.\n",
    "\n",
    "[**Osmnx**](https://geoffboeing.com/publications/osmnx-complex-street-networks/) is a package that can retrieve, construct, analyze, and visualize street networks from OpenStreetMap {cite:t}`boeing2017osmnx`. It also does some geocoding on the fly. Let's see it in action and feed in the address we obtained in the geocoding example from earlier; using the `graph_from_point` method will give us the network of streets around a location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import osmnx as ox\n",
    "\n",
    "# Fetch OSM street network from the location\n",
    "graph = ox.graph_from_point(\n",
    "    (location.latitude, location.longitude), dist=1.3e3, network_type=\"all\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned object is of `type` multidigraph, or more precisely `networkx.classes.multidigraph.MultiDiGraph`. This is an object defined in the **networkx** package, and you can find out more about its properties [here](https://networkx.org/documentation/stable/reference/classes/multidigraph.html) but what's important is that, at least here, the network's nodes are where streets meet and its edges are streets.\n",
    "\n",
    "Rather than plot this as nodes and edges, we can make a much nicer plot of the grid of streets using the below snippet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = ox.plot_graph(\n",
    "    graph, bgcolor=\"white\", node_size=0, edge_color=\"grey\", show=False, close=False\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The channel with few streets is the river Thames.\n",
    "\n",
    "We can also use **osmnx** to work out the shortest path between nodes. Let's see a practical example: going from Bank, north of the river Thames, to Borough Market (great for lunch), south of the river. First, we need the lat and long of Borough Market."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lunch_loc = geolocator.geocode(\"Borough Market\")\n",
    "lunch_loc.address"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we find the nearest node to the given locations, and ask **osmnx** to work out a route."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig = ox.nearest_nodes(graph, location.latitude, location.longitude)\n",
    "dest = ox.nearest_nodes(graph, lunch_loc.latitude, lunch_loc.longitude)\n",
    "route = ox.shortest_path(graph, orig, dest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = ox.plot_graph_route(\n",
    "    graph,\n",
    "    route,\n",
    "    route_linewidth=4,\n",
    "    node_size=0,\n",
    "    bgcolor=\"white\",\n",
    "    show=False,\n",
    "    close=False,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, you're in Borough Market but you need to decide where exactly to get your lunch from. OpenStreetMap includes information on many kinds of nearby amenities, including `bar`, `cafe`, and `pub`. We can retrieve this information and put it straight into a **geopandas** dataframe using the `geometries_from_place` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "restaurants = ox.geometries_from_place(\n",
    "    \"Borough Market\", tags={\"amenity\": \"restaurant\"}, buffer_dist=750\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_osm_geoms(osm_df, skip_names=25):\n",
    "    osm_df[\"coords\"] = (\n",
    "        osm_df[\"geometry\"].representative_point().apply(lambda x: x.coords[:][0])\n",
    "    )\n",
    "    fig, ax = ox.plot_graph(\n",
    "        graph, bgcolor=\"white\", node_size=0, edge_color=\"grey\", show=False, close=False\n",
    "    )\n",
    "    osm_df.plot(ax=ax, color=\"blue\", alpha=0.6)\n",
    "    for idx, row in osm_df[::skip_names].iterrows():\n",
    "        ax.annotate(\n",
    "            text=row[\"name\"],\n",
    "            xy=row[\"coords\"],\n",
    "            xytext=(-30, -40),\n",
    "            textcoords=\"offset points\",\n",
    "            horizontalalignment=\"center\",\n",
    "            weight=\"bold\",\n",
    "            fontsize=13,\n",
    "            color=\"black\",\n",
    "            arrowprops=dict(\n",
    "                arrowstyle=\"->\",\n",
    "                color=\"k\",\n",
    "                connectionstyle=\"arc3,rad=0.3\",\n",
    "                lw=1.5,\n",
    "            ),\n",
    "        )\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "plot_osm_geoms(restaurants)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`amenity` is just one type of geometry that can be retrieved too: others are `building`, `landuse`, and `leisure`. You can find a list of them [here](https://wiki.openstreetmap.org/wiki/Map_features). For example, here are the (many) churches within around 1km of the Bank of England:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "churches = ox.geometries_from_place(\n",
    "    \"Bank of England\", tags={\"building\": \"church\"}, buffer_dist=1.3e3\n",
    ")\n",
    "churches = churches.dropna(subset=[\"name\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_osm_geoms(churches, skip_names=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "If you know:\n",
    "\n",
    "- ✅ what a coordinate reference system (crs) is;\n",
    "- ✅ how to open and use geospatial data using **geopandas**;\n",
    "- ✅ how to create quick plots of maps;\n",
    "- ✅ what points, lines, and polygons are in a geospatial context;\n",
    "- ✅ how to perform set operations on points, lines, and polygons;\n",
    "- ✅ how to compute distances, nearest neighbours, and other geospatial properties;\n",
    "- ✅ how to aggregate, merge, and dissolve geospatial data; and\n",
    "- ✅ how to use geocoding to retrieve geospatial data\n",
    "\n",
    "then you have a good grasp of the basics of geospatial analysis!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "interpreter": {
   "hash": "c4570b151692b3082981c89d172815ada9960dee4eb0bedb37dc10c95601d3bd"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}