clip_bb_to_g.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from shapely import Point, LineString\n",
    "import osmnx as ox\n",
    "from shapely.geometry import box, Polygon, MultiPolygon\n",
    "from shapely.ops import transform\n",
    "from functools import partial\n",
    "import pyproj\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_35260/2575286281.py:12: DtypeWarning: Columns (9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  lookup = pd.read_csv('data/PCD11_OA11_LSOA11_MSOA11_LAD11_EW_LU_aligned_v2.csv')\n"
     ]
    }
   ],
   "source": [
    "# Import MSOA lookup\n",
    "msoas = gpd.read_file('data/MSOA_EngWal_Dec_2011_Generalised_ClippedEW_0/Middle_Layer_Super_Output_Areas_December_2011_Generalised_Clipped_Boundaries_in_England_and_Wales.shp').to_crs(4326).set_index('msoa11cd')\n",
    "\n",
    "# Import MSOA 2011 OD data\n",
    "od_data = pd.read_parquet('data/od_2011.parquet')\n",
    "\n",
    "#Import LSOAs\n",
    "lsoas = gpd.read_file('data/LSOA_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V4_6029841263726194941.gpkg').to_crs(4326)\n",
    "lsoas = pd.concat([lsoas, lsoas.bounds], axis=1)\n",
    "\n",
    "#Import lsoa to msoa look up\n",
    "lookup = pd.read_csv('data/PCD11_OA11_LSOA11_MSOA11_LAD11_EW_LU_aligned_v2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read geojson : 521.9917948246002\n",
      "Read Get Nodes : 156.08479642868042\n",
      "Read Get Edges : 9.695680141448975\n",
      "Construct Network : 88.07392525672913\n"
     ]
    }
   ],
   "source": [
    "# Read in 2016\n",
    "\n",
    "t0 = time.time()\n",
    "G2016 = gpd.read_file('osm2016.geojson')\n",
    "t1 = time.time()\n",
    "print('Read geojson : {}'.format(t1 - t0))\n",
    "\n",
    "t0 = time.time()\n",
    "def extract_coordinates(line_string, line_index, coord_index):\n",
    "    #line = LineString(eval(line_string))\n",
    "    return list(line_string.coords)[line_index][coord_index]\n",
    "\n",
    "# Apply the function to the DataFrame column\n",
    "G2016['coordinates_u_x'] = G2016['geometry'].apply(extract_coordinates,line_index=0,coord_index=0)\n",
    "G2016['coordinates_u_y'] = G2016['geometry'].apply(extract_coordinates,line_index=0,coord_index=1)\n",
    "\n",
    "G2016['coordinates_v_x'] = G2016['geometry'].apply(extract_coordinates,line_index=-1,coord_index=0)\n",
    "G2016['coordinates_v_y'] = G2016['geometry'].apply(extract_coordinates,line_index=-1,coord_index=1)\n",
    "\n",
    "node_set_1 = G2016[['node1','coordinates_u_x','coordinates_u_y']].drop_duplicates(subset='node1')\n",
    "node_set_2 = G2016[['node2','coordinates_v_x','coordinates_v_y']].drop_duplicates(subset='node2')\n",
    "\n",
    "node_set_1 = node_set_1.rename(columns={'node1': 'osmid','coordinates_u_x': 'x','coordinates_u_y': 'y'})\n",
    "node_set_2 = node_set_2.rename(columns={'node2': 'osmid','coordinates_v_x': 'x','coordinates_v_y': 'y'})\n",
    "\n",
    "node_df = pd.concat([node_set_1,node_set_2]).drop_duplicates(subset='osmid').set_index('osmid')\n",
    "t1 = time.time()\n",
    "print('Read Get Nodes : {}'.format(t1 - t0))\n",
    "\n",
    "del(node_set_1)\n",
    "del(node_set_2)\n",
    "\n",
    "t0 = time.time()\n",
    "G2016['key'] = 0\n",
    "edge_gdf = gpd.GeoDataFrame(G2016[['node1','node2','key','geometry','length','lts']].set_index(['node1','node2','key']),geometry=G2016[['node1','node2','key','geometry','length','lts']].set_index(['node1','node2','key'])['geometry'])\n",
    "t1 = time.time()\n",
    "print('Read Get Edges : {}'.format(t1 - t0))\n",
    "\n",
    "del(G2016)\n",
    "\n",
    "t0 = time.time()\n",
    "G = ox.graph_from_gdfs(node_df,edge_gdf,{'crs': 'epsg:4326'})\n",
    "t1 = time.time()\n",
    "print('Construct Network : {}'.format(t1 - t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = ox.graph_to_gdfs(G, edges=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbx_expansion = 0.5\n",
    "\n",
    "#Below function from ChatGPT\n",
    "#Get expanded network - method 1 km buffer\n",
    "def expand_bbox(original_bbox, expansion_distance_km=5):\n",
    "    # Create a Shapely geometry object for the original bounding box\n",
    "    original_geometry = box(*original_bbox)\n",
    "    # Define a function to project the geometry to a new coordinate reference system\n",
    "    project = partial(\n",
    "        pyproj.transform,\n",
    "        pyproj.Proj(init='epsg:4326'),  # WGS 84 coordinate reference system\n",
    "        pyproj.Proj(proj='utm', zone=33, ellps='WGS84')  # Example: UTM Zone 33\n",
    "    )\n",
    "    # Project the original geometry to the new coordinate reference system\n",
    "    projected_geometry = transform(project, original_geometry)\n",
    "    # Calculate the expansion distance in the projected coordinate system\n",
    "    expansion_distance_meters = expansion_distance_km * 1000\n",
    "    # Expand the geometry by the specified distance\n",
    "    expanded_geometry = projected_geometry.buffer(expansion_distance_meters)\n",
    "    # Project the expanded geometry back to the original coordinate reference system\n",
    "    expanded_geometry = transform(partial(pyproj.transform, pyproj.Proj(proj='utm', zone=33, ellps='WGS84'), pyproj.Proj(init='epsg:4326')), expanded_geometry)\n",
    "    # Get the coordinates of the expanded bounding box\n",
    "    expanded_bbox = expanded_geometry.bounds\n",
    "    return expanded_bbox, expanded_geometry\n",
    "\n",
    "def create_bounding_box(geometry1, geometry2):\n",
    "\n",
    "    # Calculate the union of all polygons in each multipolygon\n",
    "    union_geometry1 = geometry1.convex_hull\n",
    "    union_geometry2 = geometry2.convex_hull\n",
    "    # Calculate the union of the convex hulls of the two multipolygons\n",
    "    union_geometry = union_geometry1.union(union_geometry2)\n",
    "    # Get the bounding box of the union geometry\n",
    "    bounding_box = union_geometry.bounds\n",
    "    return bounding_box\n",
    "\n",
    "def get_centrality(edges):\n",
    "\n",
    "    edge_list = []\n",
    "    for i,e in edges.iterrows():\n",
    "        edge_list.append(tuple([str(i[0]),str(i[1]),e['length'],e['LTS']]))\n",
    "        \n",
    "    g = gt.Graph()\n",
    "\n",
    "    elength = g.new_ep(\"float\")\n",
    "    elts = g.new_ep(\"int\")\n",
    "\n",
    "    g.add_edge_list(edge_list, hashed = True,eprops=[elength, elts])\n",
    "\n",
    "    vp, ep = gt.betweenness(g,weight=elength)\n",
    "\n",
    "    lts0 = []\n",
    "    lts1 = []\n",
    "    lts2 = []\n",
    "    lts3 = []\n",
    "    lts4 = []\n",
    "\n",
    "    for e in g.edges():\n",
    "        if elts[e] == 0:\n",
    "            lts0.append(ep[e])\n",
    "        elif elts[e] == 1:\n",
    "            lts1.append(ep[e])\n",
    "        elif elts[e] == 2:\n",
    "            lts2.append(ep[e])\n",
    "        elif elts[e] == 3:\n",
    "            lts3.append(ep[e])\n",
    "        elif elts[e] == 4:\n",
    "            lts4.append(ep[e])\n",
    "\n",
    "    return lts0,lts1,lts2,lts3,lts4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1025/1176549301.py:12: DtypeWarning: Columns (9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  lookup = pd.read_csv('data/PCD11_OA11_LSOA11_MSOA11_LAD11_EW_LU_aligned_v2.csv')\n"
     ]
    }
   ],
   "source": [
    "# Import MSOA lookup\n",
    "msoas = gpd.read_file('data/MSOA_EngWal_Dec_2011_Generalised_ClippedEW_0/Middle_Layer_Super_Output_Areas_December_2011_Generalised_Clipped_Boundaries_in_England_and_Wales.shp').to_crs(4326).set_index('msoa11cd')\n",
    "\n",
    "# Import MSOA 2011 OD data\n",
    "od_data = pd.read_parquet('data/od_2011.parquet')\n",
    "\n",
    "#Import LSOAs\n",
    "lsoas = gpd.read_file('data/LSOA_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V4_6029841263726194941.gpkg').to_crs(4326)\n",
    "lsoas = pd.concat([lsoas, lsoas.bounds], axis=1)\n",
    "\n",
    "#Import lsoa to msoa look up\n",
    "lookup = pd.read_csv('data/PCD11_OA11_LSOA11_MSOA11_LAD11_EW_LU_aligned_v2.csv')\n",
    "\n",
    "#York LSOAs\n",
    "york_model = pd.read_csv('data/LSOA_york_model.csv')\n",
    "york_lsoas = list(york_model['LSOA_code'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for lsoa_id in york_lsoas:\n",
    "\n",
    "lsoa_id = york_lsoas[0]\n",
    "\n",
    "#lsoa = lsoas[lsoas['LSOA11CD'] == test_lsoa['LSOA11CD'].values[0]]\n",
    "lsoa_lookup = lookup[lookup['LSOA11CD'] == lsoa_id][:1]\n",
    "lsoa = lsoas[lsoas['LSOA11CD'] == lsoa_lookup['LSOA11CD'].values[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/shapely/geometry/polygon.py:91: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead\n",
      "  return [float(c) for c in o]\n"
     ]
    }
   ],
   "source": [
    "lsoa_bbox = box(lsoa['maxy'],lsoa['miny'],lsoa['minx'],lsoa['maxx'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "intersecting_nodes = nodes[nodes.intersects(lsoa_bbox)].index\n",
    "G_bb = G.subgraph(intersecting_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/shapely/geometry/polygon.py:91: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead\n",
      "  return [float(c) for c in o]\n",
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
      "  in_crs_string = _prepare_from_proj_string(in_crs_string)\n",
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/shapely/ops.py:276: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n",
      "  shell = type(geom.exterior)(zip(*func(*zip(*geom.exterior.coords))))\n",
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
      "  in_crs_string = _prepare_from_proj_string(in_crs_string)\n",
      "/home/chris/miniconda3/envs/ate/lib/python3.10/site-packages/shapely/ops.py:276: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n",
      "  shell = type(geom.exterior)(zip(*func(*zip(*geom.exterior.coords))))\n"
     ]
    }
   ],
   "source": [
    "expanded_bbox, expanded_geometry = expand_bbox((lsoa['minx'], lsoa['miny'], lsoa['maxx'], lsoa['maxy']), expansion_distance_km=bbx_expansion)\n",
    "intersecting_nodes = nodes[nodes.intersects(box(expanded_bbox[0],expanded_bbox[1],expanded_bbox[2],expanded_bbox[3]))].index\n",
    "G_bb_exp = G.subgraph(intersecting_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "bike_ods = od_data[(od_data['geo_code1'] == lsoa_lookup['MSOA11CD'].values[0]) & (od_data['bicycle'] > 0)][['geo_code2','bicycle']].set_index('geo_code2')\n",
    "bike_ods['geometry'] = msoas['geometry']\n",
    "bike_ods = bike_ods.dropna()\n",
    "\n",
    "origin_geom = lsoa['geometry'].values[0]\n",
    "destination_geom = msoas.loc[bike_ods['bicycle'].idxmax()]['geometry']\n",
    "\n",
    "bounding_box_od = create_bounding_box(origin_geom, destination_geom)\n",
    "\n",
    "intersecting_nodes = nodes[nodes.intersects(box(bounding_box_od[0],bounding_box_od[1],bounding_box_od[2],bounding_box_od[3]))].index\n",
    "G_bb_od = G.subgraph(intersecting_nodes)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ate",
   "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.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}