5SSG2060_Practical_6_Topology_Buffer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5SSG2060 Practical Week 6 :  Spatial Operation & Buffering\n",
    "<a href=\"#This Week's Overview\">This Week's Overview</a>\n",
    "\n",
    "<a href=\"#Learn Outcomes\">Learn Outcomes</a> \n",
    "\n",
    "<a href='#Get prepared'>Get prepared</a>\n",
    "\n",
    "<a href='#Point Data'>Point Data</a>\n",
    "  - <a href='#Aspatial inquiry and spatial visualization'>Aspatial inquiry and spatial visualization<a/>\n",
    "  - <a href='#Measure distances between points'>Measure distances between points<a/>\n",
    "    \n",
    "<a href='#Linestring Data'>Linestring Data</a>\n",
    "  - <a href='#Geo-location and spatial operation'>Geo-location and spatial operation<a/>\n",
    "  - <a href='#Spatial relationships between point and line'>Spatial relationships between point and line<a/>\n",
    "  - <a href='#Buffer'>Buffer<a/>\n",
    "    \n",
    "<a href='#Polygon Data'>Polygon Data</a>\n",
    "\n",
    "<a href='#Spatial relationships and operations'>Spatial relationships and operations<a/>\n",
    "  - <a href='#Intersections'>Intersections<a/>\n",
    "  - <a href='#Difference'>Difference<a/>\n",
    "  - <a href='#Union and Cascaded Union'>Union and Cascaded Union<a/>\n",
    "\n",
    "<a href='#Get stations and railways within City of London'>Get stations and railways within City of London</a>\n",
    "  - <a href='#Clip stations within city of london'>Clip stations within city of london<a/>\n",
    "  - <a href='#Clip railways within City of London'>Clip railways within City of London<a/>\n",
    "  - <a href='#Join the spatial data'>Join the spatial data<a/>\n",
    "\n",
    "<a href='#Task'>Task<a/>\n",
    "- <a href='#Task 1'>Task 1<a/>\n",
    "- <a href='#Task 2'>Task 2<a/>\n",
    "- <a href='#Task 3'>Task 3<a/>\n",
    "- <a href='#Task 4'>Task 4<a/>\n",
    "- <a href='#Task 5'>Task 5<a/>\n",
    "- <a href='#Task 6'>Task 6<a/>\n",
    "- <a href='#Task 7'>Task 7<a/>\n",
    "- <a href='#Task 8'>Task 8<a/>\n",
    "- <a href='#Task 9'>Task 9<a/>\n",
    "- <a href='#Task 10'>Task 10<a/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"This Week's Overview\">This Week's Overview</a>\n",
    "\n",
    "In this practical, besides of the `GeoPandas` introduced in previous weeks, we are going to practice the skills using [Shapely](https://pypi.org/project/Shapely/) in `PyPI`, for manipulation and analysis of planar geometric objects. You may want to explore its functions through [Shapely User Manual](https://shapely.readthedocs.io/en/stable/manual.html) on adding the functionality to work with planar geometric features, including the creation and manipulation of geometries such as points, polylines, and polygons; get familiar with topology on proximity (including adjacency, connectivity and containment) relationships, as well as spatial analysis functions, e.g.intersection, union, dissolve, etc.\n",
    "\n",
    "It is based on the widely deployed [GEOS](https://trac.osgeo.org/geos/) library (the engine of [PostGIS](http://postgis.refractions.net)) and Java Topology Suite  (JTS) (from which GEOS is ported) libraries. GEOS implements the OpenGIS Simple Features for SQL spatial predicate functions and spatial operators. `Shapely` is not concerned with data formats or coordinate systems, but can be readily integrated with other packages. It provides spatial analysis methods that return new geometric objects. To complete the first section on spatial operation, we will use the following datasets: London stations as Point data, London railways as Polyline data and London LSOAs as Polygon data.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"Learn Outcomes\">Learn Outcomes</a>\n",
    "In previous practicals, we had introduced main geometric objects, the spatial operations and queries on them, as well as the ways to measure distance and proximity (produce buffers), and tessellate space (Voronoi diagrams). By the end of this practical, you should be able to\n",
    "- Gain familiarity with core geometric types (Point, LineString, Polygon, and their collection types)\n",
    "- Make spatial and aspatial inqueries\n",
    "- Identify spatial relationships \n",
    "- Conduct spatial operations to visualize multi-layers "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"Get prepared\">Get prepared</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import urllib\n",
    "import zipfile\n",
    "import numpy as np\n",
    "import pysal\n",
    "from pysal.contrib import shapely_ext\n",
    "import scipy.spatial as spatial\n",
    "import geopandas as gpd\n",
    "from geopandas import GeoDataFrame\n",
    "import pandas as pd\n",
    "import shapely.geometry\n",
    "from shapely.geometry import Point\n",
    "import matplotlib as mpl\n",
    "import matplotlib.path as path\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic geometric types\n",
    "------------------------------\n",
    "For Geospatial vector data in Python, they normally have geometric features, and the 3 fundemental geometric types of objects are:\n",
    "- Point data: represents a single point in space.\n",
    "- Polyline data (\"LineString\"): represents a sequence of points that form a line.\n",
    "- Polygon data: represents a filled area.\n",
    "\n",
    "Each of them can be combined in [multi-part geometries](https://shapely.readthedocs.io/en/stable/manual.html#geometric-objects) for extensive overview.\n",
    "\n",
    "## <a id=\"Point Data\">Point Data</a>\n",
    "\n",
    "Our first primitive geometry type is **Point**, by importing from shapely\\'s `geometry` module. Let's start with a simple customized 2D point example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "\n",
    "point_test=Point(1,1)\n",
    "type(point_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We defined ($x$,$y$) for point_test as (1,1), by calling the `type` function, we can get its type as `shapely.geometry.point.Point`. For geospatial data, we normally define its geometry by latitude and longitude coordinates. So let's have a look of the london stations data and save it as a .shp file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the data dir if it doesn't exist\n",
    "if os.path.isdir('data') is not True:\n",
    "    print(\"Creating 'data' directory...\")\n",
    "    os.mkdir('data')\n",
    "\n",
    "# Configure the download\n",
    "url  = 'https://github.com/cusp-london/Spatial-Data-Analysis/blob/master/london-stations.csv?raw=true'\n",
    "path = os.path.join(\"data\",\"london-stations.csv\")\n",
    "\n",
    "# Download\n",
    "r    = urllib.request.urlretrieve(url, path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the csv file using Pandas\n",
    "stations=pd.read_csv('data/london-stations.csv') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get information on data shape and types for each column \n",
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many stations been recorded in London? For spatial analysis, we need further set up the geometry for this data using its coordinates information by calling `Shapely.geometry.Point`, and save it permanently into disk as a shapefile. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up geodataframe, initially with CRS = WGS84 (since that matches the lon and lat co-ordinates)\n",
    "crs = {'init':'epsg:4326'}\n",
    "# define the geometry of points with coordinates\n",
    "geometry = [shapely.geometry.Point(xy) for xy in zip(stations['longitude'], stations['latitude'])]\n",
    "\n",
    "stations = gpd.GeoDataFrame(stations, \n",
    "                          crs = crs, \n",
    "                          geometry = geometry)\n",
    "# convert the .csv file into .shp file \n",
    "stations.to_file(driver='ESRI Shapefile', filename='data/stations.shp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use Geopandas's GeodataFrame to define the essential elements for geospatial data, i.e. crs and geometry. A coordinate reference system (CRS) determines how the 2D coordinates of the geometry objects should be related to actual places on the earth, see https://docs.qgis.org/2.8/en/docs/gentle_gis_introduction/coordinate_reference_systems.html for further information. \n",
    "\n",
    "Now the london stations data had already saved in your data folder as shapefile, and you may want to try open it in QGIS or ArcGIS by data importing. We can get a brief idea of the layout by conducting plot operation as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot map\n",
    "fig, ax = plt.subplots(1,\n",
    "                       figsize = (8,5),\n",
    "                       dpi = 72,\n",
    "                       facecolor = 'white')\n",
    "\n",
    "ax.set_position([0,0,1,1])   # Puts axis to edge of figure\n",
    "ax.set_axis_on()            # Turns axis on\n",
    "# If you want the facecolour applies to axis area as well as bit around the outside\n",
    "# you may choose to turn the axis off\n",
    "ax.get_xaxis().set_visible(False)   # Turns the x axis off so that 'invisible' axis labels don't take up space\n",
    "ax.get_yaxis().set_visible(False)\n",
    "lims = plt.axis('equal')\n",
    "stations.plot(ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Aspatial inquiry and spatial visualization\">Aspatial inquiry and spatial visualization<a/>\n",
    "It is common to conduct aspatial inquiries when dealing with your own data. So let's start with finding the stations within Zone 1:\n",
    "1. Locate the stations within london transportation Zone 1 and get it plotted as map for information.\n",
    "2. Save stations in Zone 1 as a new shapefile.\n",
    "3. Get ideas on the differences between these two shapefiles' point features by calling `Shapely.MultiPoint`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy.where` function can help to define our targeted area \"zone value <=1\" with a new value \"zone1\" under new column \"zonename\", while the remaining stations will be valued as \"out of zone 1\". To highlight the difference, we will use different colors to plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stations['zonename']=np.where(stations['zone']<=1,'zone1','out of zone 1')\n",
    "stations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pointsPalette = {'zone1': 'orange', 'out of zone 1': 'lightblue'}\n",
    "\n",
    "fig, ax = plt.subplots(1,\n",
    "                       figsize = (8,6),\n",
    "                       dpi = 72,\n",
    "                       facecolor = 'white')\n",
    "\n",
    "for ctype, data in stations.groupby('zonename'):\n",
    "    color = pointsPalette[ctype]\n",
    "    label = ctype \n",
    "    data.plot(color=color,\n",
    "              ax=ax,\n",
    "              label=label,\n",
    "              markersize=30)\n",
    "ax.set_position([0,0,1,1])   # Puts axis to edge of figure\n",
    "ax.set_axis_on() \n",
    "ax.get_xaxis().set_visible(False)   \n",
    "ax.get_yaxis().set_visible(False)\n",
    "lims = plt.axis('equal')\n",
    "ax.set(title='Distribution of underground and overground stations in Central London')\n",
    "ax.legend(fontsize=10,\n",
    "          frameon=True,\n",
    "          loc=(.8, .8),\n",
    "          title=\"LEGEND\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 1\">Task 1<a/> \n",
    "Can you get the number of stations within Zone 1? You should get an answer around 60. Then please save these points into your directory as `stations_zone1.shp`, with CRS at EPSG: 4326."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, you may want to repeat these steps to get stations within Zone 3 or out of Zone 3. Feel free to have a try when you complete this whole practical. But now, let's continue to explore the newly generated data: list the coordinates information for these 60 stations in Zone 1. There are 2 ways to get the list, quite similar to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check output format for column geometry \n",
    "stations_zone1.geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the geometry takes format as POINT (longitude, latitude), we can use it to draw a bounding polygon for stations within zone 1 by calling Shapely's `MultiPoint`. The methods could be either `convex hull` or `envelope` as introduced in previous weeks. Try to explore both datasets and methods, and identify their differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import MultiPoint\n",
    "# the geometry for stations within zone 1 is a list of (longitude, latitude) tuples\n",
    "# the convex hull of this set of 60 points is like an elastic band stretched around the geometry\n",
    "zone1hull = MultiPoint(geometry_zone1).convex_hull\n",
    "zone1hull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The envelope is the point or smallest rectangular polygon (with sides parallel to the coordinate axes) \n",
    "# that contains the geometry\n",
    "zone1enve = MultiPoint(geometry_zone1).envelope\n",
    "zone1enve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinates of the envelope polygon above could be derived by calling bounding box `bounds` from `MultiPoint`. The bounding box is a (minx, miny, maxx, maxy) tuple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zone1bounds = MultiPoint(geometry_zone1).bounds\n",
    "zone1bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 2\">Task 2<a/>\n",
    "Please get your codes to work out the convex hull, envelope and coordinates of bounds for all stations point data. Then try to compare the results. Your codes below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the geometry for all stations \n",
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zoneenve = ???\n",
    "zoneenve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the total stations' bounds\n",
    "# your code here\n",
    "stabounds = ???\n",
    "stabounds "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Measure distances between points'>Measure distances between points<a/>\n",
    "\n",
    "In real society, it is often to measure the distances between two targeted observations, say point data. So we will call `Shapely` to measure the distances in different ways: e.g., Euclidean distance, Manhattan distance, Shortest-path distance and Travel-time distance. In this practical, we will use the common Euclidean distance introduced in lecture, to measure the distances between named stations within Zone 1 (**Westminster, Holborn, Bank**) as example. \n",
    "\n",
    "The primary idea is to get their coordinates from their longitude and latitude information. However, please keep in mind that the measurement unit should be **meters**, so we need to turn their `CRS` into National Grid Coordinates EPSG:27700 as the starting point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 3\">Task 3<a/>\n",
    "1. Reproject the **crs** for stations in Zone 1 data to EPSG 27700;\n",
    "2. Print out the **geometry** for targeted 3 stations. **Hint**: find the specific \"name\" from dataframe, then print out the geometry information for the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here\n",
    "dis_stations1_df=stations_zone1_df.???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you complete Task 3, we can define their coordinates to points, and get a rough idea of their relative distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# measure their distance using shapely\n",
    "from shapely.geometry import Point\n",
    "p1=Point(???) # Westminster\n",
    "p2=Point(???) # Holborn\n",
    "p3=Point(???) # Bank\n",
    "mp=MultiPoint([p1, p2, p3])\n",
    "mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To measure the distances between them, we can simply call shapely's `distance` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minimum distance from p1 to p2\n",
    "print(\"Distance between Westminster and Holborn is: \"+str(\"%.2f\"%p1.distance(p2))+\" meters\") \n",
    "# you may swap p1 and p2, see whether any difference\n",
    "# minimum distance from p2 to p3\n",
    "print(\"Distance between Holborn and Bank is: \"+str(\"%.2f\"%p2.distance(p3))+\" meters\")\n",
    "# minimum distance from p1 to p3\n",
    "print(\"Distance between Westminster and Bank is: \"+str(\"%.2f\"%p1.distance(p3))+\" meters\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It returns the minimum distance in float format, and there is also `hausdorff_distance` returning to the furthest distance. However, for our listed 3 points the results are the same, you are encouraged to replace the `distance` function with `hausdorff_distance` for comparison if interested; you are also welcome to choose other stations for calculation if any interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='Linestring Data'>Linestring Data</a>\n",
    "The second geometric type is **Polyline**, and defined as Linestring in Shapely by calling its `LineString` function. We will use the London Railways data for this section to get familiar with linestring's features. The data is available for downloading from OpenStreetMap as a shapefile data; to save time, I've already prepared it as zipfile in Github, so we should firstly read the London railways data into a geopandas DataFrame, and have a look of the features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure the download\n",
    "url  = 'https://github.com/cusp-london/Spatial-Data-Analysis/blob/master/railways.zip?raw=True'\n",
    "path = os.path.join(\"data\",\"railways.zip\")\n",
    "\n",
    "# Download\n",
    "r    = urllib.request.urlretrieve(url, path)\n",
    "\n",
    "# Unzip it into the data folder\n",
    "z    = zipfile.ZipFile(path)\n",
    "m    = z.extractall(\"data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "railways = gpd.read_file('data/railways/gis_osm_railways.shp')\n",
    "railways.plot(figsize=(10,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 4\">Task 4<a/>\n",
    "Write your codes below to answer questions on: \n",
    "1. The geometric type for railways.\n",
    "2. How many polylines in this dataset?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the outputs you got, we can tell that there are 12118 polylines, all of which are geometrically Linestring. \n",
    "Before we conduct any spatial operations between different datasets, it is good practice to check their **crs** information. Although we had already checked the **crs** for london stations point data, we need to keep other spatial layers (polyline or polygon) consistent with each other. So let's check the .crs attribute for railways data as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "railways.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may reach at a WGS84 (EPSG: 4326) crs for the railways data, and we need to set it to OSGB 1936(EPSG:27700) to match the stations data. So let's get the new geometry set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_railways=railways.to_crs(epsg=27700)\n",
    "new_railways.plot() # check the difference between 2 plots with varied CRS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will find the railway map has been stretched vertically from the view, and now it is ready to be used for distance meansurement as well. Now let's explore the spatial relationships between points and polylines:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Geo-location and spatial operation\">Geo-location and spatial operation<a/>\n",
    "    \n",
    "In this section, we want to select a specific railway line in London, preferrably a line we are familiar with; to find out the spatial relationships between this line and selected points (stations in real society); and measure the distance from point to the line, if the point is not on it. So we need to set up following goals to realize step by step:\n",
    "1. Find the **Circle lines** and save them together by selecting data calling `loc`.\n",
    "2. Locate a specific sub-line in **circle lines** by calling `iloc`.\n",
    "3. Get the geometric features for selected line, e.g. **length, bounds, buffer**.\n",
    "4. Explore the spatial relationship between point and linestring by calling `within` and `contains`.\n",
    "5. Measure the distance from a point to the linestring.\n",
    "\n",
    "Furthermore, if you go back to check the headers for railways data, you will find a `name` column with the searching keywords we are interested (They are: 'Circle line','Circle/District Line' and 'Circle and District Lines'). In order to make it easier to select data by label or by a conditional statment (`.loc`), we can reset the index for each polyline with `name`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we are going to search by name, so set name as the index\n",
    "new_railways.set_index('name', inplace=True)\n",
    "new_railways.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circleline=new_railways.loc[['Circle line','Circle/District Line', 'Circle and District Lines']]\n",
    "circleline.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "117 polylines have been selected based on our searching criterion, and combined into a new dataframe \"circleline\". For specific linestring, if we want to get more detailed information on its geometric coordinates, we can call `list` function to realize it. So in order to select a linestring object, we can use its row number by calling `iloc`; here I designedly choose the line at \"74\" (you will figure out my intention later), and we can get the points coordinates by calling `coords` as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the coordinates information on specific polyline\n",
    "list(circleline.geometry.iloc[74].coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many points have been used to define this line? You can count it from the output (if you like), or call another method `len` to get the count directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(circleline.geometry.iloc[74].coords) # measure the points building up this line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "YES, you are right, `len` and `length` are different! For polylines, `len` counts the number of points to define linestring, and for this polyline, there are 16 points. However, it is not a fixed number, and the number of points to define a linestring varies. Try any other rows, to see the variations of points counts. Further step, let's try to call `length` to see what we will get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circleline.geometry.iloc[74].length # measure the length of the railway line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we know the 16 points used to define specific section of circle line, we use their coordinates to define a new LineString named as circleline_74."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import LineString\n",
    "\n",
    "circleline_74=LineString([(530428.0074790481, 180413.6238164446),\n",
    " (530411.1982466972, 180376.18195359968),\n",
    " (530397.0032235021, 180344.56021560577),\n",
    " (530379.2209069991, 180274.64529610798),\n",
    " (530371.9918020066, 180229.70448672207),\n",
    " (530359.5922412437, 180162.43145132443),\n",
    " (530350.0843889426, 180104.51306821714),\n",
    " (530346.4988183748, 180071.82826697844),\n",
    " (530340.3595443233, 180026.33684885327),\n",
    " (530313.7082672049, 179845.00756705605),\n",
    " (530306.1320168389, 179804.34217539),\n",
    " (530297.7263252948, 179779.2342295017),\n",
    " (530282.1126631077, 179751.64929959673),\n",
    " (530261.5132050588, 179726.7408450343),\n",
    " (530225.912359408, 179690.44300068187),\n",
    " (530180.0387492746, 179656.06340498675)])\n",
    "circleline_74"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Upon the plotting of this polyline, we can call `.bounds` method to get the coordinates for its bounding boxes, which are the coordinates for the left-bottom and right-top."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circleline_74.bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Spatial relationships between point and line\">Spatial relationships between point and line<a/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you try to google the coordinates for 3 selected stations and the sub-section of circle line, you may find the Westminster station (point $p1$) should be on this line; so let's check it by calling either `within` or `contain`. Think of their relationship, it should be line contain point and point within line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check whether point is within the line\n",
    "p1.within(circleline_74)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check whether the line contains the point\n",
    "circleline_74.contains(p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, from the boolean results, both return **False**, indicating a not contain/within relationship between the point and polyline. But you know it is not the truth, so what is the problem? Discuss with your neighbors about it. \n",
    "\n",
    "Before we work out the solution, there is a roundabout to find the reason...........measure the distance between point and line!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1.distance(circleline_74)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that the station is 6.8 meters away from the line, emmmmm, reasonable. So reflecting from the reality, it could be systematic error generated when collecting the original data; or the measuring of coordinates take different centroid/center point. So a commonly used way to fix this could be use `buffer` function.\n",
    "### <a id=\"Buffer\">Buffer<a/>\n",
    "A buffer is a geometry with an envelope at a given distance from the object's envelope. So let's create a 10 metres (slightly larger than 6.8 meters) buffer zone for the named railway line. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create buffer of 10 metres on the linestring geometry\n",
    "circleline_74_buf=circleline_74.buffer(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinates for boundings for the buffered zone are also derivable by `bounds`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circleline_74_buf.bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the type of newly generated buffer data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(circleline_74_buf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is no longer Linestring, and now become Polygon. If we try to plot both the linestring and its buffer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import descartes\n",
    "ax = plt.gca()\n",
    "ax.add_patch(descartes.PolygonPatch(circleline_74_buf, fc='r', ec='b', alpha=0.4)) # buffer area in red with 40% transparency\n",
    "ax.add_patch(descartes.PolygonPatch(circleline_74.buffer(0.2))) # original linestring in black in the center\n",
    "minx, miny, maxx, maxy = circleline_74_buf.bounds\n",
    "ax.set_xlim(minx, maxx); ax.set_ylim(miny, maxy)\n",
    "ax.set_aspect('auto') # normally use 'equal', here emphasing the buffering\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Buffering helps us to widen the railway line by 10 meters, and we will further check whether the point locates within the buffered zone. So complete the following task by yourself.\n",
    "### <a id=\"Task 5\">Task 5<a/>\n",
    "1. Test whether Westminster station is in the buffered zone of selected section for circle line.\n",
    "2. Test whether the other point, Bank station, locates in the buffered zone for this line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! It is **True** now, which means Westminster station is actually on this line. We can try to check the containment of Bank station ($p3$) for this line as well, and do further measurement of the distance from $p3$ if not on this line. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3.within(circleline_74_buf)  # Bank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circleline_74_buf.contains(p3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3.distance(circleline_74_buf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you receive a **True** output, then try to interpret the result: \n",
    "It means ------------station is ------------ line. \n",
    "\n",
    "If you receive a **False** output, then further work out how far away from the line:\n",
    "------------station is------------meters away from the specified circle line section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Polygon Data\">Polygon Data<a/>\n",
    "\n",
    "Thus far we have encountered two types of geometries for shapefiles, namely **Point** and **Linestring**, and now we will recap on the most familiar third type, **Polygon**, with London LSOA data as the example. Copy your LDN-LSOAs.shp into \"data\" folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "lsoas=gpd.read_file('data/LDN-LSOAs.shp')\n",
    "lsoas.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you are familiar with Geodataframe, I leave it to you to get the information on data features by completing Task 6.\n",
    "### <a id=\"Task 6\">Task 6<a/>\n",
    "1. This is a polygon, please check the geometry for polygon and its **crs**.\n",
    "2. Try to measure the area and length of London LSOAs as polygon data. \n",
    "3. Get the coordinates of bounding box for London LSOAs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What have you got? They should be 4835 measures for lsoa's area and length respectively. \n",
    "\n",
    "Since you've get your task done, let's have the 3 datasets plotted together for visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot three layers in London together\n",
    "stations=gpd.read_file('data/stations.shp')\n",
    "new_stations=stations.to_crs(epsg=27700)\n",
    "\n",
    "ax = lsoas.plot(edgecolor='lightblue', facecolor='none', figsize=(10, 8))\n",
    "new_railways.plot(ax=ax, color='grey')\n",
    "new_stations.plot(ax=ax, color='red', alpha=0.5)\n",
    "ax.set(xlim=(500000, 560000), ylim=(160000, 200000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Spatial relationships and operations\">[Spatial relationships and operations](https://en.wikipedia.org/wiki/Spatial_relation)<a/>\n",
    "\n",
    "Geospatial data normally have their spatial relationships, e.g. overlap, intersect, contain, ..., etc., to describe how two or more spatial objects relate to each other. We had tried to explore the relationship between point and polyline about `contains` and `within`, and we will further explore their relationships with polygon.\n",
    "\n",
    "<img  height=\"100\" src=\"TopologicSpatialRelations2.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Intersections\">Intersections<a/>\n",
    "The spatial intersection relationship happens between polygons, but we are going to practice it simply with our defined three points and their buffers. As you might notice the distances among these points, we will define a 1000 meters' buffer for test. In addition, the library `decartes` will be used to facilitate visualization of these operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import MultiPolygon\n",
    "poly_1=p1.buffer(1000)\n",
    "poly_2=p2.buffer(1000)\n",
    "poly_3=p3.buffer(1000)\n",
    "mp_buf = MultiPolygon((poly_1, poly_2, poly_3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import descartes\n",
    "ax = plt.gca()\n",
    "ax.add_patch(descartes.PolygonPatch(poly_1, fc='b', ec='k', alpha=0.2))\n",
    "ax.add_patch(descartes.PolygonPatch(poly_2, fc='g', ec='k', alpha=0.2))\n",
    "ax.add_patch(descartes.PolygonPatch(poly_3, fc='r', ec='k', alpha=0.2))\n",
    "minx, miny, maxx, maxy = mp_buf.bounds\n",
    "ax.set_xlim(minx, maxx); ax.set_ylim(miny, maxy)\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Testing the relationships \n",
    "poly_1.intersects(poly_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting the object from this relationship \n",
    "int_1_2 = poly_1.intersection(poly_2)\n",
    "int_1_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The buffered polygon for $p1$ (blue) and $p2$ (green) are intersecting with each other, and. we plot out the object for intersected part. Now you may want to test the relationship between $p1$ and $p3$ by yourself.\n",
    "### <a id=\"Task 7\">Task 7<a/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Testing the relationships between poly_1 and poly_3\n",
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting the object from this relationship \n",
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(int_1_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Intersection` focuses on the shared part(s) between two polygons, but how if we want to find their differences?\n",
    "### <a id=\"Difference\">Difference<a/>\n",
    "As the 1000 metres' buffer zone for Point $p1$ and Point $p2$ intersect with each other, we will continue with this data for `Difference` operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the difference part in poly_1 from poly_2\n",
    "poly_1.difference(poly_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# difference between poly_1 and intersection part int_1_2\n",
    "poly_1.difference(int_1_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will find this is actually the buffered zone of $p1$, without its shared part with $p2$. The slight difference between these 2 outputs will be the outline of shared part. We could also do it in the other way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both `intersection` and `difference` are kind of \"subtraction\", how about \"addition\"?\n",
    "### <a id=\"Union and Cascaded Union\">Union and Cascaded Union<a/>    \n",
    "If we want to get the \"addition\" of two objects, we may call `union` in `shapely`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_1_2 = poly_1.union(poly_2)\n",
    "u_1_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Task 8'>Task 8<a/>\n",
    "Once you've understood the rationale of these 3 spatial operations, try to complete this task by yourself.\n",
    "- Union the difference part in poly_1 from poly_2, and the difference part in poly_2 from poly_1. \n",
    "- Further test the containment of point (531000, 180100) from your output of union."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You've got your outputs so far, however, if we change the method for union into `cascaded_union`, the results will vary. Try google the differences between these two union methods and discuss with your neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.ops import cascaded_union\n",
    "cu = cascaded_union([d_1_2, d_2_1, int_1_2])\n",
    "cu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have you spotted out, which element is the crucial reason for the different union outputs? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test the containment of test point again\n",
    "cu.contains(p_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you are happy with these operations on three main types of vector data. Let's define a specific region, which is City of London, and get the stations and railways within the area \"cut\" out.\n",
    "\n",
    "## <a id='Get stations and railways within City of London'>Get stations and railways within City of London<a/>\n",
    "To select only the railways within City of London, we could take several approaches. We have the lsoa layer for 4836 polygons, as well as the railway network layer for London with 12116 segments. We could then use the intersects method for each lsoa to test if it intersects with a particular segment of the railway network, and keep all the segments where we find an intersection with the lsoa.\n",
    "\n",
    "Although this would work, it turns out to be very inefficient requiring we compare each lsoas against each railway segments and test for an intersection. If we think about our problem from a slightly different perspective, we know that if we find a segment that intersects with a lsoa within **city_ldn**, it must, by definition, intersect with the polygon, if we had such a thing. This would substantially reduce the number of intersection tests we need to conduct. We now only need compare 1 polygon against each of the railway segments. Great.\n",
    "\n",
    "Before we start our computation, there is still one crucial thing to be solved: we do not yet have the magical **city_ldn** polygon, and this could be got by calling `Dissolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find city of london polygon\n",
    "lsoas.set_index('lsoa11nm', inplace=True)\n",
    "city_ldn = lsoas.loc[['City of London 001A','City of London 001B', 'City of London 001C', 'City of London 001E', 'City of London 001F', 'City of London 001G']]\n",
    "city_ldn.plot(edgecolor='lightblue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to dissolve all the lsoa boundaries that do not coincide with the boundary of the DataFrame\\'s geometry collection. This is done by creating a new attribute that takes on the same values for each feature, and calling the `dissovle` method with that attribute\n",
    "as the argument to the `by` option:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city_ldn['dummy'] = 1.0\n",
    "city_london = city_ldn.dissolve(by='dummy')\n",
    "city_london.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or else, we could also obtained this polygon by using the `unary_union` method of the GeoSeries, and this gives us a vector object combining all of the polygons in boundary layer. We would then have to construct a new GeoDataFrame with this as the Geometry column. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city_london_uu = city_ldn['geometry'].unary_union\n",
    "city_london_uu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Clip stations within city of london'>Clip stations within city of london<a/>\n",
    "<img  height=\"100\" src=\"clip.png\">\n",
    "When you clip a vector data set with another layer, you remove points, lines or polygons that are outside of the spatial extent of the area that you use to clip the data. (Image Source: ESRI)\n",
    "\n",
    "We can do the `clip` operation by `intersects`: We firstly use the above created unified polygon **city_london_uu** to represent City of London area covered by clip layer. Next we can use the `intersects` method to select just the points within the stations object that fall within the geometry in **city_london_uu** object. It will return a boolean mask for every point and further subset the stations object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stations_clip = new_stations[new_stations.geometry.intersects(city_london_uu)]\n",
    "stations_clip.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Clip railways within City of London'>Clip railways within City of London<a/>\n",
    "    \n",
    "In this section, we started with building a spatial index to realize our clipping goal. If you are interested to explore further on the in-depth rationale, it is really worthwhile to go through [Boeing's work](https://geoffboeing.com/2016/10/r-tree-spatial-index-python/), and [further reading](http://www.dpi.inpe.br/livros/bdados/artigos/oracle_r_tree.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_index = new_railways.sindex\n",
    "# Create a box for the initial intersection\n",
    "bbox = city_london_uu.bounds\n",
    "# Get a list of railway polylines that overlap the bounding box and subset the data to just those lines\n",
    "sidx = list(spatial_index.intersection(bbox))\n",
    "rail_sub = new_railways.iloc[sidx]\n",
    "\n",
    "# Clip the data - with these data\n",
    "railways_clip = rail_sub.copy()\n",
    "railways_clip['geometry'] = rail_sub.intersection(city_london_uu)\n",
    "\n",
    "# Return the clipped layer with no null geometry values\n",
    "railways_clip[railways_clip.geometry.notnull()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 227 segments of railways fall in City of London as shown. So let's plot the stations and railways within City of London together for a view.\n",
    "### <a id=\"Task 9\">Task 9<a/>\n",
    "Recall you plotting skillks acquired before, plot three layers on stations, railways and LSOAs for City of London together. \n",
    "\n",
    "**Hint**: assign color 'grey' for railways, 'yellow' for stations, and edgecolor 'red', facecolor 'lightblue' for lsoas. Legend should be present, and be careful about your axis setting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we go, let's check the containment of defined stations $p2$ and $p3$ in the city of london boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2.within(city_london_uu) # Holborn station"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city_london_uu.contains(p3) # Bank station"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you may find Holborn station is not in City of London, but how far is it away from the area? and one step further, will a 1000 meters' buffer zone of Holborn station happened to overlap City of London somehow? Let's try practice our techniques acquired before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "p2.distance(city_london_uu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "poly_1000=p2.buffer(1000)\n",
    "poly_1000.intersects(city_london_uu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ! (Optional) Other realization of clipping :\n",
    "If time permits, you may want to try other two ways to clip the railways in City of London: \n",
    "\n",
    "(1) use **lambdas**\n",
    "\n",
    "(2) use loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = new_railways['geometry']\n",
    "rc_routes = r[r.apply(lambda x: x.intersects(city_london.iloc[0]['geometry']))]\n",
    "rc_routes.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_routes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.gca()\n",
    "rc_routes.plot(ax=ax, edgecolor='grey')\n",
    "city_london.plot(ax=ax, color='lightblue')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do you think? **Lambdas** are handy, but tend to be a bit more difficult to read as \\\"anonymous functions\\\". So let's try another more \\\"transparent\\\" approach by using a simple loop. The rationale is to test the intersection between each railway segment and City of London, if the result returns to be true, then append the segment to a list and store the list of final segments intersecting with City of London."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rails = [] # create an empty list for further storage\n",
    "for ids, route in enumerate(rc_routes):\n",
    "    print(ids)\n",
    "    rails.append(route.intersection(city_london.iloc[0]['geometry']))\n",
    "# use list of intersection segments above into a GeoSeries:\n",
    "rc_l = gpd.GeoSeries(rails)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 8)\n",
    "ax = plt.gca()\n",
    "city_london.plot(ax=ax, color='lightblue')\n",
    "rc_l.plot(ax=ax, edgecolor='grey')\n",
    "ax.set_xlim(530900, 533900); ax.set_ylim(180500, 182200)\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(rc_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_l = gpd.GeoDataFrame({'geometry': rc_l})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_l.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city_ldn.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is more tidy than **lambada** method, and you may be interested to check the spatial relationships between rc_l and city_london then, feel free to write your own codes. \n",
    "\n",
    "**Hint**: `contains`, `crosses` and `touches`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='Join the spatial data'>Join the spatial data<a/>\n",
    "\n",
    "We had introduced you how to join or merge dataframes through functionality provided in Pandas, see [full documentation](https://pandas.pydata.org/pandas-docs/stable/merging.html)；we have seen the notion of spatial relationships between geometry objects: within, contains, intersects, ..., now we are going to transfer attributes from one layer to another based on their spatial relationship by using `geopandas.sjoin` function with following elements:\n",
    "- The GeoDataFrame to which we want add information\n",
    "- The GeoDataFrame that contains the information we want to add\n",
    "- The spatial relationship we want to use to match both datasets ('intersects', 'contains', 'within')\n",
    "- The type of join: left or inner join\n",
    "\n",
    "Since the stations and railway network clipped to the extent of City of London have prepared already, let's determine which LSOAs intersect the railway network within City of London. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spatial join, lsoas with stations\n",
    "lsoas_with_stations = gpd.sjoin(city_ldn, stations_clip, op='contains', how='inner')\n",
    "lsoas_with_stations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we scroll to the right of the DataFrame output, we see a column labeled \"index_right\". The values in this column indicate the index of the features in the right DataFrame (in our case the stations) within the feature in the current row of the left DataFrame (the lsoas).\n",
    "\n",
    "You may noticed that spatial join operation doesn't change the geometries. In another word, we are not joining respective geometries, but joining attributes based on a spatial relationship between the geometries. This means that the geometries need to at least overlap partially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsoas_with_stations.shape # check the number of rows and columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we need to further clean our joined data by dropping some useless columns, and further check the **crs** for all relevant layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=lsoas_with_stations.drop(['dummy', 'index_right', 'display_na'], axis=1)\n",
    "a.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city_ldn.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_l.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are not wrong! It should be blank, as the railway DataFrame does not have a **crs**. To correct this, we need to set it to that of the city_ldn dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc_l.crs = city_ldn.crs # create a crs for the rc_l\n",
    "rc_l = rc_l.to_crs(city_ldn.crs) # update the coordinates accordingly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All prepared, and let's join!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spatial join again\n",
    "lsoas_with_stations_railways = gpd.sjoin(a, rc_l, how='inner', op='intersects')\n",
    "lsoas_with_stations_railways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shape for joined data with lsoa, railways and stations in City of London\n",
    "lsoas_with_stations_railways.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, three layers had joined together, and you may want to save it into shapefile data for further use. Please do as you wish if interested."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Buffering\">Buffering<a/>\n",
    "\n",
    "Buffers are one of the classical GIS operations, in that with a particular geometric object as the target, we define a specific radious/distance from its envelope, then a geometry is drawn around the target with an envelope at the defined distance. These are very useful, for instance, in combination with point-in-polygon operations to calculate accessibility, catchment areas, etc. In last practical, we had actually started to use the point-in-polygon operation when getting the heatmap for London pubs, through counting the number of pubs in each LSOA. Buffer could be taken one step further, on generating guard areas that can be added for a specified distance (or a range) around points, lines or polygons. Besides, buffers could play the role as basis for specific spatial operations, e.g. overlay, through calling `buffer` in `Shapely`, and we will further use it for next week. \n",
    "\n",
    "We can use `geopandas` by simply calling the `buffer` method, passing it the the radius. Mind that the radius needs to be specified in the same units as the CRS of the geometry you are working with, hence we with generate 50 meters and 100 meters buffers respectively for each dataset to compare the effects from radius visually. Now let's try to produce buffers for **varied shapes** in City of London."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Buffers could be applied not only on **points** data, but also **polylines** and **polygons** data. So in this section, we will use the stations, railways and lsoas within the boundary of City of London as an example, to show you the buffering function for specified points, polylines and polygons. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us firstly get a 50 metres buffer around each station in City of London, and then 100 metres buffer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buf_50_s = stations_clip.buffer(50)\n",
    "buf_100_s = stations_clip.buffer(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 8)\n",
    "ax = plt.gca()\n",
    "buf_100_s.plot(ax=ax,color='lightblue', edgecolor='grey', alpha=0.2)\n",
    "buf_50_s.plot(ax=ax,color='red', edgecolor='grey', alpha=0.2)\n",
    "stations_clip.plot(ax=ax, color='yellow', edgecolor='k')\n",
    "ax.set_xlim(530900, 533900); ax.set_ylim(180500, 182200)\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The buffering operation for polyline, railways in city of london, is similar to what we've done to points data, so let's try to complete task 6.\n",
    "### <a id=\"Task 10\">Task 10<a/>\n",
    "1. Create 50 meters and 100 meters buffers for railways in City of London respectively;\n",
    "2. Plot them below following the similar methods as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here \n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buf_50_c=city_ldn.buffer(50)\n",
    "buf_100_c=city_ldn.buffer(100)\n",
    "\n",
    "ax = plt.gca()\n",
    "buf_100_c.plot(ax=ax,color='lightblue', edgecolor='grey', alpha=0.2)\n",
    "buf_50_c.plot(ax=ax,color='red', edgecolor='grey', alpha=0.2)\n",
    "city_ldn.plot(ax=ax, color='white', edgecolor='k')\n",
    "ax.set_xlim(530800, 534000); ax.set_ylim(180300, 182500)\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discuss with your neighbors: How if we choose a negative buffer distance? What will happen to the geometry? How about a zero buffer then? Have your codes here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buf_neg50_c=city_ldn.buffer(-50)\n",
    "buf_neg50_c.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buf_zero_c=city_ldn.buffer(0)\n",
    "buf_zero_c.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to read the ballpark boundary of buffers, for example the 100 meters' buffers for varied shapes, you may call `Bounding Box` to make the smallest possible rectangle that completely encloses your buffered area. The output could be just a tuple of the \"minx\", \"miny\", \"maxx\", \"maxy\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "buf_s_bb=buf_100_s.bounds\n",
    "buf_r_bb=buf_100_r.bounds\n",
    "buf_c_bb=buf_100_c.bounds\n",
    "print(buf_s_bb, buf_r_bb, buf_c_bb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Credits!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Geocomp",
   "language": "python",
   "name": "gsa2018"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}