Coverage for src/ensae_projects/datainc/data_geo_streets.py : 88%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2"""
3@file
4@brief Data related to a challenge, streets in Seattle
5"""
6import os
7import pandas
8from pyensae.datasource import download_data
9from pyensae.notebookhelper import folium_html_map
12def get_fields_description():
13 """
14 Retrieves a :epkg:`dataframe` which describes
15 the meaning of the metadata.
17 @return dataframe
18 """
19 from .seattle_streets import __file__ as local_dir
20 this = os.path.join(os.path.dirname(local_dir), "street_seattle.desc.xlsx")
21 return pandas.read_excel(this)
24def get_seattle_streets(filename=None, folder="."):
25 """
26 Retrieves processed data from
27 `Seattle Streets <https://data.seattle.gov/dataset/Street-Network-Database/
28 afip-2mzr/data)>`_.
30 @param filename local filename
31 @param folder temporary folder where to download files
32 @return shapes, records
34 The function returns a filename.
35 """
36 if filename is None:
37 names = download_data("WGS84_seattle_street.zip", whereTo=folder)
38 shp = [n for n in names if n.endswith('.shp')]
39 if len(shp) != 1:
40 from pyquickhelper.loghelper import BufferedPrint
41 buf = BufferedPrint()
42 names = download_data("WGS84_seattle_street.zip",
43 whereTo=folder, fLOG=buf.fprint)
44 raise FileNotFoundError(
45 "Unable to download data 'WGS84_seattle_street.zip' to '{0}', log={1}\nnames={2}.".format(
46 filename, str(buf), "\n".join(names)))
47 filename = shp[0]
48 elif not os.path.exists(filename):
49 raise FileNotFoundError(filename)
50 return filename
53def shapely_records(filename, **kwargs):
54 """
55 Uses `pyshp <https://pypi.python.org/pypi/pyshp/>`_ to return
56 shapes and records from shapefiles.
58 @param filename filename
59 @param kwargs addition parameter for the shapefile reader,
60 useful options car *encoding* and *encodingErrors*
61 @return shapes, records, fields
63 .. faqref::
64 :title: Fields in shapefiles
66 The following codes is usually used to retrieve
67 shapefiles:
69 ::
71 rshp = shapefile.Reader(filename)
72 shapes = rshp.shapes()
73 records = rshp.records()
75 *records* contains metadata about each shape.
76 Fields and values are not stored in a dictionary.
77 Here is a snippet of code to do so:
79 ::
81 {k[0]: v for k, v in zip(rshp.fields[1:], records[0])}
83 Here is an example of the results:
85 ::
87 {'ORD_PRE_DI': 'SW',
88 'ORD_STNAME': 'SW 149TH ST',
89 'ORD_STREET': '149TH',
90 'ORD_STRE_1': 'ST',
91 'ORD_SUF_DI': b' ',
92 'R_ADRS_FRO': 976,
93 ...
95 """
96 import shapefile
97 rshp = shapefile.Reader(filename, **kwargs)
98 shapes = rshp.shapes()
99 records = rshp.records()
100 try:
101 rshp.shp.close()
102 except IOError:
103 pass
104 try:
105 rshp.shx.close()
106 except IOError:
107 pass
108 try:
109 rshp.dbf.close()
110 except IOError:
111 pass
112 return shapes, records, rshp.fields
115def folium_html_street_map(subset, shapes, html_width=None, html_height=None,
116 color_vertices=None, **kwargs):
117 """
118 Returns a :epkg:`folium` map which represents the streets.
120 @param subset list of streets index
121 @param shapes list of shapes
122 @param kwargs extra parameters for `Map <https://github.com/python-visualization/folium/blob/master/folium/folium.py#L19>`_
123 @param html_width sent to function
124 `folium_html_map <http://www.xavierdupre.fr/app/pyensae/helpsphinx/pyensae/notebookhelper/folium_helper.html
125 #pyensae.notebookhelper.folium_helper.folium_html_map>`_
126 @param html_height sent to function
127 `folium_html_map <http://www.xavierdupre.fr/app/pyensae/helpsphinx/pyensae/notebookhelper/folium_helper.html
128 #pyensae.notebookhelper.folium_helper.folium_html_map>`_
129 @param color_vertices see below
130 @return see function
131 `folium_html_map <http://www.xavierdupre.fr/app/pyensae/helpsphinx/pyensae/notebookhelper/folium_helper.html
132 #pyensae.notebookhelper.folium_helper.folium_html_map>`_
134 if *color_vertices* is equal to `'odd'`, the function computes the degree
135 of each vertex and choose a different color for odd (yellow)
136 and even degrees (black).
137 """
138 if color_vertices == "odd":
139 count = {}
140 for edge in subset:
141 shape = shapes[edge]
142 a = (shape.points[0][0], shape.points[0][1])
143 b = (shape.points[-1][0], shape.points[-1][1])
144 count[a] = count.get(a, 0) + 1
145 count[b] = count.get(b, 0) + 1
146 color_vertices = {k: ('yellow' if v % 2 == 1 else 'black')
147 for k, v in count.items()}
148 import folium
149 map_osm = None
150 for i, index in enumerate(subset):
151 shape = shapes[index]
152 if map_osm is None:
153 if "location" not in kwargs:
154 x, y = shape.points[0]
155 map_osm = folium.Map(location=[y, x], **kwargs)
156 else:
157 map_osm = folium.Map(kwargs["location"], **kwargs)
158 map_osm.add_child(folium.PolyLine(
159 locations=[(_[1], _[0]) for _ in shape.points], weight=10))
160 if color_vertices:
161 a = (shape.points[0][0], shape.points[0][1])
162 b = (shape.points[-1][0], shape.points[-1][1])
163 c1 = color_vertices[a]
164 c2 = color_vertices[b]
165 map_osm.add_child(folium.CircleMarker(
166 [shape.points[0][1], shape.points[0][0]], popup=str((i, index)), radius=8, fill_color=c1))
167 map_osm.add_child(folium.CircleMarker(
168 [shape.points[-1][1], shape.points[-1][0]], popup=str((i, index)), radius=8, fill_color=c2))
169 else:
170 map_osm.add_child(folium.CircleMarker(
171 [shape.points[0][1], shape.points[0][0]], popup=str((i, index)), radius=8, fill_color="yellow"))
172 map_osm.add_child(folium.CircleMarker(
173 [shape.points[-1][1], shape.points[-1][0]], popup=str((i, index)), radius=8, fill_color="yellow"))
174 return folium_html_map(map_osm, width=html_width, height=html_height)
177def build_streets_vertices(edges, shapes):
178 """
179 Returns vertices and edges based on the subset of edges.
181 @param edges indexes
182 @param shapes streets
183 @return vertices, edges
185 *vertices* is a list of points.
186 *edges* is a list of `tuple(a, b)` where `a`, `b` are
187 indices refering to the array of vertices
188 """
189 points = []
190 for i in edges:
191 p = shapes[i].points
192 a, b = (p[0][0], p[0][1]), (p[-1][0], p[-1][1])
193 points.append(a)
194 points.append(b)
195 vertices = list(sorted(set(points)))
196 positions = {p: i for i, p in enumerate(vertices)}
197 new_edges = []
198 for i in edges:
199 points = shapes[i].points
200 a, b = (points[0][0], points[0][1]), (points[-1][0], points[-1][1])
201 new_edges.append((positions[a], positions[b]))
202 return vertices, new_edges
205def plot_streets_network_projection(central_longitude=0.0, min_latitude=-80.0,
206 max_latitude=84.0, globe=None,
207 latitude_true_scale=0.0):
208 """
209 Returns the default projection for @see fn plot_streets_network.
210 See `projection <https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html>`_.
211 """
212 import cartopy.crs as ccrs # pylint: disable=E0401
213 return ccrs.PlateCarree()
216def plot_streets_network(edges_index, edges, vertices, shapes, order=None,
217 color_vertices=None, color_edges=None, ax=None, **kwargs):
218 """
219 Plots the network based on :epkg:`cartopy`.
221 @param edges_index index of the edges in shapes
222 @param edges list of tuple(v1, v2) in array of vertices
223 @param vertices list of vertices coordinates
224 @param shapes streets
225 @param order list of edges composing a path (eulerian path)
226 @param color_vertices dictionary ``{ vertex_index: color }``,
227 changes the color associated to each vertex (black by default)
228 @param color_edges dictionary ``{ edges_index: color }``,
229 changes the color associated to each edge (black by default)
230 @param ax axis or None
231 @param kwargs parameter used to create the plot is ax is None
232 @return ax
234 *kwargs* may contain parameters:
235 *color_v*, *color_e*, *size_v*, *size_e*, *size_c*, *size_et*.
237 If *order* is not None, the function replaces the edge index by its position in this array.
239 if *color_vertices* is equal to `'odd'`, the function computes the degree
240 of each vertex and choose a different color for odd (yellow)
241 and even degrees (black).
242 If *ax* is predefined, it should contain the parameter::
244 import cartopy.crs as ccrs
245 import matplotlib.pyplot as plt
246 from ensae_projects.datainc.data_geo_streets import plot_streets_network_projection
247 fig = plt.figure(figsize=(7,7))
248 ax = fig.add_subplot(1, 1, 1, projection=plot_streets_network_projection())
250 The default projection is given by @see fn plot_streets_network_projection.
251 https://scitools.org.uk/cartopy/docs/v0.15/examples/hurricane_katrina.html
252 """
253 import cartopy.feature as cfeature # pylint: disable=E0401
254 import matplotlib.pyplot as plt
255 from matplotlib.lines import Line2D
256 import numpy
257 from matplotlib.path import Path
259 params = ["color_v", "color_e", "size_v", "size_e",
260 "size_c", "size_et", "projection"]
261 if ax is None:
262 options = {k: v for k, v in kwargs.items() if k not in params}
263 fig = plt.figure(**options)
264 projection = plot_streets_network_projection()
265 ax = fig.add_subplot(1, 1, 1, projection=projection)
266 else:
267 projection = kwargs.get(
268 'projection', plot_streets_network_projection())
270 try:
271 ax.add_feature(cfeature.OCEAN)
272 ax.add_feature(cfeature.COASTLINE)
273 ax.add_feature(cfeature.BORDERS, linestyle=':')
274 except AttributeError as e:
275 raise AttributeError(
276 "cartopy is not properly set up, did you use parameter projection?") from e
278 x1, x2 = min(_[0] for _ in vertices), max(_[0] for _ in vertices)
279 y1, y2 = min(_[1] for _ in vertices), max(_[1] for _ in vertices)
280 dx, dy = (x2 - x1) * 0.2, (y2 - y1) * 0.2
281 x1 -= dx
282 x2 += dx
283 y1 -= dy
284 y2 += dy
285 ax.set_extent([x1, x2, y1, y2], projection)
287 for n in edges_index:
288 sh = shapes[n]
289 geo_points = sh.points
290 lons = [_[0] for _ in geo_points]
291 lats = [_[1] for _ in geo_points]
292 ecolor = color_edges.get(
293 n, "black") if color_edges is not None else "black"
295 linewidth = kwargs.get('size_e', 2)
297 line = Line2D(lons, lats, lw=linewidth, color=ecolor, axes=ax)
298 ax.add_line(line)
300 mx, my = (lons[0] + lons[-1]) / 2, (lats[0] + lats[-1]) / 2
301 gx, gy = mx, my
302 if order is None:
303 ax.text(gx, gy, "e%d" % n, color=kwargs.get('color_e', "blue"))
304 else:
305 pos = [i + 1 for i, v in enumerate(order) if v == n]
306 if len(pos) > 0:
307 pos = [str(_) for _ in pos]
308 ax.text(gx, gy, ",".join(pos),
309 size=kwargs.get("size_et", 12),
310 color=kwargs.get('color_e', "blue"))
312 if color_vertices == "odd":
313 count = {}
314 for a, b in edges:
315 count[a] = count.get(a, 0) + 1
316 count[b] = count.get(b, 0) + 1
317 color_vertices = {k: ('yellow' if v % 2 == 1 else 'black')
318 for k, v in count.items()}
320 theta = numpy.linspace(0, 2 * numpy.pi, 100)
321 circle_verts = numpy.vstack([numpy.sin(theta), numpy.cos(theta)]).T
322 concentric_circle = Path.make_compound_path(Path(circle_verts[::-1]),
323 Path(circle_verts * 0.6))
324 for n, (a, b) in enumerate(vertices):
325 gx, gy = a, b
326 color = color_vertices.get(n, 'black') if color_vertices else 'black'
327 ax.plot(gx, gy, transform=projection,
328 marker=concentric_circle, color=color, markersize=5,
329 linestyle='')
330 ax.text(gx, gy, "v%d" % n, size=kwargs.get('size_v', 12),
331 color=kwargs.get('color_v', "red"))
332 return ax
335def connect_streets(st1, st2):
336 """
337 Tells if streets `st1`, `st2` are connected.
339 @param st1 street 1
340 @param st2 street 2
341 @return tuple or tuple (0 or 1, 0 or 1)
343 Each tuple means:
345 * 0 or 1 mean first or last extremity or the first street
346 * 0 or 1 mean first or last extremity or the second street
348 ``((0, 1),)`` means the first point of the first street is connected
349 to the second extremity of the second street.
350 """
351 a1, b1 = st1[0], st1[-1]
352 a2, b2 = st2[0], st2[-1]
353 connect = []
354 if a1 == a2:
355 connect.append((0, 0))
356 if a1 == b2:
357 connect.append((0, 1))
358 if b1 == a2:
359 connect.append((1, 0))
360 if b1 == b2:
361 connect.append((1, 1))
362 return tuple(connect) if connect else None
365def _complete_subset_streets(edges, shapes):
366 """
367 Extends a set of edges to have less extermities into the graph
368 composed by the sets of edges.
370 @param edges list of indices in shapes
371 @param shapes streets
372 @return added edges (indices)
373 """
374 sedges = set(edges)
375 extension = []
376 for i, _ in enumerate(shapes):
377 if i in sedges:
378 continue
379 add = []
380 for s in edges:
381 if s != i:
382 con = connect_streets(shapes[s].points, shapes[i].points)
383 if con is not None:
384 add.extend([_[1] for _ in con]) # pylint: disable=E1133
385 if len(set(add)) == 2:
386 extension.append(i)
387 return extension
390def _enumerate_close(lon, lat, shapes, th=None):
391 """
392 Enumerates close streets from a starting point.
394 @param lon longitude
395 @param lat latitude
396 @param shapes list of streets
397 @param th threshold or None for all
398 @return iterator on *tuple(distance, i)*
399 """
400 from shapely.geometry import Point, LineString
401 p = Point(lon, lat)
402 for i, shape in enumerate(shapes):
403 obj = LineString(shape.points)
404 d = p.distance(obj)
405 if th is None or d <= th:
406 yield d, i
409def seattle_streets_set_level(shapes, records,
410 pos=(-122.3521425, 47.6219965),
411 size=120):
412 """
413 Returns a graph of streets.
415 @param shapes list of streets
416 @param records description of each street
417 @param pos center of the graphs
418 @param size number of elements
419 @return indices of edges, edges, vertices, distances
421 The function uses the fields ``SEGMENT_TY`` to filter out train, rail ways.
422 """
423 from shapely.geometry import LineString
424 lon, lat = pos
425 closes = list(_enumerate_close(lon, lat, shapes))
426 closes.sort()
427 subset = [index for dist, index in closes[:size]]
428 newset = list(set(subset + _complete_subset_streets(subset, shapes)))
429 # we remove trains lines
430 newset = [n for n in newset if records[n][8] <= 7]
431 vertices, edges = build_streets_vertices(newset, shapes)
432 distances = [LineString(shapes[i].points).length for i in newset]
433 return newset, edges, vertices, distances
436def seattle_streets_set_small(shapes, records, size=20):
437 """
438 Returns a small graph of streets.
440 @param shapes list of streets
441 @param records description of each street
442 @param size number of elements
443 @return indices of edges, edges, vertices, distances
444 """
445 return seattle_streets_set_level(shapes, records, pos=(-122.34651954599997, 47.46947199700003), size=size)
448def seattle_streets_set_level2(shapes, records, size=120):
449 """
450 Returns a small graph of streets.
452 @param shapes list of streets
453 @param records description of each street
454 @param size number of elements
455 @return indices of edges, edges, vertices, distances
456 """
457 return seattle_streets_set_level(shapes, records, pos=(-122.3521425, 47.6219965), size=size)
460def seattle_streets_set_level3(shapes, records, size=1200):
461 """
462 Returns a small graph of streets.
464 @param shapes list of streets
465 @param records description of each street
466 @param size number of elements
467 @return indices of edges, edges, vertices, distances
468 """
469 return seattle_streets_set_level(shapes, records, pos=(-122.400931, 47.648435), size=size)