.. _citytourdatapreparationrst: ================================== Walk through all streets in a city ================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`PDF `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/challenges/city_tour/city_tour_data_preparation.ipynb|*` Preparation of the examples for the challenge: find the shortest path through a set of streets. .. code:: ipython3 import matplotlib.pyplot as plt %matplotlib inline .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: 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 ~~~~~~~~~~~~~ .. code:: ipython3 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() .. code:: ipython3 shapes[0].__dict__ .. parsed-literal:: {'shapeType': 3, 'points': [(-122.34651954599997, 47.46947199700003), (-122.34721334599999, 47.46946425100003)], 'parts': [0], 'bbox': [-122.34721334599999, 47.46946425100003, -122.346519546, 47.469471997000035]} .. code:: ipython3 {k[0]:v for k,v in zip(rshp.fields[1:], records[0])} .. parsed-literal:: {'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} .. code:: ipython3 from ensae_projects.datainc.data_geo_streets import get_fields_description get_fields_description() .. raw:: html
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 ~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 streets5 = list(zip(records[:5], shapes[:5])) streets5[2][1].points .. parsed-literal:: [(-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)] .. code:: ipython3 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%") .. raw:: html
Find connected streets ---------------------- .. code:: ipython3 street0 = streets5[0][1].points street0 .. parsed-literal:: [(-122.34651954599997, 47.46947199700003), (-122.34721334599999, 47.46946425100003)] .. code:: ipython3 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 .. parsed-literal:: [0, 107, 1670, 9989, 11274, 12783] .. code:: ipython3 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%") .. raw:: html
.. code:: ipython3 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%") .. raw:: html
Extraction of all streets in a short perimeter ---------------------------------------------- .. code:: ipython3 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] .. parsed-literal:: [(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)] .. code:: ipython3 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) .. raw:: html
.. code:: ipython3 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) .. parsed-literal:: [0, 107, 1003, 1101, 1670, 2418, 2803, 2808, 3353, 4553, 4994, 6265, 6488, 6712, 8378, 9118, 9989, 11274, 11394, 12783, 15023, 17680, 29114, 30370] .. raw:: html
.. code:: ipython3 from ensae_projects.datainc.data_geo_streets import build_streets_vertices vertices, edges = build_streets_vertices(newset, shapes) vertices[:3], edges[:3] .. parsed-literal:: ([(-122.34991548199997, 47.46763155800005), (-122.34991155699998, 47.468532819000075), (-122.349907514, 47.469446668000046)], [(10, 7), (5, 4), (4, 0)]) .. code:: ipython3 from ensae_projects.datainc.data_geo_streets import plot_streets_network plot_streets_network(newset, edges, vertices, shapes, figsize=(10,10)); .. image:: city_tour_data_preparation_24_0.png