Table Of Contents
Table Of Contents

Walk through all streets in a city

Links: notebook, city_tour_data_preparation2html.html, PDF, python, city_tour_data_preparation.slides.html, city_tour_data_preparation.slides2p.html, GitHub

Preparation of the examples for the challenge: find the shortest path through a set of streets.

import matplotlib.pyplot as plt
%matplotlib inline
from jyquickhelper import add_notebook_menu
add_notebook_menu()

Problem description

Find the shortest way going through all streets from a set of streets? This problem is known as the Route inspection problem.

Data

Seattle streets from data.seattle.gov

Read the data

import shapefile, os
if os.path.exists("Street_Network_Database/WGS84/Street_Network_Database.shp"):
    rshp = shapefile.Reader("Street_Network_Database/WGS84/Street_Network_Database.shp")
    shapes = rshp.shapes()
    records = rshp.records()
else:
    from pyensae.datasource import download_data
    files = download_data("WGS84_seattle_street.zip")
    rshp = shapefile.Reader("Street_Network_Database.shp")
    shapes = rshp.shapes()
    records = rshp.records()
shapes[0].__dict__
{'shapeType': 3,
 'points': [(-122.34651954599997, 47.46947199700003),
  (-122.34721334599999, 47.46946425100003)],
 'parts': [0],
 'bbox': [-122.34721334599999, 47.46946425100003, -122.346519546, 47.469471997000035]}
{k[0]:v for k,v in zip(rshp.fields[1:], records[0])}
{'F_INTR_ID': 21642,
 'T_INTR_ID': 21641,
 'SND_ID': 37898,
 'SND_FEACOD': 1,
 'CITYCODE': 0,
 'STNAME_ID': 2569,
 'ST_CODE': 0,
 'ARTERIAL_C': 0,
 'SEGMENT_TY': 1,
 'AGENCY_COD': 1,
 'ACCESS_COD': 1,
 'DIVIDED_CO': 1,
 'STRUCTURE_': 1,
 'LEGALLOC_C': 1,
 'VEHICLE_US': 1,
 'GIS_SEG_LE': 171.624048,
 'L_ADRS_FRO': 977,
 'L_ADRS_TO': 999,
 'R_ADRS_FRO': 976,
 'R_ADRS_TO': 998,
 'ORD_PRE_DI': 'SW',
 'ORD_STREET': '149TH',
 'ORD_STRE_1': 'ST',
 'ORD_SUF_DI': '',
 'ORD_STNAME': 'SW 149TH ST',
 'L_CITY': 'BURIEN',
 'L_STATE': 'WA',
 'L_ZIP': '98166',
 'R_CITY': 'BURIEN',
 'R_STATE': 'WA',
 'R_ZIP': '98166',
 'SNDSEG_UPD': datetime.date(2004, 5, 19),
 'COMPKEY': 0,
 'COMPTYPE': 0,
 'UNITID': '0',
 'UNITID2': '0',
 'SHAPE_Leng': 0.000693843239173}
from ensae_projects.datainc.data_geo_streets import get_fields_description
get_fields_description()
SND_FEACODE SEGMENT_TYPE AGENCY_CODE ACCESS_CODE DIVIDED_CODE STRUCTURE_TYPE ARTERIAL_CODE
0 0 Alias Street 14 (stub) 9 (stub) 6 (stub) 7 (stub) 4(stub) 0 (not an arterial)
1 1 Local Street 1 (street) any < 8 and = 10 any < 4 and = 7 any < 7 values less than 4 0 (not an arterial)
2 NaN 2 (street ramp) NaN NaN NaN NaN NaN
3 NaN 5 (alley) NaN NaN NaN NaN NaN
4 NaN 10 (dock) NaN NaN NaN NaN NaN
5 5 Major Street 1 (street) any < 8 and = 10 any < 4 and = 7 any < 7 values less than 4 1 (is an arterial)
6 NaN 2 (street ramp) NaN NaN NaN NaN NaN
7 9 Divided Highway 3 (limited access) 1 (public) 1 (open) 6 (LimitedAcccess) values less than 4 4 (highway)
8 NaN 4 (limited access ramp) 4 (WSDOT) NaN NaN NaN NaN
9 13 Interstate Freeway 3 (limited access) 4 (WSDOT) 1 (open) 6 (LimitedAcccess) values less than 4 5 (freeway)
10 NaN 4 (limited access ramp) NaN NaN NaN NaN NaN
11 77 Stairs, Trail, Walkway 6 (stairs) any <> 9 any <> 6 any any 0 (not an arterial)
12 NaN 7 (walkway) NaN NaN NaN NaN NaN
13 NaN 8 (trail) NaN NaN NaN NaN NaN
14 85 Railroad 9 (railroad) 8 (railroad) 3 (restricted) 1 (undivided) values less than 4 0 (not an arterial)
15 NaN 11 (light rail) 1 (public) NaN NaN NaN NaN
16 NaN 12 (monorail) 7 (sound transit) NaN NaN NaN NaN
17 NaN 13 (trolley) NaN NaN NaN NaN NaN

Display the streets

