import mermaid from ''; mermaid.initialize({ startOnLoad: true });
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()
Find the shortest way going through all streets from a set of streets? This problem is known as the Route inspection problem.
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 |
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%")
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%")
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));