CW2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a02b62c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import osmnx as ox\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import spaghetti\n",
    "import os\n",
    "import glob\n",
    "import networkx as nx\n",
    "from matplotlib import pyplot as plt\n",
    "from pyproj import Proj\n",
    "import smopy\n",
    "import collections\n",
    "import numpy as np\n",
    "import math\n",
    "import esda\n",
    "from haversine import haversine\n",
    "import warnings\n",
    "from shapely.geometry import Point, LineString\n",
    "from haversine import haversine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe6e2a51",
   "metadata": {},
   "source": [
    "### Select city area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9ad510c",
   "metadata": {},
   "outputs": [],
   "source": [
    "location_researched = \"Centre of Leeds, United Kingdom\"\n",
    "# buffer_dist is to specifies the distance in meters to buffer the geocoded point. \n",
    "leeds_centre = ox.geocode_to_gdf(location_researched, buffer_dist=1000)\n",
    "leeds_polygon = leeds_centre['geometry'][0]\n",
    "leeds_polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2173be24",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# get the position of the Leeds City Centre which we choose to be the original area.\n",
    "# leeds_centre.lat.values,leeds_centre.lon.values\n",
    "# The longitudes and latitudes of the central points are：\n",
    "print(\"Latitude of the central points is:\",leeds_centre.lat.values)\n",
    "print(\"longitude of the central points is:\",leeds_centre.lon.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5c0eb59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The years we used to generate the network are from 2014-2019\n",
    "year_used = ['2014','2015','2016','2017','2018','2019']# '2014','2015','2016',\n",
    "leeds_stored = []\n",
    "for i in range(len(year_used)):\n",
    "    # ISO-8859-1 encoding is a single-byte encoding, backwards compatible with ASCII, \n",
    "    # and is the encoding standard used in many European countries.\n",
    "    leeds_flag = pd.read_csv('./'+year_used[i] +'.csv',encoding='ISO-8859-1')\n",
    "    leeds_stored.append(leeds_flag)\n",
    "leeds_df = pd.concat([k for k in leeds_stored], axis=0, ignore_index=True)\n",
    "leeds_df.head(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7fb599",
   "metadata": {},
   "outputs": [],
   "source": [
    "leeds_df.dropna(subset=['Grid Ref: Easting', 'Grid Ref: Northing'], inplace=True)\n",
    "leeds_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1abece30",
   "metadata": {},
   "outputs": [],
   "source": [
    "# geometry=gpd.points_from_xy(leeds_df[\"Grid Ref: Easting\"], leeds_df[\"Grid Ref: Northing\"],crs=\"EPSG:27700\")\n",
    "leeds_gdf = gpd.GeoDataFrame(leeds_df, geometry=gpd.points_from_xy(leeds_df[\"Grid Ref: Easting\"], leeds_df[\"Grid Ref: Northing\"],crs=\"EPSG:27700\"))\n",
    "leeds_new_gdf=leeds_gdf.to_crs(epsg=4326).drop(['Grid Ref: Northing','Grid Ref: Easting'],axis=1)\n",
    "leeds_new_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27d336e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "accident_nodes = leeds_new_gdf.geometry\n",
    "# Select the nodes that are in the pre-determined range\n",
    "leeds_accident = leeds_new_gdf[accident_nodes.within(leeds_polygon)]\n",
    "len(leeds_accident)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "908653a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# in order to show the result of the accident points with the defined area, I draw them together.\n",
    "tt = leeds_accident.geometry\n",
    "box_stations = (tt.y.min(),tt.x.min(),\n",
    "       tt.y.max(), tt.x.max())\n",
    "map = smopy.Map(box_stations, z=100) # z is the zoom level, here we set it as 10, you may want to try another value.\n",
    "# map.show_ipython()\n",
    "p1 = pd.Series(leeds_polygon.exterior.coords.xy[1].tolist())\n",
    "p2 = pd.Series(leeds_polygon.exterior.coords.xy[0].tolist())\n",
    "x1, y1 = map.to_pixels(tt.y, tt.x)\n",
    "x2, y2 = map.to_pixels(p1,p2)\n",
    "ax = map.show_mpl(figsize=(12, 10))\n",
    "ax.plot(x1, y1, 'o',color='crimson',ms=4, mew=1) # 'or' means red dots/circles\n",
    "ax.plot(x2, y2, color='darkblue', linestyle='solid',linewidth=4, markersize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01bef5ab",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Download the road network data\n",
    "leeds_graph = ox.graph_from_polygon(leeds_polygon, network_type='drive')\n",
    "nx.check_planarity(leeds_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f461597d",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_nodes = leeds_graph.number_of_nodes()\n",
    "num_edges = leeds_graph.number_of_edges()\n",
    "print(\"The number of nodes in this road network is: \",num_nodes)\n",
    "print(\"The number of edges in this road network is: \",num_edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67571a0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "leeds_projection = ox.project_gdf(leeds_centre).unary_union.area\n",
    "ox.basic_stats(leeds_graph, area=leeds_projection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3489a431",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_shortest_paths = dict(nx.all_pairs_dijkstra_path_length(G, weight='length'))\n",
    "\n",
    "# Cast the keys to integers\n",
    "all_shortest_paths = {int(k):v for k,v in all_shortest_paths.items()}\n",
    "\n",
    "# Find the longest of these shortest paths\n",
    "longest_shortest_path = max([max(path_lengths.values()) for path_lengths in all_shortest_paths.values()])\n",
    "\n",
    "print(\"The spatial diameter is same as the longest shortest path in the network, which is: {:.2f} meters\".format(longest_shortest_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ebc73c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.check_planarity(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1cf6239",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_values = nx.get_node_attributes(leeds_graph, 'x')\n",
    "y_values = nx.get_node_attributes(leeds_graph, 'y')\n",
    "\n",
    "# We will work with the edges and add the missing geometries (lines denoting the roads between points)\n",
    "graph_with_geometries = list(leeds_graph.edges(data=True))\n",
    "\n",
    "# Iterate through the edges and, where missing, add a geometry attribute with the line between start and end nodes\n",
    "for e in graph_with_geometries:\n",
    "    if not 'geometry' in e[2]:\n",
    "        e[2]['geometry'] = LineString([\n",
    "            Point(x_values[e[0]], y_values[e[0]]),\n",
    "            Point(x_values[e[1]], y_values[e[1]])])\n",
    "\n",
    "graph_with_geometries[0:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09afb3c9",
   "metadata": {},
   "source": [
    "### Task2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59f31d82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will now drop the start and end nodes, as we will construct a new Spaghetti network based on the geometries of the roads\n",
    "road_lines = [x[2] for x in graph_with_geometries]\n",
    "# From this, we can construct a GeoDataFrame\n",
    "roads_geodataframe = gpd.GeoDataFrame(pd.DataFrame(road_lines))\n",
    "roads_geodataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0c15e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# From the GeoDataFrame, we can construct a network in Spaghetti from which to do point analysis\n",
    "leeds_points_graph = spaghetti.Network(in_data=roads_geodataframe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ade22b5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Check what this new network looks like by getting DataFrames for the edges and plotting them\n",
    "\n",
    "nodes_df, edges_df = spaghetti.element_as_gdf(leeds_points_graph, vertices=True, arcs=True)\n",
    "base_network = edges_df.plot(color=\"k\", zorder=0, figsize=(15, 15))\n",
    "nodes_df.plot(ax=base_network, color=\"r\", zorder=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25e2c67e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will now snap the leeds located accidents we extracted earlier, i.e. position them at the closest point on the closest road\n",
    "leeds_points_graph.snapobservations(leeds_accident, 'accidents')\n",
    "\n",
    "# We can see the difference between the original accident coordinates and their position when snapped to the road network\n",
    "print(\"observation 1\\ntrue coords:\\t%s\\nsnapped coords:\\t%s\" % (\n",
    "    leeds_points_graph.pointpatterns[\"accidents\"].points[0][\"coordinates\"],\n",
    "    leeds_points_graph.pointpatterns[\"accidents\"].snapped_coordinates[0]\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41e8c6ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the network\n",
    "base_network = edges_df.plot(color=\"k\", zorder=0, figsize =(12, 12))\n",
    "# Get a GeoDataFrame of the snapped accident locations to plot on the network image\n",
    "snapped_accidents=spaghetti.element_as_gdf(\n",
    "    leeds_points_graph, pp_name='accidents', snapped=True)\n",
    "\n",
    "# Plot these on the road network\n",
    "snapped_accidents.plot(\n",
    "    color=\"r\", marker=\"x\",\n",
    "    markersize=50, zorder=1, ax=base_network)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e0133b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the network\n",
    "base_network = edges_df.plot(color=\"k\", zorder=0, figsize =(12, 12))\n",
    "# Get a GeoDataFrame of the non-snapped (real) crime locations to plot on the net\n",
    "observed_accidents=spaghetti.element_as_gdf(\n",
    "    leeds_points_graph, pp_name='accidents', snapped=False)\n",
    "\n",
    "# Plot these on the road network\n",
    "observed_accidents.plot(\n",
    "    color=\"r\", marker=\"x\",\n",
    "    markersize=50, zorder=1, ax=base_network)\n",
    "\n",
    "plt.savefig('accidents-observations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "140747cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "kres = leeds_points_graph.GlobalAutoK(\n",
    "    leeds_points_graph.pointpatterns[\"accidents\"],\n",
    "    nsteps=100, permutations=100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30f1386b",
   "metadata": {},
   "outputs": [],
   "source": [
    "kres.lam\n",
    "kres.xaxis\n",
    "kres.observed\n",
    "kres.upperenvelope\n",
    "kres.lowerenvelope\n",
    "kres.sim\n",
    "\n",
    "print(f\"Density of points in the network (lambda): {kres.lam}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8985612",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Distances at which density is measured:\\n{kres.xaxis}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a4c8db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(kres.xaxis, kres.observed, \"b-\", label=\"Observed\")\n",
    "ax.plot(kres.xaxis, kres.upperenvelope, \"r--\", label=\"Upper\")\n",
    "ax.plot(kres.xaxis, kres.lowerenvelope, \"k--\", label=\"Lower\")\n",
    "\n",
    "ax.legend(loc=\"best\", fontsize=\"x-large\")\n",
    "ax.set_xlabel(\"Distance $(r)$\")\n",
    "ax.set_ylabel(\"$K(r)$\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad296a1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get snapped point pattern \n",
    "pointpat = leeds_points_graph.pointpatterns['accidents']\n",
    "# Get count of points per network edge: a dictionary from each edge to the crime count on that edge\n",
    "counts = leeds_points_graph.count_per_link(pointpat.obs_to_arc, graph=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e08ce315",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9c7cb14",
   "metadata": {},
   "outputs": [],
   "source": [
    "moran, yaxis_moran = leeds_points_graph.Moran('accidents',graph=True)\n",
    "moran.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58db35ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "moran, yaxis_moran = leeds_points_graph.Moran(\"accidents\")\n",
    "moran.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69101273",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the weights matrix for edges in the graph (just the adjacency matrix with 1 where edges connect at a node, 0 otherwise)\n",
    "weights = leeds_points_graph.w_network\n",
    "# Get the edges included in the weights matrix: an enumerator for a list of edges\n",
    "edges = weights.neighbors.keys()\n",
    "# Construct an array of the counts values per edge in the same order as\n",
    "# the weights matrix, with 0.0 where no counts recorded\n",
    "values = [counts[edge] if edge in counts.keys () else 0. \\\n",
    "    for index, edge in enumerate(edges)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddce7b01",
   "metadata": {},
   "outputs": [],
   "source": [
    "moran = esda.moran.Moran(values, weights)\n",
    "moran.I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19cb2fe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "moran.p_sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b76973c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "from splot.esda import moran_scatterplot, lisa_cluster, plot_moran\n",
    "figsize = (20,10)\n",
    "plot_moran(moran, zstandard=True, scatter_kwds=None, figsize=figsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "134dc88c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(moran.sim, shade=True)\n",
    "plt.vlines(moran.I, 0, 1, color='r')\n",
    "plt.vlines(moran.EI, 0,1)\n",
    "plt.xlabel(\"Moran's I\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdbd2e8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "snapped_accidents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bb67bf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_geodataframe.geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e76c22f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = snapped_accidents.geometry.x\n",
    "Y = snapped_accidents.geometry.y\n",
    "dist_point_street = ox.distance.nearest_edges(leeds_graph, X, Y, interpolate=None, return_dist=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13ac7c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_point_street"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dba68aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_accidents = len(dist_point_street)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b076210",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = [haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][0]]['y'],leeds_graph.nodes()[dist_point_street[i][0]]['x']), unit='m') \n",
    "        if haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][0]]['y'],leeds_graph.nodes()[dist_point_street[i][0]]['x']), unit='m') \n",
    "           <= haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][1]]['y'],leeds_graph.nodes()[dist_point_street[i][1]]['x']), unit='m')\n",
    "        else haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][1]]['y'],leeds_graph.nodes()[dist_point_street[i][1]]['x']), unit='m') \n",
    "        for i in range(len(dist_point_street))]\n",
    "percentage = [dist[i] / haversine((leeds_graph.nodes()[dist_point_street[i][1]]['y'],leeds_graph.nodes()[dist_point_street[i][1]]['x']), \n",
    "                              (leeds_graph.nodes()[dist_point_street[i][0]]['y'],leeds_graph.nodes()[dist_point_street[i][0]]['x']), unit='m')\n",
    "            if haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][0]]['y'],leeds_graph.nodes()[dist_point_street[i][0]]['x']), unit='m') \n",
    "               <= haversine((Y[i],X[i]), (leeds_graph.nodes()[dist_point_street[i][1]]['y'],leeds_graph.nodes()[dist_point_street[i][1]]['x']), unit='m')\n",
    "            else dist[i] / haversine((leeds_graph.nodes()[dist_point_street[i][1]]['y'],leeds_graph.nodes()[dist_point_street[i][1]]['x']), \n",
    "                                     (leeds_graph.nodes()[dist_point_street[i][0]]['y'],leeds_graph.nodes()[dist_point_street[i][0]]['x']), unit='m')\n",
    "            for i in range(len(dist_point_street))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29ceb994",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e67c5875",
   "metadata": {},
   "outputs": [],
   "source": [
    "ave_distance = sum(dist)/num_accidents\n",
    "ave_percentage = sum(percentage)/num_accidents\n",
    "print('The average distance between accident point and the intersection is: ',round(ave_distance,3),' meters') \n",
    "print('The average position percentage to intersection is: ',round(ave_percentage*100,3),' %')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0def8155",
   "metadata": {},
   "source": [
    "### Task3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b8f0961",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e40c6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.cluster import KMeans\n",
    "import distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03a8299b",
   "metadata": {},
   "outputs": [],
   "source": [
    "PLACE_OF_INTEREST = \"Leeds, United Kingdom\"\n",
    "# Get the network of streets in Soho\n",
    "leeds_graph = ox.graph_from_place(PLACE_OF_INTEREST,network_type='drive')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09b43ecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.plot_graph(leeds_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e0286e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can see that the nodes are specified by x, y, coordinates\n",
    "list(leeds_graph.nodes(data=True))[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b20bd4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# And the edges are tuples of start node, end node and a dictionary with attributes including geometry\n",
    "list(leeds_graph.edges(data=True))[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9731044",
   "metadata": {},
   "source": [
    "### Even distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31159a21",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "NUMBER_OF_SEEDS = 10\n",
    "nodes = list(leeds_graph.nodes)\n",
    "\n",
    "# Extract x,y coordinates of each node\n",
    "coordinates = [(leeds_graph.nodes[n]['x'], leeds_graph.nodes[n]['y']) for n in nodes]\n",
    "\n",
    "# Use k-means++ to get initial seed locations\n",
    "kmeans = KMeans(n_clusters=NUMBER_OF_SEEDS, init='k-means++').fit(coordinates)\n",
    "seeds = kmeans.cluster_centers_.tolist()\n",
    "all_nodes = list(leeds_graph.nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef68399",
   "metadata": {},
   "outputs": [],
   "source": [
    "seeds = [ox.distance.nearest_nodes(leeds_graph, longitude, latidude) for longitude, latidude in seeds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7450540",
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = {seed: nx.single_source_dijkstra_path_length(\n",
    "    leeds_graph, seed, weight='length') for seed in seeds}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b3d67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nearest_from_list(node_distances):\n",
    "    return sorted(node_distances, key=lambda node_length: node_length[1])[0] \\\n",
    "         if len(node_distances) > 0 else None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9845d05b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nearest_seed(node):\n",
    "    seed_distances = [(seed, distances[seed][node]) \\\n",
    "        for seed in seeds if node in distances[seed]]\n",
    "    return nearest_from_list(seed_distances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cdfae84",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nearest_for_edge(edge):\n",
    "    nearest_to_ends_all = [nearest_seed(edge[0]), nearest_seed(edge[1])]\n",
    "    nearest_to_ends = [distance for distance in nearest_to_ends_all if distance]\n",
    "    return nearest_from_list(nearest_to_ends)\n",
    "colours = ox.plot.get_colors(NUMBER_OF_SEEDS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52135237",
   "metadata": {},
   "outputs": [],
   "source": [
    "def colour_for_seed_distance(seed):\n",
    "    return colours[seeds.index(seed[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9c20538",
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_nearest_seeds = [nearest_for_edge(edge) for edge in leeds_graph.edges]\n",
    "# Note that edges not connected to a seed shown in black, so invisible on black background\n",
    "edge_colours = [colour_for_seed_distance(seed) if seed else 'k' for seed in edge_nearest_seeds]\n",
    "# For the road network nodes, we want the seeds to be coloured red and the non-seed nodes to be coloured white.\n",
    "node_colours = ['r' if node in seeds else 'none' for node in all_nodes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c9f7316",
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.plot.plot_graph(leeds_graph, edge_color = edge_colours, node_color = node_colours, bgcolor = 'k', save = True, filepath = 'nvd.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21585db1",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_subgraphs = {}\n",
    "for seed in seeds:\n",
    "    seed_nodes = [node for node in list(leeds_graph.nodes()) if nearest_seed(node)!=None and nearest_seed(node)[0] == seed]\n",
    "    seed_edges = [edge for edge in list(leeds_graph.edges()) if nearest_for_edge(edge)!=None and nearest_for_edge(edge)[0] == seed]\n",
    "    seed_subgraphs[seed] = leeds_graph.subgraph(seed_nodes + [n for e in seed_edges for n in e])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bff2504",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_subgraphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b492bff",
   "metadata": {},
   "outputs": [],
   "source": [
    "for seed, subgraph in seed_subgraphs.items():\n",
    "    print(f\"Seed {seed}:\")\n",
    "    print(f\"Number of nodes: {subgraph.number_of_nodes()}\")\n",
    "    print(f\"Number of edges: {subgraph.number_of_edges()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b8be1fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.plot.plot_graph(seed_subgraphs[393344811])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c0e3d36",
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.plot.plot_graph(seed_subgraphs[384962202])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4467733",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = seed_subgraphs[1955838372]\n",
    "H = nx.Graph(G)  # 将多重图 G 转换成普通图 H\n",
    "cycles = nx.cycle_basis(H) \n",
    "for cycle in cycles:\n",
    "    subgraph = H.subgraph(cycle)\n",
    "    length = sum([subgraph[u][v]['length'] for u, v in subgraph.edges()])\n",
    "    print(f\"Cycle: {cycle}, Length: {length}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aafcb36",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cycle in cycles:\n",
    "    temp_subgraph = H.subgraph(cycle)\n",
    "    length = sum([temp_subgraph[u][v]['length'] for u, v in temp_subgraph.edges()])\n",
    "    if length<=44000 and length>=40000:\n",
    "        # print(n)\n",
    "        print(\"The length of the cycle is: \",length)\n",
    "        break\n",
    "print(temp_subgraph)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df26d719",
   "metadata": {},
   "source": [
    "### seed_subgraphs[393344811]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66f85c6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = seed_subgraphs[393344811]\n",
    "H = nx.Graph(G)  # 将多重图 G 转换成普通图 H\n",
    "cycles = nx.cycle_basis(H) \n",
    "for cycle in cycles:\n",
    "    subgraph = H.subgraph(cycle)\n",
    "    length = sum([subgraph[u][v]['length'] for u, v in subgraph.edges()])\n",
    "    print(f\"Cycle: {cycle}, Length: {length}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "345031ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cycle in cycles:\n",
    "    temp_subgraph = H.subgraph(cycle)\n",
    "    length = sum([temp_subgraph[u][v]['length'] for u, v in temp_subgraph.edges()])\n",
    "    if length<=42500 and length>=41500:\n",
    "        # print(n)\n",
    "        print(\"The length of the cycle is: \",length)\n",
    "        break\n",
    "print(temp_subgraph)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11ff3bff",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(250,250))\n",
    "pos = nx.spring_layout(temp_subgraph, seed=42)\n",
    "# nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='r')\n",
    "\n",
    "nx.draw_networkx_nodes(temp_subgraph, pos=pos, node_color='blue', node_size=500)\n",
    "nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='red', width=5)\n",
    "\n",
    "plt.title(f\"Cycle: {cycle}, Length: {length}\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "930e1ebd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(350, 350))\n",
    "\n",
    "# draw the entire graph H\n",
    "pos = nx.spring_layout(H, seed=42)\n",
    "nx.draw_networkx_nodes(H, pos=pos, node_color='blue', node_size=500)\n",
    "nx.draw_networkx_edges(H, pos=pos, edge_color='gray', width=5)\n",
    "\n",
    "# draw the subgraph temp_subgraph on top of the entire graph\n",
    "nx.draw_networkx_nodes(temp_subgraph, pos=pos, node_color='red', node_size=500)\n",
    "nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='red', width=10)\n",
    "\n",
    "plt.title(f\"Cycle: {cycle}, Length: {length}\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8c8326c",
   "metadata": {},
   "source": [
    "### seed_subgraphs[381915141]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "401e971c",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = seed_subgraphs[381915141]\n",
    "H = nx.Graph(G)  # 将多重图 G 转换成普通图 H\n",
    "cycles = nx.cycle_basis(H) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d78ba497",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cycle in cycles:\n",
    "    subgraph = H.subgraph(cycle)\n",
    "    length = sum([subgraph[u][v]['length'] for u, v in subgraph.edges()])\n",
    "    print(f\"Cycle: {cycle}, Length: {length}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "317c7bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "n=0\n",
    "for cycle in cycles:\n",
    "    n=n+1\n",
    "    temp_subgraph = H.subgraph(cycle)\n",
    "    length = sum([temp_subgraph[u][v]['length'] for u, v in temp_subgraph.edges()])\n",
    "    if length<=42500 and length>=41500:\n",
    "        print(\"The length of the cycle is: \",length)\n",
    "        break\n",
    "    else: \n",
    "        temp_subgraph == None\n",
    "        length == 0\n",
    "        \n",
    "print(temp_subgraph)\n",
    "print(length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "659bcbe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(250,250))\n",
    "pos = nx.spring_layout(temp_subgraph, seed=42)\n",
    "# nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='r')\n",
    "\n",
    "nx.draw_networkx_nodes(temp_subgraph, pos=pos, node_color='blue', node_size=500)\n",
    "nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='red', width=5)\n",
    "\n",
    "plt.title(f\"Cycle: {cycle}, Length: {length}\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "447b2dd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(350, 350))\n",
    "\n",
    "# draw the entire graph H\n",
    "pos = nx.spring_layout(H, seed=42)\n",
    "nx.draw_networkx_nodes(H, pos=pos, node_color='blue', node_size=500)\n",
    "nx.draw_networkx_edges(H, pos=pos, edge_color='gray', width=5)\n",
    "\n",
    "# draw the subgraph temp_subgraph on top of the entire graph\n",
    "nx.draw_networkx_nodes(temp_subgraph, pos=pos, node_color='red', node_size=500)\n",
    "nx.draw_networkx_edges(temp_subgraph, pos=pos, edge_color='red', width=10)\n",
    "\n",
    "plt.title(f\"Cycle: {cycle}, Length: {length}\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94db62a",
   "metadata": {},
   "source": [
    "### Question 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2550473d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "NUMBER_OF_SEEDS = 4\n",
    "nodes = list(leeds_graph.nodes)\n",
    "\n",
    "# Set the random seed\n",
    "np.random.seed(76)\n",
    "\n",
    "# Extract x,y coordinates of each node\n",
    "coordinates = [(leeds_graph.nodes[n]['x'], leeds_graph.nodes[n]['y']) for n in nodes]\n",
    "\n",
    "# Use k-means++ to get initial seed locations\n",
    "kmeans = KMeans(n_clusters=NUMBER_OF_SEEDS, init='k-means++').fit(coordinates)\n",
    "seeds = kmeans.cluster_centers_.tolist()\n",
    "all_nodes = list(leeds_graph.nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "822977f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "seeds = [ox.distance.nearest_nodes(leeds_graph, longitude, latidude) for longitude, latidude in seeds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a950562f",
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = {seed: nx.single_source_dijkstra_path_length(\n",
    "    leeds_graph, seed, weight='length') for seed in seeds}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de48ff4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nearest_from_list(node_distances):\n",
    "    return sorted(node_distances, key=lambda node_length: node_length[1])[0] \\\n",
    "         if len(node_distances) > 0 else None\n",
    "def nearest_seed(node):\n",
    "    seed_distances = [(seed, distances[seed][node]) \\\n",
    "        for seed in seeds if node in distances[seed]]\n",
    "    return nearest_from_list(seed_distances)\n",
    "\n",
    "def nearest_for_edge(edge):\n",
    "    nearest_to_ends_all = [nearest_seed(edge[0]), nearest_seed(edge[1])]\n",
    "    nearest_to_ends = [distance for distance in nearest_to_ends_all if distance]\n",
    "    return nearest_from_list(nearest_to_ends)\n",
    "colours = ox.plot.get_colors(NUMBER_OF_SEEDS)\n",
    "\n",
    "def colour_for_seed_distance(seed):\n",
    "    return colours[seeds.index(seed[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d58ebe59",
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_nearest_seeds = [nearest_for_edge(edge) for edge in leeds_graph.edges]\n",
    "# Note that edges not connected to a seed shown in black, so invisible on black background\n",
    "edge_colours = [colour_for_seed_distance(seed) if seed else 'k' for seed in edge_nearest_seeds]\n",
    "# For the road network nodes, we want the seeds to be coloured red and the non-seed nodes to be coloured white.\n",
    "node_colours = ['r' if node in seeds else 'none' for node in all_nodes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa785023",
   "metadata": {},
   "outputs": [],
   "source": [
    "ox.plot.plot_graph(leeds_graph, edge_color = edge_colours, node_color = node_colours, bgcolor = 'k', save = True, filepath = 'nvd.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96d152a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_subgraphs = {}\n",
    "for seed in seeds:\n",
    "    seed_nodes = [node for node in list(leeds_graph.nodes()) if nearest_seed(node)!=None and nearest_seed(node)[0] == seed]\n",
    "    seed_edges = [edge for edge in list(leeds_graph.edges()) if nearest_for_edge(edge)!=None and nearest_for_edge(edge)[0] == seed]\n",
    "    seed_subgraphs[seed] = leeds_graph.subgraph(seed_nodes + [n for e in seed_edges for n in e])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "138cfa70",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_subgraphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd245e72",
   "metadata": {},
   "outputs": [],
   "source": [
    "for seed, subgraph in seed_subgraphs.items():\n",
    "    print(f\"Seed {seed}:\")\n",
    "    print(f\"Number of nodes: {subgraph.number_of_nodes()}\")\n",
    "    print(f\"Number of edges: {subgraph.number_of_edges()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f78aa96f",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = seed_subgraphs[1984729084]\n",
    "H = nx.Graph(G)  # 将多重图 G 转换成普通图 H\n",
    "cycles = nx.cycle_basis(H) \n",
    "for cycle in cycles:\n",
    "    subgraph = H.subgraph(cycle)\n",
    "    length = sum([subgraph[u][v]['length'] for u, v in subgraph.edges()])\n",
    "    print(f\"Cycle: {cycle}, Length: {length}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a5e7e86",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cycle in cycles:\n",
    "    temp_subgraph = H.subgraph(cycle)\n",
    "    length = sum([temp_subgraph[u][v]['length'] for u, v in temp_subgraph.edges()])\n",
    "    if length<=42500 and length>=41500:\n",
    "        # print(n)\n",
    "        print(\"The length of the cycle is: \",length)\n",
    "        break\n",
    "print(temp_subgraph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "348d6945",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b62afb09",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}