streets5 = list(zip(records[:5], shapes[:5]))
streets5[2][1].points
[(-122.31554790099995, 47.511287922000065),
 (-122.31553241799998, 47.51120351700007),
 (-122.31552978999997, 47.51118938700006),
 (-122.31546052299996, 47.51092530900007),
 (-122.31537415499997, 47.510596031000034),
 (-122.31534125099995, 47.51046084500007),
 (-122.31532328399999, 47.51032437400005)]
import folium
from random import randint
from pyensae.notebookhelper import folium_html_map

c = streets5[0][1]
map_osm = folium.Map(location=[c.bbox[1], c.bbox[0]], zoom_start=9)
for rec, shape in streets5:
    d = {k[0]: v for k,v in zip(rshp.fields[1:], rec)}
    map_osm.add_child(folium.Marker([shape.points[0][1], shape.points[0][0]], popup=str(d["ORD_STNAME"])))
    map_osm.add_child(folium.PolyLine(locations=[[_[1], _[0]] for _ in shape.points], weight=10))
folium_html_map(map_osm, width="60%")

Find connected streets

street0 = streets5[0][1].points
street0
[(-122.34651954599997, 47.46947199700003),
 (-122.34721334599999, 47.46946425100003)]
def connect_streets(st1, st2):
    a1, b1 = st1[0], st1[-1]
    a2, b2 = st2[0], st2[-1]
    connect = []
    if a1 == a2:
        connect.append((0, 0))
    if a1 == b2:
        connect.append((0, 1))
    if b1 == a2:
        connect.append((1, 0))
    if b1 == b2:
        connect.append((1, 1))
    return tuple(connect) if connect else None

neighbours = []
for i, street in enumerate(shapes):
    points = street.points
    con = connect_streets(street0, points)
    if con:
        neighbours.append(i)

neighbours
[0, 107, 1670, 9989, 11274, 12783]
import folium
from pyensae.notebookhelper import folium_html_map
c = shapes[neighbours[0]]
map_osm = folium.Map(location=[c.bbox[1], c.bbox[0]], zoom_start=15)
points = set()
for index in neighbours:
    rec, shape = records[index], shapes[index]
    corners = [(_[1], _[0]) for _ in shape.points]
    map_osm.add_child(folium.PolyLine(locations=corners, weight=10))
    points |= set([corners[0], corners[-1]])
for x, y in points:
    map_osm.add_child(folium.Marker((x, y), popup=str(index)))
folium_html_map(map_osm, width="50%")
c = shapes[neighbours[0]]
map_osm = folium.Map(location=[c.bbox[1], c.bbox[0]], zoom_start=15)
points = set()
for index in neighbours:
    rec, shape = records[index], shapes[index]
    corners = [(_[1], _[0]) for _ in shape.points]
    map_osm.add_child(folium.PolyLine(locations=corners, weight=10))
    points |= set([corners[0], corners[-1]])
for x, y in points:
    map_osm.add_child(folium.CircleMarker((x, y), popup=str(index), radius=8, fill_color="yellow"))
folium_html_map(map_osm, width="50%")

Extraction of all streets in a short perimeter

from shapely.geometry import Point, LineString

def enumerate_close(x, y, shapes, th=None):
    p = Point(x,y)
    for i, shape in enumerate(shapes):
        obj = LineString(shape.points)
        d = p.distance(obj)
        if th is None or d <= th:
            yield d, i

x, y = shapes[0].points[0]
closes = list(enumerate_close(x, y, shapes))
closes.sort()
closes[:10]
[(0.0, 0),
 (0.0, 1670),
 (0.0, 9989),
 (0.0, 12783),
 (0.0006938432391730961, 107),
 (0.0006938432391730961, 11274),
 (0.0009050591972649863, 9118),
 (0.0009122767287444535, 6488),
 (0.001006818609911273, 1101),
 (0.001006818609911273, 2808)]
import folium
from ensae_projects.datainc.data_geo_streets import folium_html_street_map
folium_html_street_map([_[1] for _ in closes[:20]], shapes, html_width="50%", zoom_start=15)
def complete_subset_streets(subset, shapes):
    extension = []
    for i, shape in enumerate(shapes):
        add = []
        for s in subset:
            to = shapes[s]
            if s != i:
                con = connect_streets(shapes[s].points, shapes[i].points)
                if con is not None:
                    add.extend([_[1] for _ in con])
        if len(set(add)) == 2:
            extension.append(i)
    return extension

subset = [index for dist, index in closes[:20]]
newset = set(subset + complete_subset_streets(subset, shapes))

print(list(sorted(newset)))
folium_html_street_map(newset, shapes, html_width="50%", zoom_start=15)
[0, 107, 1003, 1101, 1670, 2418, 2803, 2808, 3353, 4553, 4994, 6265, 6488, 6712, 8378, 9118, 9989, 11274, 11394, 12783, 15023, 17680, 29114, 30370]
from ensae_projects.datainc.data_geo_streets import build_streets_vertices
vertices, edges = build_streets_vertices(newset, shapes)
vertices[:3], edges[:3]
([(-122.34991548199997, 47.46763155800005),
  (-122.34991155699998, 47.468532819000075),
  (-122.349907514, 47.469446668000046)],
 [(10, 7), (5, 4), (4, 0)])
from ensae_projects.datainc.data_geo_streets import plot_streets_network
plot_streets_network(newset, edges, vertices, shapes, figsize=(10,10));
../_images/city_tour_data_preparation_24_0.png