5SSG2060_Practical_4_PPA_PointDistribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5SSG2060_Week4_Point Pattern Analysis (PPA)  I\n",
    "<a href=\"#This Week's Overview\">This Week's Overview</a>\n",
    "\n",
    "<a href=\"#Learning Outcomes\">Learning Outcomes</a>\n",
    "\n",
    "<a href=\"#Get prepared\">Get prepared</a>\n",
    "\n",
    "<a href=\"#What is Point Pattern Analysis (PPA)\">What is Point Pattern Analysis (PPA)</a>\n",
    "\n",
    "<a href=\"#Data processing\">Data processing</a>\n",
    "- <a href=\"#Generate random point data\">Generate random point data</a>\n",
    "- <a href=\"#Generating CSR\">Generating CSR</a>\n",
    "- <a href=\"#Convex hull\">Convex hull</a>\n",
    "\n",
    "<a href=\"#Centrography analysis\">Centrography analysis</a>\n",
    "- <a href=\"#Mean center $(x_{mc},y_{mc})$\">Mean center $(x_{mc},y_{mc})$</a> \n",
    "- <a href=\"#Weighted mean center $(x_{wmc},y_{wmc})$\">Weighted mean center $(x_{wmc},y_{wmc})$</a> \n",
    "- <a href=\"#Manhattan Median $(x_{mm},y_{mm})$\">Manhattan Median $(x_{mm},y_{mm})$</a>\n",
    "- <a href=\"#Euclidean Median $(x_{em},y_{em})$\">Euclidean Median $(x_{em},y_{em})$</a> \n",
    "\n",
    "<a href=\"#Dispersion and Orientation\">Dispersion and Orientation</a>\n",
    "- <a href=\"#Standard Distance & Standard Distance Circle\">Standard Distance & Standard Distance Circle</a>\n",
    "- <a href=\"#Standard Deviational Ellipse\">Standard Deviational Ellipse</a>\n",
    "- <a href=\"#Shape Analysis\">Shape Analysis</a>\n",
    "\n",
    "<a href=\"#Nearest Neighbors and Statistics\">Nearest Neighbors and Statistics</a>\n",
    "- <a href=\"#Nearest Neighbors\">Nearest Neighbors</a>\n",
    "- <a href=\"#Mean Nearest Neighbor Distance Statistics\">Mean Nearest Neighbor Distance Statistics</a>\n",
    "- <a href=\"#Nearest Neighbor Distance Functions\">Nearest Neighbor Distance Functions</a>\n",
    "- <a href=\"#Interevent Distance Functions\">Interevent Distance Functions</a>\n",
    "\n",
    "<a href=\"#K Nearest Neighbors (KNN)\">K Nearest Neighbors (KNN)</a>\n",
    "- <a href=\"#Weights\">Weights</a>\n",
    "  - <a href=\"#k-nearest neighbor weights\">k-nearest neighbor weights</a>\n",
    "  - <a href=\"#Distance band weights\">Distance band weights</a>\n",
    "  - <a href=\"#Kernel weights\">Kernel weights</a>\n",
    "\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 the following two weeks, we will develop our capability in using `pointpats` from `PySAL` to realize spatial point pattern analysis, visualization and exploration. Tasks for this week will mainly focus on visualizing spatial point patterns and exploring *distance-based* spatial point patterns analysis; while on Week 5 we will spend most of the time exploring *density-based* spatial point patterns analysis methods.\n",
    "\n",
    "In this practical, besides of the first section on CSR process using simulated points, all rest tasks will rely on London pubs data to develop your understanding on:\n",
    "\n",
    "- Point processing.\n",
    "- Centrography and visualization.\n",
    "- \n",
    "- Distance based methods.\n",
    "\n",
    "\n",
    "## <a id=\"Learning Outcomes\">Learning Outcomes</a>\n",
    "\n",
    "By the end of this practical you should be able to:\n",
    "- explore how random numbers and random points can be generated\n",
    "- appreciate the nature of uniform random points and how they are distributed \n",
    "- measure the centrality and the spatial dispersion of random points\n",
    "- familiarise ourselves with distance based point pattern analysis methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Get prepared\">Get prepared</a>\n",
    "We will import the general libraries at the beginning of this practical, and leave specific ones for corresponding sections for clarification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import urllib\n",
    "import zipfile\n",
    "import re\n",
    "\n",
    "import libpysal as lps\n",
    "from shapely import geometry as sgeom\n",
    "from shapely.geometry import Point\n",
    "import descartes as des\n",
    "import pysal as ps\n",
    "import pointpats \n",
    "from pointpats import PointPattern\n",
    "from pysal.contrib import shapely_ext\n",
    "from pointpats import PoissonPointProcess as ppp_csr\n",
    "from pointpats import as_window\n",
    "from pysal.contrib.viz import mapping as maps # For maps.plot_choropleth\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.collections as mplc\n",
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"What is Point Pattern Analysis (PPA)\">What is Point Pattern Analysis (PPA)</a>\n",
    "\n",
    "What are points in space? Should we simply take them as locations of interest? We could normally categorize points into 2 groups as **Event Points** and **Arbitrary Points**. \n",
    "- *Event Points* are locations where something of interest has occurred, hence could represent [a wide variety of phenomena](https://en.wikipedia.org/wiki/Point_process). \n",
    "- *Arbitrary Points* are locations where the phenomena of interest has not been observed, such as the so-called \"empty space\" or \"regular\" points. \n",
    "\n",
    "As our practicals mostly looking at real life data in our society, so you will be given **Event Points** to observe their collective spatial pattern characteristics. Let's recall the concept of `complete spatial randomness (CSR)` introduced in lecture: if there are a series of *point locations*, $(p_1, p_2, \\ldots, p_n)$ in a two-dimensional study region $\\Re$, then for point location $i$, its spatial identity should be $p_i = (x_i, y_i)$. \n",
    "\n",
    "If it is a `CSR` process, then the points are independent from one another (generated from stochastic process, and with constant probability); \n",
    "\n",
    "If the point patterns are not under `CSR` (which is the normal case), they could be affected by **First order effects** and **Second order effects** covered in lecture. \n",
    "- *First order effects*: Non-CSR point process' intensity pattern driven by underlying covariate.\n",
    "- *Second order effects*: Non-CSR point process due to interaction and dependence between events in space. \n",
    "\n",
    "So to further investigate the process' spatial structure, detect the points' deviations from `CSR` and test the point pattern statistics, we will use the main mathods from point pattern analysis module `pointpats` in `PySAL` as below:\n",
    "1. Point Processing\n",
    "2. Centrography and Visualization\n",
    "3. Nearest Neighbors\n",
    "4. Other Distance Based Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Data processing\">Data processing</a>\n",
    "\n",
    "Set up a **\"data\"** folder in your directory for this notebook, and copy the data we've used for last week (LDN-LSOAs shapefile) into your data folder. \n",
    "\n",
    "### <a id=\"Generate random point data\">Generate random point data</a>\n",
    "In order to simulate the `CSR` process, let's start with generating some random series of data by calling `random` in `numpy`. It will get you 100 random data between 0 and 1, with 6 decimals. Upon the execution of the code, you will find a .txt file named as \"random_points\" generated in your local directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test randomness\n",
    "import random\n",
    "for i in range(100):\n",
    "    print ('%.6f' % random.random(), '%.6f' % random.random())    \n",
    "file = open('random_points.txt', 'w')\n",
    "\n",
    "for i in range(100):\n",
    "    file.write('%.6f' % random.random())\n",
    "    file.write('%.6f' % random.random())\n",
    "    file.write(\"\\n\")\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the random data to generate random points \n",
    "df = pd.DataFrame(np.random.rand(50, 2), columns=['a', 'b'])\n",
    "\n",
    "# plot out the random points\n",
    "df.plot(kind='scatter', x='a', y='b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id=\"Task 1\">Task 1</a>  \n",
    "\n",
    "1. Let us think of the mean centre of the random points. Can you calculate the spatial mean? If the points are uniformly randomly distributed, between (0,0) and (1,1), where would the spatial mean most likely be found? Why?\n",
    "\n",
    "2. Repeat producing another set of random points, calculate the spatial mean, and see how they change their position from one set to another. Is there a way to make the spatial mean more stable?\n",
    "\n",
    "3. Just as we do in ordinary statistics, we measure the scatteredness of observations on a single variable about their mean by the standard deviation. In other words, we can measure the degree of spatial dispersion of a point pattern by its standard distance using square root of the average square of the distances from every point to the spatial mean. Try calculating this and produce the circle of influence around the spatial mean.\n",
    "\n",
    "** If you feel any stuck in this task, don't be \"panic\" at all, just keep continuing your practical and you will get the hints/answers on your halfway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your answers here or turn back later"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Generating CSR\">Generating CSR</a>\n",
    "\n",
    "Now let's move one step back into our London reality, to try to simulate a 1000-point dataset within London boundary from a `CSR (complete spatial randomness)` process. From your lecture, you may realize that `CSR` events follow a homogeneous Poisson Process over the study region $\\Re$. Hence, point pattern is considered to be number of events occurring in sub-regions A, of the whole study region $\\Re$.\n",
    "\n",
    "<img src=\"Picture1.png\" style=\"width: 100px;\"/>\n",
    "\n",
    "where Y(A) is the number of events occurring in the area A, following a Poisson distribution with mean λA.\n",
    "<img src=\"Picture2.png\" style=\"width: 200px;\"/>\n",
    "\n",
    "So given N events in A, they are (1) independent random sample from a uniform distribution with equal probability of occurring at any position, indicating no first order effects; (2) independent of the position of any other, implying no spatial interaction with one another.\n",
    "\n",
    "In the following section, we will call `PoissonPointProcess` in `pointpats` to generate the dataset. It is similar with the previous practice on getting random point data: we generate events with x coordinates from a uniform distribution on $(x_1,x_2,...,x_i)$ and y coordinates from a uniform distribution on $(y_1,y_2,...,y_i)$. So a real-world observed pattern of points can be compared with the simulated ones based on CSR, in order to assess whether observed patterns are regular, clustered or random distribution.\n",
    "\n",
    "For more references, you may find the `pointpats` documents at <https://pointpats.readthedocs.io/en/latest/index.html>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open London lsoa polygon shapefile\n",
    "lsoas = ps.open ('data/LDN-LSOAs.shp') \n",
    "# define the polygon shapes from London lsoa data\n",
    "polys = [shp for shp in lsoas] \n",
    "# Create the exterior polygons for Greater London from the union of the polygon shapes\n",
    "boundary = shapely_ext.cascaded_union(polys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# generate 1000 points following csr process\n",
    "# define its region boundary as london boundary\n",
    "pp = ppp_csr(as_window(boundary), 1000, 1, asPP=True).realizations[0] \n",
    "# You may find realizations at index 0, which means first realized point pattern\n",
    "# plot your point pattern out\n",
    "pp.plot(window=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 2\">Task 2</a>\n",
    "Download the London Pubs shapefile data from my github into your data folder, and get it ready for use. \n",
    "Hint: open point data from shapefile by calling `PointPattern` function, you need to pay attention to its array nature. This data will be used for following tasks, with coordinates for UTM zone 30U."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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_Pubs.zip?raw=true'\n",
    "path = os.path.join(\"data\",\"London_Pubs.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": [
    "f = ps.open('data/London_Pubs/London_Pubs.shp')\n",
    "# get the points\n",
    "pp_pubs = PointPattern(np.asarray([pnt for pnt in f]))\n",
    "f.close() # why we need to call f.close()?\n",
    "# get the points plotted\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following, we will use the 1000 CSR points as example to illustrate various spatial point pattern characteristics, and let you to complete corresponding tasks by using London pubs data.\n",
    "### <a id=\"Convex Hull\">Convex Hull</a>\n",
    "The 1000 points following CSR process were now plotted within the London boundary, now let's try something different. Have you spotted any differences and why? What is `convex hull`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp.plot(window=True, hull=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Convex hull`  is the set of all convex combinations of the points. Does it look like an envelope? That's why it is also called `convex envelope` or `convex closure` in an affine space over the points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 3\">Task 3</a>\n",
    "Get the convex hull for London pubs data below, and compare the output with plot above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Centrography analysis\">Centrography analysis</a>\n",
    "\n",
    "Once we were given a point data set, normally we want to first detect where the center is. Hence the center point of two-dimensional point data distribution pattern could be termed as `Central Tendency`, and could be measured by calling these four listed functions in order. \n",
    "1. `mean_center`:  the mean center of the unmarked point pattern.\n",
    "2. `weighted_mean_center`:  the weighted mean center of the marked point pattern.\n",
    "3. `manhattan_median`:  the manhattan median\n",
    "4. `euclidean_median`:  the Euclidean median\n",
    "\n",
    "As they have their respective pros and cons, so for your own project data in the future, the appropriate measures should be selected carefully according to your specific objective and corresponding data description. But now let's just have a general view of all of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import four centrography analysis functions \n",
    "from pointpats.centrography import mean_center, weighted_mean_center, manhattan_median, euclidean_median"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id=\"Mean center $(x_{mc},y_{mc})$\">Mean center $(x_{mc},y_{mc})$</a> \n",
    "Calculate the mean value for $x$ and $y$ respectively.\n",
    "\n",
    "$$x_{mc}=\\frac{1}{n} \\sum^n_{i=1}x_i$$\n",
    "$$y_{mc}=\\frac{1}{n} \\sum^n_{i=1}y_i$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the mean centre for a series of point events, based on their (x,y) coordinates.\n",
    "mc = mean_center(pp.points)\n",
    "#plot point pattern \"pp\", convex hull, and Minimum Bounding Rectangle\n",
    "pp.plot(title='Centers',  hull=True , window=True )\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center') # r indicates the color for marker, ^ represents icon shape\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to get the $x,y$ information for Mean Centre?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Weighted mean center $(x_{wmc},y_{wmc})$\">Weighted mean center $(x_{wmc},y_{wmc})$</a> \n",
    "\n",
    "Weighted mean center is for marked point patterns, with specification for the weight for each event point.\n",
    "\n",
    "$$x_{wmc}=\\sum^n_{i=1} \\frac{w_i x_i}{\\sum^n_{i=1}w_i}$$\n",
    "$$y_{wmc}=\\sum^n_{i=1} \\frac{w_i y_i}{\\sum^n_{i=1}w_i}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the weights\n",
    "weights = np.arange(1000)\n",
    "wmc = weighted_mean_center(pp.points, weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you try to get the $x,y$ information for Weighted Mean Centre?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp.plot(window=True) #use method \"plot\" to visualize point pattern\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center') \n",
    "plt.plot(wmc[0], wmc[1], 'gd', label='Weighted Mean Center')\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Manhattan Median $(x_{mm},y_{mm})$\">Manhattan Median $(x_{mm},y_{mm})$</a> \n",
    "\n",
    "We are familiar with Median in one-dimensional space as the middle number separating the dataset by half. However in two-dimensional space, Manhattan median is defined as a tuple whose first element is the median of $x$ coordinates, and second element is the median of $y$ coordinates, which finds the location minimized the absolute distance to all the event points.  \n",
    "\n",
    "If you can recall your exercise on `Numpy` last semester, `PySAL` handle the Manhattan median the same way as numpy.median: return the average of the two middle values.\n",
    "\n",
    "$$min  f(x_{mm},y_{mm})= \\sum^n_{i=1}(|x_i-x_{mm}|+|y_i-y_{mm}|)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For your understanding of the median_centre calculation, we may try to interpret its function by calling `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def median_center(points, crit=0.0001):\n",
    "    points = np.asarray(points)\n",
    "    x0, y0 = points.mean(axis=0)\n",
    "    dx = np.inf\n",
    "    dy = np.inf\n",
    "    iteration = 0\n",
    "    while np.abs(dx) > crit or np.abs(dy) > crit:\n",
    "        xd = points[:, 0] - x0\n",
    "        yd = points[:, 1] - y0\n",
    "        d = np.sqrt(xd*xd + yd*yd)\n",
    "        w = 1./d\n",
    "        w = w / w.sum()\n",
    "        x1 = w * points[:, 0]\n",
    "        x1 = x1.sum()\n",
    "        y1 = w * points[:, 1]\n",
    "        y1 = y1.sum()\n",
    "        dx = x1 - x0\n",
    "        dy = y1 - y0\n",
    "        iteration +=1 \n",
    "        print(x0, x1, dx, dy, d.sum(), iteration)\n",
    "        x0 = x1\n",
    "        y0 = y1\n",
    "               \n",
    "    return x1, y1\n",
    "\n",
    "median_center(pp.points, crit=.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may already have read the $x,y$ information for Manhattan median centre for this output. Is it largely different from the mean center and weighted mean center coordinates?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Manhattan Median calculation\n",
    "mm = manhattan_median(pp.points)\n",
    "# plot\n",
    "pp.plot(window=True)\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center')\n",
    "plt.plot(wmc[0], wmc[1], 'gd', label='Weighted Mean Center')\n",
    "plt.plot(mm[0], mm[1], 'yv', label='Manhattan Median') # \"y\" means yellow, \"v\" defines the icon's shape\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Euclidean Median $(x_{em},y_{em})$\">Euclidean Median $(x_{em},y_{em})$</a> \n",
    "\n",
    "Euclidean Median is the location with minimum sum of the Euclidean distances to all points. It has been widely applied in location allocation/selection problems in the real world, and is an optimization solution in nature. In this practical, we will use first iterative algorithm (Kuhn and Kuenne, 1962) to approximate Euclidean Median.\n",
    "\n",
    "$$min  f(x_{em},y_{em})= \\sum^n_{i=1} \\sqrt{(x_i-x_{em})^2+(y_i-y_{em})^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Euclidean Median calculation\n",
    "em = euclidean_median(pp.points)\n",
    "pp.plot(window=True, hull=True)\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center')\n",
    "plt.plot(wmc[0], wmc[1], 'gd', label='Weighted Mean Center')\n",
    "plt.plot(mm[0], mm[1], 'yv', label='Manhattan Median')\n",
    "plt.plot(em[0], em[1], 'm+', label='Euclidean Median')\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 4\">Task 4</a>\n",
    "Use the London Pubs data to calculate its **mean center**, **weighted mean center**, **Manhattann median**, **Euclidean median** respectively, get the results plotted out in one figure; list their $(x,y)$ coordinates information as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here for plot\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here for point coordinates\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Dispersion and Orientation\">Dispersion and Orientation</a>\n",
    "\n",
    "### <a id=\"Standard Distance & Standard Distance Circle\">Standard Distance & Standard Distance Circle</a>\n",
    "\n",
    "Standard distance is closely related to standard deviation of a dataset, and measures the dispersion of the events from their mean center $(x_m,y_m)$; hence a summary circle (standard distance circle) for the point pattern could be plotted based on the measurements, centered at $(x_m,y_m)$ with radius $SD$.\n",
    "\n",
    "$$SD = \\displaystyle \\sqrt{\\frac{\\sum^n_{i=1}(x_i-x_{m})^2}{n} + \\frac{\\sum^n_{i=1}(y_i-y_{m})^2}{n}}$$\n",
    "\n",
    "We will call the `std_distance` function to calculate the standard distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pointpats.centrography import std_distance,ellipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stdd = std_distance(pp.points)\n",
    "stdd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circle1=plt.Circle((mc[0], mc[1]),stdd,color='g')\n",
    "ax = pp.plot(get_ax=True, title='Standard Distance Circle', window=True)\n",
    "ax.add_artist(circle1)\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center')\n",
    "ax.set_aspect('equal')\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 5\">Task 5</a>\n",
    "Calculate the standard distance for London pubs and get the result plotted in yellow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id=\"Standard Deviational Ellipse\">Standard Deviational Ellipse</a>\n",
    "\n",
    "The aforementioned `standard distance circle` measures dispersion using a sole parameter $SD$; whilst standard deviational ellipse measures dispersion and trend in two dimensions: angle of rotation $\\theta$, dispersion along major axis $s_x$ and dispersion along minor axis $s_y$:\n",
    "\n",
    "* Major axis defines the direction of maximum spread in the distribution. $s_x$ is the semi-major axis (half the length of the major axis):\n",
    "\n",
    "$$ s_x = \\displaystyle \\sqrt{\\frac{2(\\sum_{i=1}^n (x_i-\\bar{x})\\cos(\\theta) - \\sum_{i=1}^n (y_i-\\bar{y})\\sin(\\theta))^2}{n-2}}$$\n",
    "\n",
    "* Minor axis defines the direction of minimum spread and is orthogonal to major axis. $s_y$ is the semi-minor axis (half the length of the minor axis):\n",
    "\n",
    "$$ s_y = \\displaystyle \\sqrt{\\frac{2(\\sum_{i=1}^n (x_i-\\bar{x})\\sin(\\theta) - \\sum_{i=1}^n (y_i-\\bar{y})\\cos(\\theta))^2}{n-2}}$$\n",
    "\n",
    "* The ellipse is rotated clockwise through an angle $\\theta$:\n",
    "\n",
    "$$\\theta = \\displaystyle \\arctan{\\{ (\\sum_i(x_i-\\bar{x})^2-\\sum_i(y_i-\\bar{y})^2) + \\frac{[(\\sum_i(x_i-\\bar{x})^2-\\sum_i(y_i-\\bar{y})^2)^2 + 4(\\sum_i(x-\\bar{x})(y_i-\\bar{y}))^2]^\\frac{1}{2}}{2\\sum_i(x-\\bar{x})(y_i-\\bar{y})}\\}}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sx, sy, theta = ellipse(pp.points)\n",
    "sx, sy, theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_degree = np.degrees(theta) #need degree of rotation to plot the ellipse\n",
    "theta_degree # The Standard Deviational Ellipse for the point pattern is rotated clockwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Ellipse\n",
    "from pylab import figure, show,rand\n",
    "fig = figure()\n",
    "e = Ellipse(xy=mean_center(pp.points), width=sx*2, height=sy*2, angle=-theta_degree) #angle is rotation in degrees (anti-clockwise)\n",
    "ax = pp.plot(get_ax=True, title='Standard Distance Circle', window=True)\n",
    "ax.add_artist(e)\n",
    "e.set_clip_box(ax.bbox)\n",
    "e.set_facecolor([0.4,0,0])\n",
    "e.set_edgecolor([0,0,0])\n",
    "ax.set_xlim(500000,560000)\n",
    "ax.set_ylim(155000,210000)\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center')\n",
    "plt.legend(numpoints=1)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 6\">Task 6</a>\n",
    "Present your work on exploring the standard deviation of london pubs dataset. Change the color for mean center into yellow. Hint: you need to consider to change the xlim and ylim."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Shape Analysis\">Shape Analysis</a>\n",
    "\n",
    "At the beginning of this practical, we've tried to define `hull` by caluculating the convex hull of point pattern. The full function is actually a measure for *shape analysis* of point pattern. So we are going to explore two methods: `Convex Hull (Hull)` and `Minimum Bounding Rectange (mbr)`, as follows.\n",
    "\n",
    "1. [Convex Hull](https://en.wikipedia.org/wiki/Convex_hull) is the smallest convex set that contains a point pattern *pp*, and could get realized through calling function **hull**.\n",
    "\n",
    "2. [Minimum Bounding Rectangle (Box)](https://en.wikipedia.org/wiki/Minimum_bounding_rectangle) is the same as the minimum bounding Rectangle of its convex hull, which is bigger than convex hull, by calling **mbr** function to calculate the vertices' values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import relevant methods\n",
    "from pointpats.centrography import hull, mbr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the points defining convex hull\n",
    "hull(pp.points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the points defining minimum bounding rectangle\n",
    "mbr(pp.points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Could you write the four vertices of the minimum bounding rectangle then?\n",
    "\n",
    "(???, ???), (???, ???), (???, ???), (???, ???)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot both Convex hull and Minimum Bounding Rectangele once\n",
    "# your code here\n",
    "pp.plot(title='Centers',  hull=True , window=True )\n",
    "plt.plot(mc[0], mc[1], 'r^', label='Mean Center')\n",
    "plt.plot(wmc[0], wmc[1], 'yd', label='Weighted Mean Center')\n",
    "plt.plot(mm[0], mm[1], 'g*', label='Manhattan Median')\n",
    "plt.plot(em[0], em[1], 'm+', label='Euclidean Median')\n",
    "plt.legend(numpoints=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Rectangle\n",
    "ax=pp.plt(get_ax=True, title ='Convex hull and MBR', hull=True)\n",
    "mbr_pts=mbr(pp.points)\n",
    "mbr_patch=Rectangle((mbr_pts[0], mbr_ptc[1], mbr_pts[2]-mbr_pts[0], \n",
    "                     mbr_pts[3]-mbr_pts[1], alpha=.2))\n",
    "ax.add_patch(mbr_patch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have you found it quite identical to what  we've got for \"Centrography analysis\"? Think about the reason and discuss with your neighbours. Now let's play with different combination of the \"blocks\".\n",
    "### <a id=\"Task 7\">Task 7</a>\n",
    "Please get your code below to plot the london pubs' mean center and the standard distance circle, with convex hull presented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"Nearest Neighbors and Statistics\">Nearest Neighbors and Statistics</a>\n",
    "\n",
    "Neighbors-based methods are known as non-generalizing machine learning methods, and there are many learning routines relying on nearest neighbors at their core. The working mechanism is that, they try to “remember” all of its training data (possibly transformed into a fast indexing structure such as a `Ball Tree` or `KD Tree`), build up the distance metrics through searching the distances between targeted point and its neighbors, and get the results of nearest neighbor point to each targeted point, as well as the corresponding distance. It can be realized through calling `NearestNeighbors` in `sklearn.neighbors`, with input data as either `NumPy` arrays or scipy.sparse matrices (the latter doesn't work for our data), and outputs as indices for the nearest neighbor to each point, together with respective distance. \n",
    "\n",
    "We will use the London pubs data to get nearest neighbors and the relevant distance based statistical methods.\n",
    "* [Nearest Neighbors](#Nearest-Neighbors)\n",
    "* [Mean Nearest Neighbor Distance Statistics](#Mean-Nearest-Neighbor-Distance-Statistics)\n",
    "* [Nearest Neighbor Distance Functions](#Nearest-Neighbor-Distance-Functions)\n",
    "* [Interevent Distance Functions](#Interevent-Distance-Functions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Nearest Neighbors\">Nearest Neighbors</a>\n",
    "One can use the `KDTree` or `BallTree` classes directly to find nearest neighbors. This is the functionality wrapped by the NearestNeighbors class. The Ball Tree and KD Tree have the same interface; we’ll show an example of using the `BallTree` here, and let you do the task by using `KDTree`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import NearestNeighbors\n",
    "import fiona\n",
    "from shapely.geometry import shape\n",
    "\n",
    "pp_pubs_2 = fiona.open('data/London_pubs/London_pubs.shp') # read the london pubs point shapefile data\n",
    "geoms = [ shape(feat['geometry']) for feat in pp_pubs_2 ] # read the geometry\n",
    "pp_pubs_arrays = [ np.array((geom.xy[0][0], geom.xy[1][0])) for geom in geoms ] # turn points into point array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbrs = NearestNeighbors(n_neighbors=2, algorithm='ball_tree').fit(pp_pubs_arrays) # use ball_tree\n",
    "distances, indices = nbrs.kneighbors(pp_pubs_arrays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indices "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the query set matches the training set, the nearest neighbor of each point is the point itself, at a distance of zero.\n",
    "It is also possible to efficiently produce a sparse graph showing the connections between neighboring points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbrs.kneighbors_graph(pp_pubs_arrays).toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id=\"Task 8\">Task 8</a>\n",
    "Can you get the nearest neighbors for London pubs by using KDTree? Your code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indices_kd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances_kd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are there significant differences spotted? For KDTree realization, you could also get the result through calling `KDTree` from `sklearn.neighbors`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another KDTree function\n",
    "from sklearn.neighbors import KDTree\n",
    "kdt = KDTree(pp_pubs_arrays, leaf_size=30, metric='euclidean')\n",
    "kdt.query(pp_pubs_arrays, k=2, return_distance=False)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Mean Nearest Neighbor Distance Statistics\">Mean Nearest Neighbor Distance Statistics</a>\n",
    "\n",
    "$$\\bar{d}_{min}=\\frac{1}{n} \\sum_{i=1}^n d_{min}(s_i)$$\n",
    "\n",
    "It is the average of all the points and corresponding distances demonstrated by Clark and Evans(1954), as a normal distribution under null hypothesis (CSR process). The distances between the nearest neighbor(s) $N(p)$ and the point $p$ is nearest neighbor distance for $p$, which meet the condition\n",
    "$$d_{p,N(p)} \\leq d_{p,j} \\forall j \\in S - p$$\n",
    "\n",
    "Hence the test statistics could be used to determine whether the point pattern is CSR, cluster or regular spatial process. For example, we can call the method `knn` to find $k$ nearest neighbors for each pub point in the point pattern *pp_pubs*. The first array is the most nearest neighbor for each pub point, the second array is the distance between each pub point and its nearest neighboring pub."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one nearest neighbor\n",
    "pp_pubs.knn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# two nearest neighbors\n",
    "pp_pubs.knn(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.max_nnd # Maximum nearest neighbor distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.min_nnd # Minimum nearest neighbor distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.mean_nnd # mean nearest neighbor distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_pubs.nnd # Nearest neighbor distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 9\">Task 9</a>\n",
    "Can you try to think of another way to calculate the mean nearest neighbor distance?\n",
    "**Hint**: recap on the definition of Mean, it should be the Sum divided by the total counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id=\"Nearest Neighbor Distance Functions\">Nearest Neighbor Distance Functions</a>\n",
    "\n",
    "Nearest neighbour distance distribution functions of a point process are cumulative distribution functions of several kinds -- $G, F, J$. By comparing the distance function of the observed point pattern with that of the point pattern from a CSR process, we are able to infer whether the underlying spatial process of the observed point pattern is CSR or not for a given confidence level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.spatial\n",
    "from pointpats.distance_statistics import G, F, J, K, L, Genv, Fenv, Jenv, Kenv, Lenv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $G$ function - event-to-event\n",
    "\n",
    "For a given distance $d$, $G(d)$ is the proportion of nearest neighbor distances that are less than $d$.\n",
    "$$G(d) = \\sum_{i=1}^n \\frac{ \\phi_i^d}{n}$$\n",
    "\n",
    "$$ \n",
    "\\phi_i^d =\n",
    " \\begin{cases}\n",
    "    1       & \\quad \\text{if } d_{min}(s_i)<d \\\\\n",
    "    0       & \\quad \\text{otherwise } \\\\\n",
    "  \\end{cases}\n",
    "$$\n",
    "\n",
    "If the underlying point process is a CSR process, $G$ function has an expectation of:\n",
    "$$\n",
    "G(d) = 1-e(-\\lambda \\pi d^2)\n",
    "$$\n",
    "If the $G$ function plot is above the expectation then it is a clustering distribution, whilst dispersive distribution with a curve below the expected plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp = G(pp, intervals=20)\n",
    "gp.plot() # plot empirical function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp1 = G(pp_pubs, intervals=20)\n",
    "gp1.plot() # plot empirical function "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orange curve is the function's expection, and the blue curve above is the real value, which indicating a clustering point pattern with its above position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.plot(qq=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp1.plot(qq=True) # quantile-quantile plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the q-q plot the CSR function is represented with a diagonal line. So obviously the $G$ increases very slowly at small distances and the line is above the expected value for a CSR process. We might think that the underlying spatial process is clustering. However, this visual inspection is not enough for a final conclusion.  In [Simulation Envelopes](#Simulation-Envelopes), we are going to demonstrate how to simulate data under CSR many times and construct the $95\\%$ simulation envelope for $G$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp1.d # distance domain sequence (corresponding to the x-axis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp1.G #cumulative nearest neighbor distance distribution over d (corresponding to the y-axis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $F$ function - \"point-event\" \n",
    "\n",
    "When the number of events in a point pattern is small, $G$ function is rough (see the $G$ function plot for the 12 size point pattern above). One way to get around this is to turn to $F$ funtion where a given number of randomly distributed points are generated in the domain and the nearest event neighbor distance is calculated for each point. The cumulative distribution of all nearest event neighbor distances is called $F$ function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1 = F(pp_pubs, intervals=100) # The default is to randomly generate 100 points.\n",
    "fp1.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1.plot(qq=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1 = F(pp_pubs, intervals=500) #We can increase the number of intervals to make F more smooth.\n",
    "fp1.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1.plot(qq=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Which function is comparatively smoother so far?\n",
    "$F$ function is more smooth than $G$ function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $J$ function - a combination of \"event-event\" and \"point-event\"\n",
    "\n",
    "$J$ function is defined as follows:\n",
    "\n",
    "$$J(d) = \\frac{1-G(d)}{1-F(d)}$$\n",
    "\n",
    "If $J(d)<1$, the underlying point process is a cluster point process; if $J(d)=1$, the underlying point process is a random point process; otherwise, it is a regular point process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jp1 = J(pp_pubs, intervals=100)\n",
    "jp1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above figure, we can observe that $J$ function is obviously below the $J(d)=1$ horizontal line. It is approaching finity with nearest neighbor distance increasing. We might tend to conclude that the underlying point process is a CSR one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Interevent Distance Functions\">Interevent Distance Functions</a>\n",
    "\n",
    "Nearest neighbor distance functions consider only the nearest neighbor distances, \"event-event\", \"point-event\" or the combination. Thus, distances to higer order neighbors are ignored, which might reveal important information regarding the point process. Interevent distance functions $K$ function, is proposed to consider distances between all pairs of event points. Similar to $G$, $F$ and $J$ functions, $K$ function is also cumulative distribution function.\n",
    "\n",
    "#### $K$ function - \"interevent\"\n",
    "\n",
    "Given distance $d$, $K(d)$ is defined as:\n",
    "$$K(d) = \\frac{\\sum_{i=1}^n \\sum_{j=1}^n \\psi_{ij}(d)}{n \\hat{\\lambda}}$$\n",
    "\n",
    "where\n",
    "$$ \n",
    "\\psi_{ij}(d) =\n",
    " \\begin{cases}\n",
    "    1       & \\quad \\text{if } d_{ij}<d \\\\\n",
    "    0       & \\quad \\text{otherwise } \\\\\n",
    "  \\end{cases}\n",
    "$$\n",
    "\n",
    "$\\sum_{j=1}^n \\psi_{ij}(d)$ is the number of events within a circle of radius $d$ centered on event $s_i$ .\n",
    "\n",
    "Still, we use CSR as the benchmark (null hypothesis) and see how the $K$ funtion estimated from the observed point pattern deviate from that under CSR, which is $K(d)=\\pi d^2$. $K(d)<\\pi d^2$ indicates that the underlying point process is a regular point process. $K(d)>\\pi d^2$ indicates that the underlying point process is a cluster point process. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kp1 = K(pp_pubs)\n",
    "kp1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is obvious that, as distance increases, the K function increases dramatically, indicating a more deviation from CSR process. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From above, we may conclude that the principle behind nearest neighbor methods is distance based (most commonly use Euclidean distance), to find a specified **number** of \"training\" samples closest to the target point (which is \"unseen\" in real case), and try to predict the label for the target from samples. So the **number** here could be a constant defined by us (for example, K nearest neighbors), or varied number based on local points density (for example, radius-based neighbors). So in the following section, we will get an idea on K Nearest Neighbors statistics, but leave the plot visualization part to next week.\n",
    "\n",
    "## <a id=\"K Nearest Neighbors (KNN)\">K Nearest Neighbors (KNN)</a>\n",
    "\n",
    "**KNN** is a non parametric and instance-based algorithm, and is normally used in a supervised learning setting for classification purpose. With a dataset of training observations $(x,y)$, where x denotes the feature and y denotes the target, we want to capture the relationship between x and y for prediction purpose.\n",
    "\n",
    "This algorithm is similarly defined according to a distance metric between two data points, which we are familiar with, and choosing the popular Euclidean distance by\n",
    "\n",
    "\n",
    "$$d(x,x') = \\sum^n_{i=1} \\sqrt{(x_1-x'_1)^2+(x_2-x'_2)^2+......+(x_n-x'_n)^2}$$\n",
    "\n",
    "\n",
    "The K-nearest neighbor algorithm  was mostly used for classification, given a positive integer $K$, to search for the major vote between the $K$ most similar instances to a given “unseen” observation $x$. It firstly runs through the whole dataset computing a similarity metric $d$ between $x$ and each training observation (the samples); a set consisting of $K$ closest points to $x$ will be get; then it estimates the conditional probability for each class, which is the proportion of points in the generated set with corresponding class label; at last, the observation $x$ gets assigned to the class with largest probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Weights\">Weights</a>\n",
    "\n",
    "PySAL weights can be constructed easily from many ways, i.e. array, shapefile and dataframe. For this practical, we are going to construct the weights directly from pandas and geopandas using the .from_shapefile and .from_dataframe method. Although we will still use points data for this practical, but it is better to know that if it is polygon data, then centroid points will be used by default.\n",
    "\n",
    "#### <a id=\"k-nearest neighbor weights\">k-nearest neighbor weights</a>\n",
    "\n",
    "The neighbors for a given observations can be defined using a k-nearest neighbor criterion, to create nearest neighbor weights matrix based on k nearest neighbors' distances. Based on our dataset, we can either build a knn weights directly from the london pubs shapefile, or implement it through dataframe. Please refer to http://pysal.readthedocs.io/en/latest/users/tutorials/weights.html#k-nearest-neighbor-weights for further information. Let's take $k=8$ as an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wknn_shp = ps.weights.KNN.from_shapefile(('data/london_pubs/london_pubs.shp'), k=8)\n",
    "wknn_shp.neighbors[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wdf = gpd.read_file('data/london_pubs/london_pubs.shp')\n",
    "wknn_df = ps.weights.KNN.from_dataframe(wdf, k=8)\n",
    "wknn_df.neighbors[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will get the identical series of 8 nearest neighbours from both functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <a id=\"Distance band weights\">Distance band weights</a>\n",
    "\n",
    "Knn weights ensure that all observations have the same number of neighbors. The distance bands weights define the neighbor set for each spatial unit as those other units falling within a threshold distance of the focal unit, so the number of neighbours is very likely to vary across observations with distance band weights. They can be generated for shapefiles and arrays of points with the starting point at determining the minimum nearest neighbour distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh = ps.min_threshold_dist_from_shapefile('data/london_pubs/london_pubs.shp')\n",
    "print (thresh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with this threshold in hand, the distance band weights are obtained as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wband = ps.weights.DistanceBand.from_shapefile('data/london_pubs/london_pubs.shp', threshold=thresh, binary=True)\n",
    "wband.min_neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Could you please explain the output here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wband.histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(wband.neighbors[0]) == set([1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(wband.neighbors[1]) == set([3,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wband = ps.weights.DistanceBand.from_dataframe(wdf, threshold=thresh, binary=True)\n",
    "wband.max_neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <a id=\"Kernel weights\">Kernel weights</a>\n",
    "\n",
    "This is a combination of distance-based thresholds together with continuously valued weights supported through kernel weights. The bandwidth attribute plays the role of the distance threshold with kernel weights, while the form of the kernel function determines the distance decay in the derived continuous weights (e.g.,‘triangular’,’uniform’,’quadratic’,’epanechnikov’,’quartic’,’bisquare’,’gaussian’). All kernel methods also support construction from shapefiles with Kernel.from_shapefile, and from dataframes with Kernel.from_dataframe.\n",
    "\n",
    "In the following example, we will use the default bandwidth setting as fixed across the observations.However you may set it into adaptive bandwidth by turning off the default 'fixed' upon your requirements.\n",
    "More details on kernel weights can be found in Kernel (http://pysal.readthedocs.io/en/latest/library/weights/Distance.html#pysal.weights.Distance.Kernel)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kw = ps.weights.Kernel.from_dataframe(wdf, fixed = True, function = 'gaussian')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kw.weights[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kw.bandwidth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"Task 10\">Task 10</a>\n",
    "\n",
    "Summarize the aforementioned 3 .from_dataframe methods below, and compare their outputs using the LDN-LSOAs dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "lsoadf = gpd.read_file('data/LDN-LSOAs.shp')\n",
    "???\n",
    "\n",
    "print(wknn_lsoa.neighbors[0])\n",
    "print(kw_lsoa.neighbors[0])\n",
    "print(dbb_lsoa.neighbors[0])\n",
    "print(dbc_lsoa.neighbors[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Credits\n",
    "You've done almost 99.9%! Let's have something Fun, to get a simple heatmap of London Pubs for your reference.\n",
    "\n",
    "To realize the heatmap function for today, we scheduled to use `GeoDataFrame` to get point data geometries, then to realize the visualization of heatmap by defining a heatmap function, which utilizes `Numpy`'s 2D historgram binning, smoothing from `Scipy` and `matplotlib`'s plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use geopandas to read london pubs shapefile data\n",
    "# your code here\n",
    "ldn_pubs=gpd.read_file('data/london_pubs/london_pubs.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define heatmap function\n",
    "\n",
    "from scipy import ndimage\n",
    "def heatmap(d, bins=(50,50), smoothing=2, cmap='jet'):\n",
    "    def getx(pt):\n",
    "        return pt.coords[0][0]     \n",
    "\n",
    "    def gety(pt):\n",
    "        return pt.coords[0][1]\n",
    "\n",
    "    x = list(d.geometry.apply(getx)) # get the x value from geometry, replace ???\n",
    "    y = list(d.geometry.apply(gety)) # get the y value from geometry, replace ???\n",
    "    heatmap, xedges, yedges = np.histogram2d(y, x, bins=bins) # call numpy's 2D histogram function\n",
    "    extent = [yedges[0], yedges[-1], xedges[-1], xedges[0]]\n",
    "\n",
    "    logheatmap = np.log(heatmap)\n",
    "    logheatmap[np.isneginf(logheatmap)] = 0\n",
    "    logheatmap = ndimage.filters.gaussian_filter(logheatmap, smoothing, mode='nearest') # scipy's ndimage\n",
    "    \n",
    "    plt.imshow(logheatmap, cmap=cmap, extent=extent)  # plot the heatmap\n",
    "    plt.colorbar()\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(ldn_pubs, bins=50, smoothing=2)\n",
    "# try change the value of bins and smoothing, what are the differences?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next week, we will mainly look at density based point pattern analysis (density estimates, KDE, etc.), and Nearest Neighbors Classification realization by using London pubs data and Airbnb listings data, so please get them at hand then.\n",
    "\n",
    "## Take home task (Optional):\n",
    "\n",
    "This is handy for you to understand what is heatmap, but the map quality is still not good, do you have any idea in mind for better visualization? We may then use it for next week. \n",
    "\n",
    "**Hint: Interactive Map we've learned!**\n",
    "\n",
    "Since we are going to realize KDE heatmap visualization next week, it is better to have an idea of KDE utilization in the GIS application you are familiar with. So if you are interested to have a try, it is recommended to realize the KDE Heatmap visualization for London Pubs data in ArcGIS or QGIS for visualization. Try to compare these methods and interpret what the heatmap telling us about?\n",
    "\n",
    "For ArcGIS user: \n",
    "-  in ArcMap, then go to find the Spatial Analyst Tools -> Density -> Kernel Density and carry out KDE under ArcToolBox.\n",
    "\n",
    "For QGIS user:\n",
    "- find Heatmap Plugin, or download and install the plug-in from (http://docs.qgis.org/2.0/en/docs/user_manual/plugins/plugins_heatmap.html)\n",
    "- View -> Toolbars -> Raster -> Heatmap\n",
    "- set up radius, further select colormap from Layer properties."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 1
}