%matplotlib inline
import geopandas
import matplotlib.pyplot as plt
from rivus.utils.pandashp import match_vertices_and_edges
from rivus.utils import shapelytools
from shapely.geometry import (box, LineString, MultiLineString, MultiPoint, Point, Polygon)
#import skeletrontools
fpath = "./data/haag15/{}".format
streets_filename = fpath('streets.shp')
edge_filename = fpath('edge.shp')
vertex_filename = fpath('vertex.shp')
streets = geopandas.read_file(streets_filename)
streets = streets.to_crs(epsg=32632) # EPSG:32632 == UTM Zone 32N (Germany!)
streets.head(10)
geometry | id | name | osm_id | other_tags | type | |
---|---|---|---|---|---|---|
0 | LINESTRING (734017.7521000658 5341822.15924835... | 1 | None | 4050378 | "maxspeed"=>"100","ref"=>"B 15","surface"=>"as... | primary |
1 | LINESTRING (736871.4876152406 5339170.68830887... | 1 | Hauptstraße | 4050621 | "maxspeed"=>"50" | secondary |
2 | LINESTRING (735627.0591133552 5338650.80274698... | 1 | Tannenstraße | 4050749 | None | residential |
3 | LINESTRING (733629.1316536779 5341877.42621907... | 1 | None | 9660295 | None | residential |
4 | LINESTRING (735678.8561374607 5339439.84541021... | 1 | None | 9931785 | None | service |
5 | LINESTRING (735532.1246327882 5339407.67760010... | 1 | Gerberstraße | 9931784 | None | residential |
6 | LINESTRING (734226.9737598393 5341329.16469117... | 1 | None | 11538813 | "surface"=>"unpaved","tracktype"=>"grade3" | track |
7 | LINESTRING (736712.5506254345 5340021.60657790... | 1 | None | 12956380 | None | unclassified |
8 | LINESTRING (736725.7023614591 5339111.27784710... | 1 | None | 12957323 | None | footway |
9 | LINESTRING (736659.5552210254 5339085.82769116... | 1 | Marktplatz | 12957325 | None | residential |
# filter away roads by type
road_types = ['motorway', 'motorway_link', 'primary', 'primary_link',
'secondary', 'secondary_link', 'tertiary', 'tertiary_link',
'residential', 'living_street', 'service', 'unclassified']
streets = streets[streets['type'].isin(road_types)]
streets.head(10)
geometry | id | name | osm_id | other_tags | type | |
---|---|---|---|---|---|---|
0 | LINESTRING (734017.7521000658 5341822.15924835... | 1 | None | 4050378 | "maxspeed"=>"100","ref"=>"B 15","surface"=>"as... | primary |
1 | LINESTRING (736871.4876152406 5339170.68830887... | 1 | Hauptstraße | 4050621 | "maxspeed"=>"50" | secondary |
2 | LINESTRING (735627.0591133552 5338650.80274698... | 1 | Tannenstraße | 4050749 | None | residential |
3 | LINESTRING (733629.1316536779 5341877.42621907... | 1 | None | 9660295 | None | residential |
4 | LINESTRING (735678.8561374607 5339439.84541021... | 1 | None | 9931785 | None | service |
5 | LINESTRING (735532.1246327882 5339407.67760010... | 1 | Gerberstraße | 9931784 | None | residential |
7 | LINESTRING (736712.5506254345 5340021.60657790... | 1 | None | 12956380 | None | unclassified |
9 | LINESTRING (736659.5552210254 5339085.82769116... | 1 | Marktplatz | 12957325 | None | residential |
10 | LINESTRING (736700.4427397507 5339080.82016933... | 1 | Marktplatz | 12957327 | None | residential |
11 | LINESTRING (736722.7920502152 5339158.32395773... | 1 | Zeno-Kern-Straße | 12957778 | None | residential |
#skeleton = skeletrontools.skeletonize(streets,
# buffer_length=30,
# dissolve_length=15,
# simplify_length=30)
Create simple MultiLineString while I cannot install Skeletron
streets['geometry'][0]
ls1 = LineString([streets['geometry'][0].coords[0], streets['geometry'][0].coords[-1]])
ls2 = LineString([streets['geometry'][1].coords[0], streets['geometry'][1].coords[-1]])
MultiLineString([ls1, ls2])
streets['geometry'][0].bounds[2:]
(734570.5783283255, 5341822.159248358)
streets['geometry'][0].coords[-1]
(734570.5783283255, 5341197.624399373)
skeleton = MultiLineString([(line.coords[0], line.coords[-1]) for line in streets['geometry'] if isinstance(line, LineString)])
#skeleton = shapelytools.one_linestring_per_intersection(skeleton)
skeleton
This is a quite a strong distortion in production either skeletron should be resurrected, or an alternative should be found.
skeleton = shapelytools.snappy_endings(skeleton, max_distance=100)
skeleton[:3]
[<shapely.geometry.linestring.LineString at 0x11fc5a90>, <shapely.geometry.linestring.LineString at 0x11fae2b0>, <shapely.geometry.linestring.LineString at 0x11fe0400>]
skeleton = shapelytools.one_linestring_per_intersection(skeleton)
skeleton
skeleton = shapelytools.prune_short_lines(skeleton, min_length=55)
geopandas.GeoDataFrame(geometry=skeleton, crs=streets.crs).plot()
<matplotlib.axes._subplots.AxesSubplot at 0x11fd35c0>
# convert back to GeoPandas
edge = geopandas.GeoDataFrame(geometry=skeleton, crs=streets.crs)
edge = edge.to_crs(epsg=4326) # World Geodetic System (WGS84)
edge['Edge'] = edge.index
# derive vertex points
vertices = shapelytools.endpoints_from_lines(edge.geometry)
vertex = geopandas.GeoDataFrame(geometry=vertices, crs=edge.crs)
vertex['Vertex'] = vertex.index
vertex.head(10)
geometry | Vertex | |
---|---|---|
0 | POINT (12.18014738230067 48.1652135423499) | 0 |
1 | POINT (12.1883682 48.1663995) | 1 |
2 | POINT (12.19606713957078 48.15712345702804) | 2 |
3 | POINT (12.18007501864393 48.16522824633694) | 3 |
4 | POINT (12.1756668 48.1557562) | 4 |
5 | POINT (12.198809 48.158848) | 5 |
6 | POINT (12.15573914996715 48.18062360000686) | 6 |
7 | POINT (12.19349565563646 48.16369785138707) | 7 |
8 | POINT (12.18566450292734 48.15921207188419) | 8 |
9 | POINT (12.187404 48.15355780000001) | 9 |
match_vertices_and_edges(vertex, edge)
edge.head(10)
geometry | Edge | Vertex1 | Vertex2 | |
---|---|---|---|---|
0 | LINESTRING (12.13807510000005 48.1870267, 12.1... | 0 | 62 | 121 |
1 | LINESTRING (12.1394549 48.1855127, 12.14056790... | 1 | 62 | 156 |
2 | LINESTRING (12.14056790847939 48.1869023644199... | 2 | 19 | 62 |
3 | LINESTRING (12.1435042 48.1871218, 12.1449267 ... | 3 | 19 | 40 |
4 | LINESTRING (12.14527093227269 48.1866579163458... | 4 | 19 | 61 |
5 | LINESTRING (12.14527093227269 48.1866579163458... | 5 | 19 | 63 |
6 | LINESTRING (12.1486945 48.1864823, 12.14799653... | 6 | 63 | 159 |
7 | LINESTRING (12.1486945 48.1864823, 12.15573914... | 7 | 6 | 63 |
8 | LINESTRING (12.15573914996715 48.1806236000068... | 8 | 6 | 132 |
9 | LINESTRING (12.15573914996715 48.1806236000068... | 9 | 6 | 163 |
# drop loops
edge = edge[~(edge['Vertex1'] == edge['Vertex2'])]
edge.head(10)
geometry | Edge | Vertex1 | Vertex2 | |
---|---|---|---|---|
0 | LINESTRING (12.13807510000005 48.1870267, 12.1... | 0 | 62 | 121 |
1 | LINESTRING (12.1394549 48.1855127, 12.14056790... | 1 | 62 | 156 |
2 | LINESTRING (12.14056790847939 48.1869023644199... | 2 | 19 | 62 |
3 | LINESTRING (12.1435042 48.1871218, 12.1449267 ... | 3 | 19 | 40 |
4 | LINESTRING (12.14527093227269 48.1866579163458... | 4 | 19 | 61 |
5 | LINESTRING (12.14527093227269 48.1866579163458... | 5 | 19 | 63 |
6 | LINESTRING (12.1486945 48.1864823, 12.14799653... | 6 | 63 | 159 |
7 | LINESTRING (12.1486945 48.1864823, 12.15573914... | 7 | 6 | 63 |
8 | LINESTRING (12.15573914996715 48.1806236000068... | 8 | 6 | 132 |
9 | LINESTRING (12.15573914996715 48.1806236000068... | 9 | 6 | 163 |
# if there are >1 edges that connect the same vertex pair, drop one of them
edge = edge.drop_duplicates(subset=['Vertex1', 'Vertex2'])
edge.head(20)
geometry | Edge | Vertex1 | Vertex2 | |
---|---|---|---|---|
0 | LINESTRING (12.13807510000005 48.1870267, 12.1... | 0 | 62 | 121 |
1 | LINESTRING (12.1394549 48.1855127, 12.14056790... | 1 | 62 | 156 |
2 | LINESTRING (12.14056790847939 48.1869023644199... | 2 | 19 | 62 |
3 | LINESTRING (12.1435042 48.1871218, 12.1449267 ... | 3 | 19 | 40 |
4 | LINESTRING (12.14527093227269 48.1866579163458... | 4 | 19 | 61 |
5 | LINESTRING (12.14527093227269 48.1866579163458... | 5 | 19 | 63 |
6 | LINESTRING (12.1486945 48.1864823, 12.14799653... | 6 | 63 | 159 |
7 | LINESTRING (12.1486945 48.1864823, 12.15573914... | 7 | 6 | 63 |
8 | LINESTRING (12.15573914996715 48.1806236000068... | 8 | 6 | 132 |
9 | LINESTRING (12.15573914996715 48.1806236000068... | 9 | 6 | 163 |
10 | LINESTRING (12.1572016 48.18036679999999, 12.1... | 10 | 125 | 132 |
11 | LINESTRING (12.1604319 48.1661044, 12.15676416... | 11 | 29 | 160 |
12 | LINESTRING (12.15777632897053 48.1808760262701... | 12 | 18 | 132 |
13 | LINESTRING (12.16223287484881 48.1705480751595... | 13 | 59 | 120 |
14 | LINESTRING (12.1612266 48.17295960000001, 12.1... | 14 | 60 | 163 |
15 | LINESTRING (12.1604319 48.1661044, 12.1623088 ... | 15 | 55 | 160 |
16 | LINESTRING (12.1628662 48.1655761, 12.1604319 ... | 16 | 55 | 161 |
17 | LINESTRING (12.1612266 48.17295960000001, 12.1... | 17 | 59 | 163 |
18 | LINESTRING (12.16223287484881 48.1705480751595... | 18 | 30 | 59 |
19 | LINESTRING (12.1624808 48.16757969999999, 12.1... | 19 | 83 | 161 |
... wirte to file
fig, axes = plt.subplots(1, 2, figsize=(12,8))
streets.plot(ax=axes[0], color='black')
edge.to_crs(crs=streets.crs).plot(ax=axes[1], color='blue')
<matplotlib.axes._subplots.AxesSubplot at 0x12519278>