01_make_network_graph.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "boston-metro\n"
     ]
    }
   ],
   "source": [
    "import copy\n",
    "import urllib.request\n",
    "import random\n",
    "import folium\n",
    "import xml.sax\n",
    "from math import radians, cos, sin, asin, sqrt\n",
    "from pathlib import Path\n",
    "import networkx as nx\n",
    "from pprint import pprint\n",
    "\n",
    "\n",
    "from vars import DATASET\n",
    "print(DATASET)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": [
     0,
     20,
     31
    ]
   },
   "outputs": [],
   "source": [
    "def haversine_distance(lon1, lat1, lon2, lat2, unit_m=True):\n",
    "    \"\"\"\n",
    "    Calculate the great circle distance between two points\n",
    "    on the earth (specified in decimal degrees)\n",
    "    default unit : km\n",
    "    \"\"\"\n",
    "    # convert decimal degrees to radians\n",
    "    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])\n",
    "\n",
    "    # haversine formula\n",
    "    dlon = lon2 - lon1\n",
    "    dlat = lat2 - lat1\n",
    "    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2\n",
    "    c = 2 * asin(sqrt(a))\n",
    "    r = 6371  # Radius of the Earth in kilometers. Use 3956 for miles\n",
    "    if unit_m:\n",
    "        r *= 1000\n",
    "    return c * r\n",
    "\n",
    "\n",
    "class Node(object):\n",
    "    def __init__(self, id, lon, lat):\n",
    "        self.id = id\n",
    "        self.lon = lon\n",
    "        self.lat = lat\n",
    "        self.tags = {}\n",
    "\n",
    "    def __str__(self):\n",
    "        return \"Node (id : %s) lon : %s, lat : %s \"%(self.id, self.lon, self.lat)\n",
    "\n",
    "\n",
    "class Way(object):\n",
    "    def __init__(self, id, osm):\n",
    "        self.osm = osm\n",
    "        self.id = id\n",
    "        self.nds = []\n",
    "        self.tags = {}\n",
    "\n",
    "    def split(self, dividers):\n",
    "        # slice the node-array using this nifty recursive function\n",
    "        def slice_array(ar, dividers):\n",
    "            for i in range(1,len(ar)-1):\n",
    "                if dividers[ar[i]]>1:\n",
    "                    left = ar[:i+1]\n",
    "                    right = ar[i:]\n",
    "\n",
    "                    rightsliced = slice_array(right, dividers)\n",
    "\n",
    "                    return [left]+rightsliced\n",
    "            return [ar]\n",
    "\n",
    "        slices = slice_array(self.nds, dividers)\n",
    "\n",
    "        # create a way object for each node-array slice\n",
    "        ret = []\n",
    "        i = 0\n",
    "        for slice in slices:\n",
    "            littleway = copy.copy(self)\n",
    "            littleway.id += \"-%d\" % i\n",
    "            littleway.nds = slice\n",
    "            ret.append(littleway)\n",
    "            i += 1\n",
    "\n",
    "        return ret\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "class OSM(object):\n",
    "    def __init__(self, osm_xml_data, is_xml_string=True):\n",
    "        \"\"\" File can be either a filename or stream/file object.\n",
    "\n",
    "        set `is_xml_string=False` if osm_xml_data is a filename or a file stream.\n",
    "        \"\"\"\n",
    "        nodes = {}\n",
    "        ways = {}\n",
    "\n",
    "        superself = self\n",
    "\n",
    "        class OSMHandler(xml.sax.ContentHandler):\n",
    "            @classmethod\n",
    "            def setDocumentLocator(self, loc):\n",
    "                pass\n",
    "\n",
    "            @classmethod\n",
    "            def startDocument(self):\n",
    "                pass\n",
    "\n",
    "            @classmethod\n",
    "            def endDocument(self):\n",
    "                pass\n",
    "\n",
    "            @classmethod\n",
    "            def startElement(self, name, attrs):\n",
    "                if name == 'node':\n",
    "                    self.currElem = Node(attrs['id'], float(attrs['lon']), float(attrs['lat']))\n",
    "                elif name == 'way':\n",
    "                    self.currElem = Way(attrs['id'], superself)\n",
    "                elif name == 'tag':\n",
    "                    self.currElem.tags[attrs['k']] = attrs['v']\n",
    "                elif name == 'nd':\n",
    "                    self.currElem.nds.append(attrs['ref'])\n",
    "\n",
    "            @classmethod\n",
    "            def endElement(self, name):\n",
    "                if name == 'node':\n",
    "                    nodes[self.currElem.id] = self.currElem\n",
    "                elif name == 'way':\n",
    "                    ways[self.currElem.id] = self.currElem\n",
    "\n",
    "            @classmethod\n",
    "            def characters(self, chars):\n",
    "                pass\n",
    "\n",
    "        if is_xml_string:\n",
    "            xml.sax.parseString(osm_xml_data, OSMHandler)\n",
    "        else:\n",
    "            with open(osm_xml_data, mode='r') as f:\n",
    "                xml.sax.parse(f, OSMHandler)\n",
    "\n",
    "        self.nodes = nodes\n",
    "        self.ways = ways\n",
    "\n",
    "        # count times each node is used\n",
    "        node_histogram = dict.fromkeys(self.nodes.keys(), 0)\n",
    "        for way in self.ways.values():\n",
    "            if len(way.nds) < 2:  # if a way has only one node, delete it out of the osm collection\n",
    "                del self.ways[way.id]\n",
    "            else:\n",
    "                for node in way.nds:\n",
    "                    node_histogram[node] += 1\n",
    "\n",
    "        # use that histogram to split all ways, replacing the member set of ways\n",
    "        new_ways = {}\n",
    "        for id, way in self.ways.items():\n",
    "            split_ways = way.split(node_histogram)\n",
    "            for split_way in split_ways:\n",
    "                new_ways[split_way.id] = split_way\n",
    "        self.ways = new_ways\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_osm(osm_xml_data, is_xml_string=True, only_roads=True):\n",
    "    \"\"\"Read graph in OSM format from file specified by name or by stream object.\n",
    "    Parameters\n",
    "    ----------\n",
    "    filename_or_stream : filename or stream object\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    G : Graph\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    >>> G=nx.read_osm(nx.download_osm(-122.33,47.60,-122.31,47.61))\n",
    "    >>> import matplotlib.pyplot as plt\n",
    "    >>> plt.plot([G.node[n]['lat']for n in G], [G.node[n]['lon'] for n in G], 'o', color='k')\n",
    "    >>> plt.show()\n",
    "    \"\"\"\n",
    "    osm = OSM(osm_xml_data, is_xml_string=is_xml_string)\n",
    "    G = nx.DiGraph()\n",
    "\n",
    "    ## Add ways\n",
    "    for w in osm.ways.values():\n",
    "        if only_roads and 'highway' not in w.tags:\n",
    "            continue\n",
    "          \n",
    "        if ('oneway' in w.tags):\n",
    "            if (w.tags['oneway'] == 'yes'):\n",
    "                # ONLY ONE DIRECTION\n",
    "                nx.add_path(G, w.nds, id=w.id, tags=w.tags)\n",
    "            else:\n",
    "                # BOTH DIRECTION\n",
    "                nx.add_path(G, w.nds, id=w.id, tags=w.tags)\n",
    "                nx.add_path(G, w.nds[::-1], id=w.id, tags=w.tags)        \n",
    "                \n",
    "        else:\n",
    "            # BOTH DIRECTION\n",
    "            nx.add_path(G, w.nds, id=w.id, tags=w.tags)\n",
    "            nx.add_path(G, w.nds[::-1], id=w.id, tags=w.tags)\n",
    "#             nx.add_path(G, w.nds, id=w.id)\n",
    "#             nx.add_path(G, w.nds[::-1], id=w.id)\n",
    "\n",
    "    # Complete the used nodes' information\n",
    "    coordinates_map = {}\n",
    "    for n_id in G.nodes():\n",
    "        n = osm.nodes[n_id]\n",
    "        G.nodes[n_id]['lat'] = n.lat\n",
    "        G.nodes[n_id]['lon'] = n.lon\n",
    "        G.nodes[n_id]['id'] = n.id\n",
    "        G.nodes[n_id]['tags'] = n.tags\n",
    "        coordinates_map[n_id] = (n.lon, n.lat)\n",
    "\n",
    "    # Estimate the length of each way\n",
    "    for u, v, d in G.edges(data=True):\n",
    "        # Give a realistic distance estimation (neither EPSG nor projection nor reference system are specified)\n",
    "        distance = haversine_distance(G.nodes[u]['lon'], G.nodes[u]['lat'], G.nodes[v]['lon'], G.nodes[v]['lat'], unit_m=True)\n",
    "\n",
    "        G.add_weighted_edges_from([(u, v, distance)], weight='havlen')\n",
    "\n",
    "#     G = nx.relabel_nodes(G, coordinates_map)\n",
    "    return G\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "918800"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = read_osm(\"./data/osm/\" + DATASET + \".osm\", is_xml_string=False)\n",
    "# make it undirected, because directed just isn't worth the trouble yet\n",
    "# g = g.to_undirected()\n",
    "len(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write to Disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = nx.read_gpickle(\"./data/01_\" + DATASET + \".gpickle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "918800"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chunks(lst, n):\n",
    "    \"\"\"Yield successive n-sized chunks from lst.\"\"\"\n",
    "    for i in range(0, len(lst), n):\n",
    "        yield lst[i:i + n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_chunks = list(chunks(list(g.nodes()), 300000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['./data/sylvia_exports/01_01_boston-metro_SYLVIA_NODES.csv',\n",
       " './data/sylvia_exports/01_02_boston-metro_SYLVIA_NODES.csv',\n",
       " './data/sylvia_exports/01_03_boston-metro_SYLVIA_NODES.csv',\n",
       " './data/sylvia_exports/01_04_boston-metro_SYLVIA_NODES.csv']"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filenames = ['01_01_', '01_02_', '01_03_', '01_04_']\n",
    "filenames = ['./data/sylvia_exports/' + name + DATASET + '_SYLVIA_NODES.csv' for name in filenames]\n",
    "filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# length of filenames must match length of chunks\n",
    "len(node_chunks) == len(filenames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, filename in enumerate(filenames):\n",
    "    f = open(filename, 'a')\n",
    "    chunk = node_chunks[idx]\n",
    "    for node in chunk:\n",
    "        data = g.nodes()[node]\n",
    "        nid, lat, lon = data['id'], data['lat'], data['lon']\n",
    "        f.write(str(nid) + ',' + str(lat) + ',' + str(lon) + '\\n')\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "expected an indented block (<ipython-input-44-74e4423515e2>, line 8)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-44-74e4423515e2>\"\u001b[0;36m, line \u001b[0;32m8\u001b[0m\n\u001b[0;31m    data = g.nodes()[node]\u001b[0m\n\u001b[0m       ^\u001b[0m\n\u001b[0;31mIndentationError\u001b[0m\u001b[0;31m:\u001b[0m expected an indented block\n"
     ]
    }
   ],
   "source": [
    "# Export for Sylvia\n",
    "## nodes\n",
    "# f = open(\"./data/sylvia_exports/01_\" + DATASET + \"_unpruned_nodes_for_sylvia.csv\", \"a\")\n",
    "# f.write(\"id,lat,lon\\n\")\n",
    "# for chunk in node_chunks:\n",
    "#     for node in chunk:\n",
    "#         data = g.nodes()[node]\n",
    "#         nid, lat, lon = data['id'], data['lat'], data['lon']\n",
    "#         f.write(str(nid) + ',' + str(lat) + ',' + str(lon) + '\\n')\n",
    "# f.close()\n",
    "\n",
    "## edges\n",
    "f = open(\"./data/sylvia_exports/01_\" + DATASET + \"_unpruned_edges_for_sylvia.csv\", \"a\")\n",
    "f.write(\"id1,id2\\n\")\n",
    "\n",
    "for edge in g.edges():\n",
    "    f.write(edge[0] + ',' + edge[1] + '\\n')\n",
    "f.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.write_gpickle(g, \"./data/01_\" + DATASET + \".gpickle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export to graphml. graphml export can't handle some datatypes (including dicts and None)\n",
    "# so delete them\n",
    "h = g.copy()\n",
    "\n",
    "for n in h.nodes():\n",
    "    h.nodes[n].pop('tags', None)\n",
    "for e in g.edges():\n",
    "    h.edges[e].pop('tags', None)\n",
    "nx.write_graphml(h, \"./data/01_\" + DATASET + \".graphml\")\n",
    "\n",
    "# clear up memory\n",
    "h = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cruft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dictionary of Tags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "tags = {}\n",
    "osm = OSM('./' + DATASET + '.osm', is_xml_string=False)\n",
    "\n",
    "for w in osm.ways.values():\n",
    "    for t in w.tags:\n",
    "        if tags.get(t):\n",
    "            tags[t] = tags[t] + 1\n",
    "        else:\n",
    "            tags[t] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'attribution': 146,\n",
       " 'condition': 141,\n",
       " 'highway': 205,\n",
       " 'lanes': 142,\n",
       " 'massgis:way_id': 142,\n",
       " 'maxspeed': 107,\n",
       " 'name': 169,\n",
       " 'oneway': 108,\n",
       " 'source': 154,\n",
       " 'surface': 126,\n",
       " 'width': 131,\n",
       " 'cycleway:right': 15,\n",
       " 'massgis:ref': 16,\n",
       " 'parking:lane:left': 10,\n",
       " 'lit': 3,\n",
       " 'foot': 12,\n",
       " 'access': 13,\n",
       " 'building': 643,\n",
       " 'addr:housenumber': 508,\n",
       " 'addr:street': 507,\n",
       " 'amenity': 17,\n",
       " 'brand': 5,\n",
       " 'brand:wikidata': 4,\n",
       " 'brand:wikipedia': 4,\n",
       " 'dispensing': 3,\n",
       " 'drive_through': 3,\n",
       " 'healthcare': 3,\n",
       " 'opening_hours': 6,\n",
       " 'payment:cash': 3,\n",
       " 'payment:visa': 3,\n",
       " 'phone': 7,\n",
       " 'website': 10,\n",
       " 'wheelchair': 3,\n",
       " 'addr:city': 10,\n",
       " 'addr:postcode': 11,\n",
       " 'building:levels': 8,\n",
       " 'operator': 3,\n",
       " 'power': 3,\n",
       " 'ref': 2,\n",
       " 'substation': 1,\n",
       " 'voltage': 3,\n",
       " 'addr:state': 5,\n",
       " 'shop': 2,\n",
       " 'denomination': 3,\n",
       " 'religion': 3,\n",
       " 'area': 8,\n",
       " 'created_by': 6,\n",
       " 'leisure': 20,\n",
       " 'massgis:ARTICLE97': 8,\n",
       " 'massgis:ASSESS_ACR': 8,\n",
       " 'massgis:ATT_DATE': 8,\n",
       " 'massgis:DCAM_ID': 8,\n",
       " 'massgis:DEED_ACRES': 8,\n",
       " 'massgis:EOEAINVOLV': 8,\n",
       " 'massgis:FEE_OWNER': 8,\n",
       " 'massgis:FEESYM': 8,\n",
       " 'massgis:FY_FUNDING': 8,\n",
       " 'massgis:LEV_PROT': 8,\n",
       " 'massgis:OS_DEED_BO': 8,\n",
       " 'massgis:OS_DEED_PA': 8,\n",
       " 'massgis:OS_ID': 8,\n",
       " 'massgis:OWNER_ABRV': 8,\n",
       " 'massgis:OWNER_TYPE': 8,\n",
       " 'massgis:POLY_ID': 8,\n",
       " 'massgis:PRIM_PURP': 8,\n",
       " 'massgis:PUB_ACCESS': 8,\n",
       " 'massgis:SITE_NAME': 8,\n",
       " 'massgis:TOWN_ID': 8,\n",
       " 'owner': 8,\n",
       " 'ownership': 8,\n",
       " 'protected': 8,\n",
       " 'wikidata': 10,\n",
       " 'cuisine': 1,\n",
       " 'outdoor_seating': 1,\n",
       " 'takeaway': 1,\n",
       " 'landuse': 11,\n",
       " 'massgis:SOURCE_MAP': 1,\n",
       " 'building:units': 1,\n",
       " 'ele': 2,\n",
       " 'email': 2,\n",
       " 'gnis:county_id': 2,\n",
       " 'gnis:created': 2,\n",
       " 'gnis:feature_id': 2,\n",
       " 'gnis:state_id': 2,\n",
       " 'addr:unit': 1,\n",
       " 'barrier': 11,\n",
       " 'sport': 9,\n",
       " 'fee': 4,\n",
       " 'parking': 4,\n",
       " 'park_ride': 3,\n",
       " 'supervised': 3,\n",
       " 'service': 9,\n",
       " 'sidewalk': 27,\n",
       " 'bicycle': 6,\n",
       " 'hazmat': 6,\n",
       " 'cables': 2,\n",
       " 'circuits': 2,\n",
       " 'frequency': 2,\n",
       " 'location': 2,\n",
       " 'capacity': 1,\n",
       " 'covered': 1,\n",
       " 'residential': 1,\n",
       " 'playground': 2,\n",
       " 'footway': 10,\n",
       " 'crossing': 6,\n",
       " 'brand:website': 1,\n",
       " 'contact:facebook': 1,\n",
       " 'contact:instagram': 1,\n",
       " 'contact:twitter': 1,\n",
       " 'network': 1,\n",
       " 'old_name': 1,\n",
       " 'operator:type': 1,\n",
       " 'payment:credit_cards': 1,\n",
       " 'payment:debit_cards': 1,\n",
       " 'payment:mobile_app': 1,\n",
       " 'start_date': 1,\n",
       " 'type': 1,\n",
       " 'wikipedia': 1,\n",
       " 'public_transport:version': 1,\n",
       " 'route': 1,\n",
       " 'day_off': 1,\n",
       " 'day_on': 1,\n",
       " 'hour_off': 1,\n",
       " 'hour_on': 1,\n",
       " 'restriction': 1,\n",
       " 'fixme': 1,\n",
       " 'restriction:hgv': 1}"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bikepoints",
   "language": "python",
   "name": "bikepoints"
  },
  "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.7.0"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}