Skip to content
Snippets Groups Projects
Commit f615d10d authored by Jonas Wassmer's avatar Jonas Wassmer
Browse files

Merge branch 'main' of... 

Merge branch 'main' of https://gitlab.pik-potsdam.de/jonaswa/resilience-of-transportation-infrastructure-networks-to-road-failures
parents b0885796 30927903
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1930 deletions
src/__pycache__/flooding.cpython-310.pyc deleted 100644 → 0
+ 0
− 0
View file @ b0885796
File deleted
src/__pycache__/gravity_model.cpython-310.pyc deleted 100644 → 0
+ 0
− 0
View file @ b0885796
File deleted
src/__pycache__/gravity_model.cpython-39.pyc deleted 100644 → 0
+ 0
− 0
View file @ b0885796
File deleted
src/__pycache__/population.cpython-310.pyc deleted 100644 → 0
+ 0
− 0
View file @ b0885796
File deleted
src/__pycache__/population.cpython-39.pyc deleted 100644 → 0
+ 0
− 0
View file @ b0885796
File deleted
src/old/DaganzoModule.py deleted 100644 → 0
+ 0
− 218
View file @ b0885796
import networkx as nx
import osmnx as ox
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import numpy as np
import warnings
from mpl_toolkits.axes_grid1 import make_axes_locatable
def map_highway_to_number_of_lanes(highway_type):
if highway_type == "motorway" or highway_type == "trunk":
return 4
elif highway_type == "primary":
return 3
elif (
highway_type == "secondary"
or highway_type == "motorway_link"
or highway_type == "trunk_link"
or highway_type == "primary_link"
):
return 2
else:
return 1
# def set_number_of_lanes(G):
# edges = ox.graph_to_gdfs(G, nodes=False)
# lanes = {e : map_highway_to_number_of_lanes(highway_type) for e, highway_type in edges['highway'].items()}
# nx.set_edge_attributes(G, lanes, 'lanes')
# return G
def set_number_of_lanes(G):
edges = ox.graph_to_gdfs(G, nodes=False)
edges["lanes"] = [
np.mean(list(map(float, v[0])))
if type(v[0]) == list
else float(v[0])
if type(v[0]) == str
else map_highway_to_number_of_lanes(v[1])
if np.isnan(v[0])
else v[0]
for k, v in edges[["lanes", "highway"]].iterrows()
]
nx.set_edge_attributes(G, edges["lanes"], "lanes")
return G
def set_actual_speed_diff(G, load_kw="load", alpha=50.0, kph_min=10):
loads = nx.get_edge_attributes(G, load_kw)
lengths = nx.get_edge_attributes(G, "length")
lanes = nx.get_edge_attributes(G, "lanes")
speed_lim = {e: v for e, v in nx.get_edge_attributes(G, "speed_kph").items()}
speeds = {
e: np.divide(alpha * lanes[e] * lengths[e], loads[e]) - 9
for e in G.edges(keys=True)
}
actual_speeds = {
e: kph_min if v < kph_min else speed_lim[e] if v > speed_lim[e] else v
for e, v in speeds.items()
}
speed_diff = {e: (speed_lim[e] - v) for e, v in actual_speeds.items()}
if load_kw == "load":
nx.set_edge_attributes(G, actual_speeds, "actual_speed")
nx.set_edge_attributes(G, speed_diff, "speed_diff")
elif load_kw == "load_flood":
nx.set_edge_attributes(G, actual_speeds, "actual_speed_flood")
nx.set_edge_attributes(G, speed_diff, "speed_diff_flood")
else:
print("wrong load_kw. Returning G")
return G
def set_daganzo_velocities(G, alpha, kph_min):
G = set_number_of_lanes(G)
loads = nx.get_edge_attributes(G, "load")
lengths = nx.get_edge_attributes(G, "length") # length in m
lanes = nx.get_edge_attributes(G, "lanes")
speed_lim = {e: v for e, v in nx.get_edge_attributes(G, "speed_kph").items()}
velocities = {
e: speed_lim[e]
if loads[e] == 0
else 60
* 60
/ 1000
* (alpha * np.divide(lanes[e] * lengths[e], loads[e] * 2) - 5 / 2)
for e in G.edges(keys=True)
}
daganzo_velo = {
e: kph_min
if v < kph_min
else speed_lim[e]
if v > speed_lim[e]
else np.round(v, 1)
for e, v in velocities.items()
}
nx.set_edge_attributes(G, daganzo_velo, "effective_velocity")
# def actual_speed(G, alpha = 50., kph_min = 10):
# loads = nx.get_edge_attributes(G, 'load')
# lengths = nx.get_edge_attributes(G, 'length')
# lanes = nx.get_edge_attributes(G, 'lanes')
# speed_lim = { e : v for e, v in nx.get_edge_attributes(G, 'speed_kph').items()}
# speeds = { e : alpha*lanes[e]*lengths[e]/loads[e] for e in G.edges(keys=True) }
#
# actual_speeds = { e : kph_min if v < kph_min else speed_lim[e] if v > speed_lim[e] else v for e, v in speeds.items() }
#
# nx.set_edge_attributes(G, actual_speeds, 'actual_speed')
# return G
def set_speed_diff(g):
speed = nx.get_edge_attributes(g, "actual_speed")
speed_kph = nx.get_edge_attributes(g, "speed_kph")
speed_diff = {k: speed_kph[k] - v for k, v in speed.items()}
nx.set_edge_attributes(g, speed_diff, "speed_diff")
return g
def split_graph_attributes(
G, attr="speed_diff", vmin=0, vmax=0, cmapsteps=8, cmap="viridis"
):
edges = ox.graph_to_gdfs(G, nodes=False)
vals = pd.Series(nx.get_edge_attributes(G, attr))
cmap = plt.cm.get_cmap(cmap).copy()
cmaplist = [cmap(i) for i in range(cmap.N)]
cmap_discrete = mpl.colors.LinearSegmentedColormap.from_list(
"viridis_discrete", cmaplist, cmap.N
)
cmap_discrete.set_under("dimgrey")
cmap_discrete.set_over("lightgrey")
# define the bins and normalize
if vmin == 0 and vmax == 0:
bounds = np.ceil(np.linspace(min(vals), max(vals), cmapsteps))
else:
bounds = np.ceil(np.linspace(vmin, vmax, cmapsteps))
norm = mpl.colors.BoundaryNorm(bounds, cmap_discrete.N)
scalar_mapper = plt.cm.ScalarMappable(cmap=cmap_discrete, norm=norm)
ec = dict(vals.map(scalar_mapper.to_rgba))
nx.set_edge_attributes(G, ec, "ec")
# split graph according to attributes
# create generators
edg_2 = (
edge
for edge, hw in edges["highway"].items()
if hw == "motorway"
or hw == "trunk"
or hw == "primary"
or hw == "motorway_link"
or hw == "trunk_link"
or hw == "primary_link"
)
edg_1 = (
edge
for edge, hw in edges["highway"].items()
if hw == "secondary"
or hw == "secondary_link"
or hw == "tertiary"
or hw == "tertiary_link"
)
G2 = G.edge_subgraph(edg_2)
G1 = G.edge_subgraph(edg_1)
warnings.filterwarnings("ignore", category=DeprecationWarning)
fig, ax = ox.plot_graph(
G,
figsize=(20, 12),
edge_color=pd.Series(nx.get_edge_attributes(G, "ec")),
edge_linewidth=0.2,
node_size=0,
show=False,
)
ox.plot_graph(
G2,
ax=ax,
figsize=(20, 12),
edge_color=pd.Series(nx.get_edge_attributes(G2, "ec")),
edge_linewidth=3,
node_size=0,
show=False,
)
ox.plot_graph(
G1,
ax=ax,
figsize=(20, 12),
edge_color=pd.Series(nx.get_edge_attributes(G1, "ec")),
edge_linewidth=1,
node_size=0,
show=False,
)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=-3)
cb = fig.colorbar(scalar_mapper, cax=cax)
cb.set_label(r"$\Delta V_{i,j}$ [km/h]", color="w")
cb.ax.yaxis.set_tick_params(color="w")
plt.setp(plt.getp(cb.ax.axes, "yticklabels"), color="w")
cb.ax.tick_params(width=1.0)
src/old/EffectiveVelocities.py deleted 100644 → 0
+ 0
− 234
View file @ b0885796
# %%
import networkx as nx
import numpy as np
import osmnx as ox
import warnings
import matplotlib.pyplot as plt
import matplotlib as mpl
def map_highway_to_number_of_lanes(highway_type):
if highway_type == "motorway" or highway_type == "trunk":
return 4
elif highway_type == "primary":
return 3
elif (
highway_type == "secondary"
or highway_type == "motorway_link"
or highway_type == "trunk_link"
or highway_type == "primary_link"
):
return 2
else:
return 1
def set_number_of_lanes(G):
edges = ox.graph_to_gdfs(G, nodes=False)
edges["lanes"] = [
np.mean(list(map(int, v[0])))
if type(v[0]) == list
else int(v[0])
if type(v[0]) == str
else map_highway_to_number_of_lanes(v[1])
if np.isnan(v[0])
else v[0]
for k, v in edges[["lanes", "highway"]].iterrows()
]
nx.set_edge_attributes(G, edges["lanes"], "lanes")
return G
def peak_loads(G, N_road):
Gc = G.copy()
Gc = set_number_of_lanes(Gc)
pop = sum(nx.get_node_attributes(Gc, "population").values())
loads = nx.get_edge_attributes(Gc, "load")
tot_load = sum(loads.values())
peak_load = {k: N_road * pop * (v / tot_load) for k, v in loads.items()}
nx.set_edge_attributes(Gc, peak_load, "peak_load")
return Gc
def peak_loads_and_capacities(G, N_road):
Gc = G.copy()
Gc = set_number_of_lanes(Gc)
pop = sum(nx.get_node_attributes(Gc, "population").values())
loads = nx.get_edge_attributes(Gc, "load")
lanes = nx.get_edge_attributes(Gc, "lanes")
lanes_mod = dict((k, max(v - 1, 1)) for k, v in lanes.items())
lengths = nx.get_edge_attributes(Gc, "length")
tot_load = sum(loads.values())
peak_load = {k: N_road * pop * (v / tot_load) for k, v in loads.items()}
capacities = {
e: lengths[e] * lanes_mod[e] / (2 * 5 * 1000 / 60 / 60 + 5)
for e in Gc.edges(keys=True)
}
nx.set_edge_attributes(Gc, peak_load, "peak_load")
nx.set_edge_attributes(Gc, capacities, "capacity")
return Gc
def reroute_overloaded_roads(Gc, N_road):
G = Gc.copy()
G = peak_loads_and_capacities(G, N_road)
peak_loads = nx.get_edge_attributes(G, "peak_load")
max_caps = nx.get_edge_attributes(G, "capacity")
if sum(peak_loads.values()) > sum(max_caps.values()):
print(
"Warning: Total loads greater than capacity. Will not converge returning without rerouting."
)
return G
over_idxs = [e for e in G.edges(keys=True) if (peak_loads[e] - max_caps[e]) > 0]
nx.set_edge_attributes(G, peak_loads, "spillover_load")
while len(over_idxs) > 0:
for e in over_idxs:
u, v, k = e
capacity = G[u][v][k]["capacity"]
diff_load = G[u][v][k]["spillover_load"] - capacity
G[u][v][k]["spillover_load"] = capacity # diff_load
pred_edges = G.in_edges(u)
# if len(pred_edges) == 0:
# warnings.warn(f'No predecessor edges to {e} with spillover {diff_load}. Load can not be rerouted. This may impact the results.')
for n, m in pred_edges:
G[n][m][0]["spillover_load"] += diff_load / len(pred_edges)
peak_loads = nx.get_edge_attributes(G, "spillover_load")
over_idxs = [
e for e in G.edges(keys=True) if (peak_loads[e] - max_caps[e]) > 1e-5
]
return G
""""
def reroute_overloaded_roads(G, N_road):
Gc = peak_loads_and_capacities(G, N_road)
peak_load = nx.get_edge_attributes(Gc, 'peak_load')
max_caps = nx.get_edge_attributes(Gc, 'capacity')
if sum(peak_load.values()) > sum(max_caps.values()):
print('Warning: Total loads greater than capacity. Will not converge returning without rerouting.')
return Gc
over_idxs = [ e for e in Gc.edges(keys=True) if (peak_load[e]-max_caps[e]) > 1e-7 ]
nx.set_edge_attributes(Gc, peak_load, 'spillover_load')
for u,v,k in over_idxs:
diff_load = Gc[u][v][k]['spillover_load'] - Gc[u][v][k]['capacity']
pred_edges = nx.bfs_edges(Gc, u, reverse=True)
for (m, n) in pred_edges:
max_comp = Gc[n][m][0]['capacity'] - Gc[n][m][0]['spillover_load']
if max_comp >= diff_load:
Gc[n][m][0]['spillover_load'] += diff_load
Gc[u][v][k]['spillover_load'] -= diff_load
break
elif (max_comp > 0) & (max_comp < diff_load):
Gc[n][m][0]['spillover_load'] += max_comp
Gc[u][v][k]['spillover_load'] -= max_comp
diff_load -= max_comp
spillover_load = nx.get_edge_attributes(Gc, 'spillover_load')
max_caps = nx.get_edge_attributes(Gc, 'capacity')
over_spills = [ (spillover_load[e]-max_caps[e]) for e in Gc.edges(keys=True) if (spillover_load[e]-max_caps[e]) > 1e-7 ]
if len(over_spills)>0:
print('Warning: not all overloaded roads could be rerouted.')
print('Overload:', sum(over_spills))
return Gc
"""
def effective_spillover_velocities(G, N_road):
x_veh = 5
t_react = 2
Gc = reroute_overloaded_roads(G, N_road)
spillover_loads = nx.get_edge_attributes(Gc, "spillover_load")
lanes = nx.get_edge_attributes(Gc, "lanes")
lengths = nx.get_edge_attributes(Gc, "length")
speed_lims = nx.get_edge_attributes(Gc, "speed_kph")
lanes_mod = dict((k, max(v - 1, 1)) for k, v in lanes.items())
eff_velos = {
e: (
lengths[e] * lanes_mod[e] / (spillover_loads[e] * t_react) - x_veh / t_react
)
* 60
* 60
/ 1000
for e in Gc.edges(keys=True)
}
eff_velos = dict(
(e, speed_lims[e]) if v > speed_lims[e] else (e, np.round(v, 1))
for e, v in eff_velos.items()
)
rr_travel_time = {
e: np.round(np.divide(lengths[e] / 1000 * 60 * 60, eff_velos[e]), 1)
for e in Gc.edges(keys=True)
}
nx.set_edge_attributes(Gc, rr_travel_time, "spillover_travel_time")
nx.set_edge_attributes(Gc, eff_velos, "spillover_velocity")
return Gc
def effective_velocities(G, N_road):
x_veh = 5
t_react = 2
Gc = peak_loads_and_capacities(G, N_road)
peak_loads = nx.get_edge_attributes(Gc, "peak_load")
lanes = nx.get_edge_attributes(Gc, "lanes")
lengths = nx.get_edge_attributes(Gc, "length")
speed_lims = nx.get_edge_attributes(Gc, "speed_kph")
lanes_mod = dict((k, max(v - 1, 1)) for k, v in lanes.items())
eff_velos = {
e: (lengths[e] * lanes_mod[e] / (peak_loads[e] * t_react) - x_veh / t_react)
* 60
* 60
/ 1000
for e in Gc.edges(keys=True)
}
eff_velos = dict(
(e, speed_lims[e])
if v > speed_lims[e]
else (e, 5)
if v < 5
else (e, np.round(v, 1))
for e, v in eff_velos.items()
)
travel_time = {
e: np.round(np.divide(lengths[e] / 1000 * 60 * 60, eff_velos[e]), 1)
for e in Gc.edges(keys=True)
}
nx.set_edge_attributes(Gc, travel_time, "effective_travel_time")
nx.set_edge_attributes(Gc, eff_velos, "effective_velocity")
return Gc
# %%
src/old/EmergencyModule.py deleted 100644 → 0
+ 0
− 317
View file @ b0885796
import osmnx as ox
import networkx as nx
import pandas as pd
# import geopandas as gpd
import numpy as np
# import matplotlib.pyplot as plt
import os
from pathlib import Path
from shapely.strtree import STRtree
# import src.PopulationModule as pm
import src.GravityModule2 as gm
def add_emergencies(G, north, east, south, west, date, key, tag):
nodes, edges = ox.graph_to_gdfs(G)
# check if emergency-filtered file osm exists already
if tag == "*":
filetag = "all"
else:
filetag = tag
if not Path(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf"
).is_file():
os.system(f"mkdir data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}")
os.system(
f"osmium extract -b {west},{south},{east},{north} data/historic-data.nosync/{date}/germany.osm.pbf\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf --overwrite"
)
if not Path(
f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm"
).is_file():
os.system(
f"osmium tags-filter data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf\
{key}={tag}\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm.pbf --overwrite"
)
os.system(
f"osmium cat data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm.pbf -o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm.bz2 --overwrite"
)
os.system(
f"bzip2 -d data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm.bz2 -f"
)
# get fire-stations/hospitals
df = ox.geometries_from_xml(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{key}-{filetag}.osm"
)
if len(df) == 0:
print(f"No {key}_{tag} in specified area.")
return
emergencies = df[~df[key].isna()]
# set true/false for emergency nodes
if key not in nodes.columns:
nodes[key] = np.full(len(nodes), False)
emerg_geos = emergencies.geometry
node_tree = STRtree(nodes.geometry)
for j, geo in enumerate(emerg_geos):
# find nearest emegreny to road i
emergency_point = geo.centroid
idx = node_tree.nearest(emergency_point)
# set emergency
nodes.iloc[idx, nodes.columns.get_loc(key)] = emergencies[key].iloc[j]
nx.set_node_attributes(G, nodes[key], name=key)
def emergency_graph_from_osmfile(
osm_file, emergency, west, south, east, north, date=None, roads="all"
):
# check if unfiltered osm-file with respected bounding box exists
G = gm.graph_from_osmfile(
osm_file, west, south, east, north, date=date, roads=roads
)
# check if emergency-filtered file osm exists already
if not Path(
f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm"
).is_file():
os.system(
f"osmium tags-filter data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf\
amenity={emergency}\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{emergency}.osm.pbf --overwrite"
)
os.system(
f"osmium cat data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{emergency}.osm.pbf -o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{emergency}.osm.bz2 --overwrite"
)
os.system(
f"bzip2 -d data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{emergency}.osm.bz2 -f"
)
nodes, edges = ox.graph_to_gdfs(G)
# get fire-stations/hospitals
df = ox.geometries_from_xml(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{emergency}.osm"
)
emergencies = df[df["amenity"] == emergency]
# set true/false for emergency nodes
nodes["emergency"] = np.full(len(nodes), False)
emerg_geos = emergencies.geometry
node_tree = STRtree(nodes.geometry)
index_by_id = dict((id(pt), i) for i, pt in enumerate(nodes.geometry))
for geo in emerg_geos:
# find nearest emegreny to road i
emergency_point = geo.centroid
nearest_node = node_tree.nearest(emergency_point)
# find index and osmid
idx = index_by_id[id(nearest_node)]
# n_osmid = nodes.iloc[idx].name
# set emergency to True
nodes.iloc[idx, nodes.columns.get_loc("emergency")] = True
G = ox.graph_from_gdfs(nodes, edges)
# impute speed on all edges missing data
G = ox.add_edge_speeds(G)
# calculate travel time (seconds) for all edges
G = ox.add_edge_travel_times(G)
# population and population per box
# city_df = pd.read_csv("data/city_pop_clean.csv", index_col=0)
pop = pm.pop_guess_from_poi(north, south, east, west)
# add bbox to graph
# boxcount = pm.nodes_in_grid_bbox(G, 50)
R = gm.residential_graph(date, north, south, west, east)
# G = pm.population_from_voronoi_polys(G, pop, boxcount)
G = pm.population_from_voronoi_polys(G, R, pop, box_bins=10)
G.bbox = [north, west, south, east]
return G
def set_shortest_path_to_emergency(G, weight="travel_time"):
nodes = ox.graph_to_gdfs(G, edges=False)
path_length_emerg = nx.multi_source_dijkstra_path_length(
G, list(nodes[nodes["emergency"]].index), weight=weight
)
path_length_emerg = {k: v / 60 for k, v in path_length_emerg.items()}
# if 'flood' in weight:
# nx.set_node_attributes(G, path_length_emerg, 'shortest_path_to_emergency_flood')
# else:
# nx.set_node_attributes(G, path_length_emerg, 'shortest_path_to_emergency')
return path_length_emerg
# def add_emergencies_to_gdfs(nodes, edges, emergency_df):
# """
# DEPRECATED!!!
# Add fire-stations or hospitals to the nodes and edges.
# Split existing edges connecting the new emergency roads in a smart way.
# May only return an UNDIRECTED graph, as only one road is split.
# """
#
# #warnings.filterwarnings("ignore", category=FutureWarning) #ignore deprecation warnings
#
# nodes.loc[:, 'ref'] = 0
# #edges_updated = edges.copy()
#
# edg_tree = STRtree(edges.geometry)
# index_by_id = dict((id(pt), i) for i, pt in enumerate(edges.geometry))
#
#
# for i, v in emergency_df.iterrows():
# #find nearest emegreny to road i
# emergency_point = v.geometry.centroid
# nearest_road = edg_tree.nearest(emergency_point)
# #find index and edges u,v,k
# idx = index_by_id[id(nearest_road)]
# e_idx = edges.iloc[idx].name
#
# #split road
# p, q = nearest_points(nearest_road, emergency_point)
#
# #new road to emergency
# emergency_linestring = LineString([p, q])
#
# #update indices
# p_id = id(p)
# q_id = id(q)
# pi_id = (e_idx[0], p_id, 0)
# pj_id = (p_id, e_idx[1], 0)
# pq_id = (p_id, q_id, 0)
#
# #add emergency-linestring to df
# edg_e = edges.iloc[idx].copy()
# edg_e.name = pq_id
# edg_e['osmid'] = id(emergency_linestring)
# edg_e['length'] *= emergency_linestring.length / nearest_road.length
# edg_e['geometry'] = emergency_linestring
# #edg_e['from'] = q_id
# #edg_e['to'] = p_id
#
# #update nodes
# nd_p = nodes.iloc[0].copy()
# nd_q = nd_p.copy()
#
# nd_p.name = p_id
# nd_q.name = q_id
#
# nd_p['x'] = p.x
# nd_p['y'] = p.y
# nd_p['geometry'] = p
#
# nd_q['x'] = q.x
# nd_q['y'] = q.y
# nd_q['ref'] = 'emergency'
# nd_q['geometry'] = q
#
# nodes = nodes.append([nd_p, nd_q])
#
# #split edge
# #allow small inaccuracies at split
# buff = p.buffer(0.000001)
# road_split = split(nearest_road, buff)
#
# #check if road was split, or if nearest point to emergency is at one endpoint of the road
# #this can be done by comparing the max split lenth to the nearest road_geometry length.
# #if road was split, update dataframes
# split_lengths = (road.length for road in road_split)
# split_compare = abs(max(split_lengths) - nearest_road.length)
# if split_compare > 1e-5:
# first_seg, buff_seg, last_seg = road_split
# line = LineString(list(first_seg.coords) + list(p.coords) + list(last_seg.coords))
# a, b = split(line, p)
#
# #update dataframes with split edges
# edg_a = edges.iloc[idx].copy()
# edg_a.name = pi_id
# edg_a['osmid'] = id(a)
# edg_a['length'] *= a.length/(a.length+b.length)
# edg_a['geometry'] = a
#
# edg_b = edges.iloc[idx].copy()
# edg_b.name = pj_id
# edg_b['osmid'] = id(b)
# edg_b['length'] *= b.length/(a.length+b.length)
# edg_b['geometry'] = b
#
# edges = edges.append([edg_a, edg_b, edg_e])
#
# #drop old edge
# edges = edges.drop(e_idx)
#
# return nodes, edges
# def prepare_emergency_graph(osm_file, emergency, west, south, east, north, GCC=True):
# #warnings.filterwarnings("ignore", category=FutureWarning) #ignore deprecation warnings
#
# #check if unfiltered osm-file with respected bounding box exists
# if not Path(f"data/cache/osmfiles/{north}-{south}-{west}-{east}/nofilter.osm.pbf").is_file():
# os.system(f"mkdir data/cache/osmfiles/{north}-{south}-{west}-{east}")
# os.system(f"osmium extract -b {west},{south},{east},{north} {osm_file} -o data/cache/osmfiles/{north}-{south}-{west}-{east}/nofilter.osm.pbf --overwrite")
#
# #check if highway-filtered file osm exists already
# if not Path(f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm").is_file():
# os.system(f"osmium tags-filter data/cache/osmfiles/{north}-{south}-{west}-{east}/nofilter.osm.pbf w/highway=motorway,trunk,primary,secondary,tertiary,unclassified,residential,motorway_link,trunk_link,primary_link,secondary_link,tertiary_link,living_street -o data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm.pbf --overwrite")
# os.system(f"osmium cat data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm.pbf -o data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm.bz2 --overwrite")
# os.system(f"bzip2 -d data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm.bz2 -f")
#
#
# #prepare graph from osmfile
# G = ox.graph_from_xml(f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-highway.osm", retain_all=True, simplify=True)
#
# if GCC:
# G = ox.utils_graph.get_largest_component(G, strongly=True)
#
#
# #check if emergency file exists alreadt
# if not Path(f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm").is_file():
#
# os.system(f"osmium tags-filter data/cache/osmfiles/{north}-{south}-{west}-{east}/nofilter.osm.pbf emergency={emergency} -o data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm.pbf --overwrite")
# os.system(f"osmium cat data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm.pbf -o data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm.bz2 --overwrite")
# os.system(f"bzip2 -d data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm.bz2 -f")
#
# #get df with emergencies
# emergency_df = ox.geometries_from_xml(f"data/cache/osmfiles/{north}-{south}-{west}-{east}/filter-{emergency}.osm")
#
#
# U = ox.get_undirected(G)
#
#
# nodes, edges = ox.graph_to_gdfs(U, fill_edge_geometry=True)
#
# #add emergency nodes to graph and split nereast edges to link them to existing graph
# nodes, edges = add_emergencies_to_gdfs(nodes, edges, emergency_df)
#
# G = ox.graph_from_gdfs(nodes, edges)
# G = ox.get_undirected(G)
#
# # impute speed on all edges missing data
# G = ox.add_edge_speeds(G)
# # calculate travel time (seconds) for all edges
# G = ox.add_edge_travel_times(G)
# return G
src/old/GermanMobiltyPanel.py deleted 100644 → 0
+ 0
− 59
View file @ b0885796
#%%
import pandas as pd
import numpy as np
import os
from scipy.optimize import curve_fit
#path = 'data/MOP-data/mobility/'
def read_PT_kmPKW_files(path, attr):
data_points = np.array([])
frames = []
for subdir, dirs, files in os.walk(path):
for file in files:
if ('PT' in file) and ('.csv' in file):
PT = pd.read_csv(os.path.join(subdir, file), sep=';')
PT_clean = PT[(PT['AnzPKW'] != 0) & (PT['dau_PKW'] > 0) & (PT['km_PKW'] > 0)]
frames.append(PT_clean)
PTS = pd.concat(frames)
if attr == 'travel_time':
attr = 'dau_PKW'
elif attr == 'length':
attr = 'km_PKW'
q = PTS[attr].quantile(0.75)
filtered_PT = PTS[PTS[attr] < q]
data_points = (filtered_PT[attr]/filtered_PT['AnzPKW']).to_numpy()
return data_points
#%%
def exp_func(x, a, k):
return a * np.exp(-k * x)
def lin_func(x, m, b):
return m * x + b
def mobility_fit_params(path, attr, bincount=250):
if attr == 'length':
multiplier = 1000 #get m from km
elif attr == 'travel_time':
multiplier = 60 #get sec from min
data_points = read_PT_kmPKW_files(path, attr)*multiplier
data_pkw, bins = np.histogram(data_points, bins=bincount, density=True)
binscenters = np.array([0.5 * (bins[i] + bins[i+1]) for i in range(len(bins)-1)])
max_ind = np.argmax(data_pkw)
popt_exp, pcov_exp = curve_fit(exp_func, binscenters[max_ind:], data_pkw[max_ind:], p0=[0.05, 0.001])
popt_lin, pcov_lin = curve_fit(lin_func, binscenters[:max_ind], data_pkw[:max_ind], p0=[0.01, 0.01])
return binscenters[max_ind], popt_exp, popt_lin
def mobility_func(xs, max_bin, popt_exp, popt_lin):
if hasattr(xs, "__len__"):
return [ exp_func(x, *popt_exp) if x>max_bin else lin_func(x, *popt_lin) for x in xs ]
elif xs>max_bin:
return exp_func(xs, *popt_exp)
else:
return lin_func(xs, *popt_lin)
\ No newline at end of file
src/old/GravityModule2.py deleted 100644 → 0
+ 0
− 157
View file @ b0885796
# %%
import networkx as nx
import osmnx as ox
import os
import numpy as np
import pickle
import time
import multiprocessing
from pathlib import Path
from tqdm import tqdm
import src.GermanMobiltyPanel as gmp
def flows_od(pop_o, pops, dists, mobility_fit):
"""Returns dict with flows from origin (o) to all destinations (d). {d : f_o}
pop_o : population at origin o
pops : population dictionary
dists : distances from o to destinations
"""
denominator = sum(mobility_fit(dist) * pops[l] for (l, dist) in dists.items())
return {
d: pop_o * pops[d] * mobility_fit(dist) / denominator
for (d, dist) in dists.items()
}
def flows(G, pops, mobility_fit, path_attr, cutoff=None):
"""Returns generator with flows from all origins to all destinations. (o : {d} : f_o})
G : road graph
path_attr : edge weight for shortest path algorithm.
"""
return (
(flows_od(pops[o], pops, dists, mobility_fit), paths)
for o, (dists, paths) in nx.all_pairs_dijkstra(
G, cutoff=cutoff, weight=path_attr
)
)
def compute_loads_OLD(
graph, path_attr="travel_time", disable_tqdm=False, cutoff="default"
):
"""Return graph object with commuter loads as edge weight attr on all edges."""
G = graph.copy()
max_bin, popt_exp, popt_lin = gmp.mobility_fit_params(
"data/MOP-data/mobility/", path_attr, bincount=250
)
mobility_fit = (
lambda x: gmp.exp_func(x, *popt_exp)
if x > max_bin
else gmp.lin_func(x, *popt_lin)
)
L = dict(zip(G.edges(keys=True), np.zeros(len(G.edges()))))
pops = nx.get_node_attributes(G, "population")
if cutoff == "default":
if path_attr == "length":
cutoff = 60 * 1000 # 60 km
elif path_attr == "travel_time":
cutoff = 60 * 60 # 60 mins
if not disable_tqdm:
print("Computing loads...")
for fo, paths in tqdm(
flows(G, pops, mobility_fit, path_attr, cutoff=cutoff),
total=len(G.nodes()),
disable=disable_tqdm,
):
for d, path in paths.items():
f_od = fo[d]
for i, j in zip(path[:-1], path[1:]):
L[(i, j, 0)] += f_od
nx.set_edge_attributes(G, L, "load")
return G
def compute_loads(graph, weight="travel_time", cutoff="default", disable_tqdm=False):
"""Return graph object with commuter loads as edge weight attr on all edges."""
G = graph.copy()
max_bin, popt_exp, popt_lin = gmp.mobility_fit_params(
"data/MOP-data/mobility/", weight, bincount=250
)
mobility_fit = (
lambda x: gmp.exp_func(x, *popt_exp)
if x > max_bin
else gmp.lin_func(x, *popt_lin)
)
L = dict(zip(G.edges(keys=True), np.zeros(len(G.edges()))))
populations = nx.get_node_attributes(G, "population")
if cutoff == "default":
if weight == "length":
cutoff = 60 * 1000 # 60 km
elif weight == "travel_time":
cutoff = 60 * 60 # 60 mins
if not disable_tqdm:
print("Computing loads...")
for o, (dist, path) in tqdm(
nx.all_pairs_dijkstra(G, cutoff=cutoff, weight=weight),
total=len(G.nodes()),
disable=disable_tqdm,
):
No = populations[o]
denominator = sum(
mobility_fit(dist) * populations[l] for (l, dist) in dist.items()
)
for d, path_od in path.items():
Nd = populations[d]
dist_od = dist[d]
fod = No * Nd * mobility_fit(dist_od) / denominator
L.update(
{
(i, j, 0): L[(i, j, 0)] + fod
for i, j in zip(path_od[:-1], path_od[1:])
}
)
# for i, j in zip(path_od[:-1], path_od[1:]):
# L[(i, j, 0)] += fod
nx.set_edge_attributes(G, L, "load")
return G
def commuter_loads(G, path_attr="travel_time", cache=True, cutoff="default"):
"""Return graph object with commuter loads as edge weight attr on all edges."""
if cache:
hash = nx.weisfeiler_lehman_graph_hash(
ox.get_digraph(G),
edge_attr=path_attr,
node_attr="population",
iterations=10,
digest_size=32,
)
path = f"data/cache/load-files/{hash}.pkl"
if os.path.isfile(path):
print(
"The load has already been calculated and saved. Returning this output."
)
with open(path, "rb") as f:
load = pickle.load(f)
nx.set_edge_attributes(G, load, "load")
# return G
else:
G = compute_loads(G, weight=path_attr, cutoff=cutoff, disable_tqdm=False)
load = nx.get_edge_attributes(G, "load")
with open(path, "wb") as f:
pickle.dump(load, f)
# return G
else:
G = compute_loads(G, weight=path_attr, cutoff=cutoff, disable_tqdm=False)
return G
src/old/ParallelLoad.py deleted 100644 → 0
+ 0
− 173
View file @ b0885796
# %%
from multiprocess import Pool
import multiprocess
import networkx as nx
import numpy as np
import time
import os
import osmnx as ox
import pickle
from src import GermanMobiltyPanel as gmp
# %%
def parallel_load_with_cache(
graph,
cpu_cores=5,
cache=True,
jobs_per_cpu=5,
cutoff="default",
weight="travel_time",
):
if cache:
hash = nx.weisfeiler_lehman_graph_hash(
ox.get_digraph(graph),
edge_attr=weight,
node_attr="population",
iterations=10,
digest_size=32,
)
path = f"data/cache/load-files/{hash}.pkl"
if os.path.isfile(path):
print("Load already computed. Returning saved values.")
with open(path, "rb") as f:
load = pickle.load(f)
nx.set_edge_attributes(graph, load, "load")
else:
graph = parallel_load(
graph,
cpu_cores=cpu_cores,
jobs_per_cpu=jobs_per_cpu,
weight=weight,
cutoff=cutoff,
)
load = nx.get_edge_attributes(graph, "load")
with open(path, "wb") as f:
pickle.dump(load, f)
else:
graph = parallel_load(
graph,
cpu_cores=cpu_cores,
jobs_per_cpu=jobs_per_cpu,
weight=weight,
cutoff=cutoff,
)
return graph
def parallel_load(
graph, cpu_cores=5, jobs_per_cpu=5, cutoff="default", weight="travel_time"
):
global load_to_orig
if cutoff == "default":
if weight == "length":
cutoff = 60 * 1000 # 60 km
elif weight == "travel_time":
cutoff = 60 * 60 # 60 mins
max_bin, popt_exp, popt_lin = gmp.mobility_fit_params(
"data/MOP-data/mobility/", "travel_time", bincount=250
)
mobility_fit = (
lambda x: gmp.exp_func(x, *popt_exp)
if x > max_bin
else gmp.lin_func(x, *popt_lin)
)
population = nx.get_node_attributes(graph, "population")
zeros = np.zeros(len(graph.edges), dtype=np.float64)
def load_to_orig(origin_list):
# print(multiprocessing.current_process())
loaddict = dict(zip(graph.edges(keys=True), zeros))
for o in origin_list:
dist, pathes = nx.single_source_dijkstra(
graph, o, weight=weight, cutoff=cutoff
)
No = population[o]
denominator = sum(
mobility_fit(dist) * population[l] for (l, dist) in dist.items()
)
for d, path in pathes.items():
Nd = population[d]
fod = No * Nd * mobility_fit(dist[d]) / denominator
for i, j in zip(path[:-1], path[1:]):
loaddict[i, j, 0] += fod
data = np.array(list(loaddict.values()))
return data
# cpu_count() will report /all/ CPUs on the node
# This is not what we should use.
ncpus = multiprocess.cpu_count()
print("detected {} cores".format(ncpus))
# This will report the number of CPU cores SLURM
# has allocated us. This is the correct number to
# pass to Pool()
try:
ncpus = int(os.environ["SLURM_JOB_CPUS_PER_NODE"])
# ncpus = int(os.environ.get("SLURM_CPUS_PER_TASK"))
print("my Slurm allocation is {} cores".format(ncpus))
except KeyError:
ncpus = cpu_cores
print(f"Not running under Slurm, setting ncpus to {ncpus}")
node_arr = np.array(list(graph.nodes()))
input_list_len = ncpus * jobs_per_cpu
vertices_per_proc = int(np.floor(len(node_arr) / input_list_len))
node_list_proc = node_arr[: input_list_len * vertices_per_proc].reshape(
input_list_len, vertices_per_proc
)
# the last column can be shorter then the rest (will be used first and used for mp eta)
if len(node_arr[input_list_len * vertices_per_proc :]) > 0:
nodes_last_col = node_arr[input_list_len * vertices_per_proc :]
else:
nodes_last_col = node_list_proc[-1]
# eta
start = time.time()
load_last_col = load_to_orig(nodes_last_col)
end = time.time()
est_time = round((end - start) * len(graph.nodes()) / (ncpus * len(nodes_last_col)))
print(f"Estimated time order: O({est_time} s)")
# start multiprocessing
start = time.time()
with Pool(processes=ncpus) as pool:
# future = pool.map(parallel_loads, list(G.nodes()))
future = pool.map_async(load_to_orig, node_list_proc[:-1])
load_arr = np.sum(future.get(), axis=0)
nx.set_edge_attributes(
graph,
dict(zip(graph.edges(keys=True), load_arr + load_last_col)),
"load",
)
# loadp = nx.get_edge_attributes(G, "parallel_load")
end = time.time()
print("Time:", round(end - start, 1), "seconds")
return graph
# %%
""""
west, east = 6.75, 7.35
south, north = 50.35, 50.65
date = "2021-01-02"
region = rn.RoadNetwork(north, east, south, west, date, "all")
region.loads("travel_time")
H = region.graph
H = parallel_load(H, 40, cutoff=None)
pload = nx.get_edge_attributes(H, "parallel_load")
load = nx.get_edge_attributes(H, "load")
# %%
#plt.scatter(load.values(), pload.values())
# %%
"""
src/old/PopulationFromRaster.py deleted 100644 → 0
+ 0
− 202
View file @ b0885796
# %%
import networkx as nx
import warnings
from rasterstats import zonal_stats
import osmnx as ox
import numpy as np
import rioxarray as riox
from rasterio.enums import Resampling
import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
from shapely.geometry import MultiPoint, Point, Polygon
import geopandas as gpd
def box_to_poly(long0, lat0, lat1, long1):
return Polygon([[long0, lat0], [long1, lat0], [long1, lat1], [long0, lat1]])
def voronoi_finite_polygons_2d(vor, radius=None):
if vor.points.shape[1] != 2:
raise ValueError("Requires 2D input")
new_regions = []
new_vertices = vor.vertices.tolist()
center = vor.points.mean(axis=0)
if radius is None:
radius = vor.points.ptp().max()
# Construct a map containing all ridges for a given point
all_ridges = {}
for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
all_ridges.setdefault(p1, []).append((p2, v1, v2))
all_ridges.setdefault(p2, []).append((p1, v1, v2))
# Reconstruct infinite regions
for p1, region in enumerate(vor.point_region):
vertices = vor.regions[region]
if all(v >= 0 for v in vertices):
# finite region
new_regions.append(vertices)
continue
# reconstruct a non-finite region
ridges = all_ridges[p1]
new_region = [v for v in vertices if v >= 0]
for p2, v1, v2 in ridges:
if v2 < 0:
v1, v2 = v2, v1
if v1 >= 0:
# finite ridge: already in the region
continue
# Compute the missing endpoint of an infinite ridge
t = vor.points[p2] - vor.points[p1] # tangent
t /= np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # normal
midpoint = vor.points[[p1, p2]].mean(axis=0)
direction = np.sign(np.dot(midpoint - center, n)) * n
far_point = vor.vertices[v2] + direction * radius
new_region.append(len(new_vertices))
new_vertices.append(far_point.tolist())
# sort region counterclockwise
vs = np.asarray([new_vertices[v] for v in new_region])
c = vs.mean(axis=0)
angles = np.arctan2(vs[:, 1] - c[1], vs[:, 0] - c[0])
new_region = np.array(new_region)[np.argsort(angles)]
# finish
new_regions.append(new_region.tolist())
return new_regions, np.asarray(new_vertices)
def compoute_voronoi_polys_of_nodes(nodes):
vor_nodes = nodes.copy()
points = np.array(list(zip(np.array(vor_nodes.x), np.array(vor_nodes.y))))
vor = Voronoi(points)
regions, vertices = voronoi_finite_polygons_2d(vor)
pts = MultiPoint([Point(i) for i in points])
mask = pts.convex_hull
voronoi_polys = gpd.GeoSeries([Polygon(vertices[region]) for region in regions])
voronoi_polys = voronoi_polys.intersection(mask)
vor_nodes["voronoi"] = list(voronoi_polys)
vor_nodes = vor_nodes.set_geometry("voronoi").set_crs(nodes.crs)
return vor_nodes
def clip_and_upsample_to_gdf(voronoi_nodes, raster_path):
# Read raster
raster = riox.open_rasterio(raster_path)
# update crs
crs = raster.spatial_ref.crs_wkt
voronoi_nodes = voronoi_nodes.to_crs(crs)
# clip to convex hull of gdf
clipped = raster.rio.clip_box(*voronoi_nodes.unary_union.bounds)
clipped = clipped.rio.clip([voronoi_nodes.unary_union.convex_hull])
# compute updscale factors
y_res = int(
np.floor(
min(
voronoi_nodes.voronoi.bounds["maxy"]
- voronoi_nodes.voronoi.bounds["miny"]
)
/ 2
)
)
x_res = int(
np.floor(
min(
voronoi_nodes.voronoi.bounds["maxx"]
- voronoi_nodes.voronoi.bounds["minx"]
)
/ 2
)
)
# upsample raster
up_sampled = clipped.rio.reproject(
clipped.rio.crs,
resolution=(y_res, x_res),
resampling=Resampling.sum,
nodata=np.nan,
)
return up_sampled
def pop_in_polygon(raster, geom):
clipped_poly = raster.rio.clip([geom], drop=True)
clipped_poly = clipped_poly.where(clipped_poly != clipped_poly.rio.nodata)
return float(clipped_poly.sum(skipna=True).data)
def population_from_raster_to_gdf(
gdf,
raster_path="/Users/jonas/Documents/PHD/osm-population/data/GHS/GHS_POP_P2030_GLOBE_R2022A_54009_100_V1_0.tif",
):
org_crs = gdf.crs
upsampled_clip = clip_and_upsample_to_gdf(gdf, raster_path)
# update crs
crs = upsampled_clip.spatial_ref.crs_wkt
gdf = gdf.to_crs(crs)
warnings.filterwarnings("ignore")
stats = zonal_stats(
gdf["voronoi"],
upsampled_clip.data[0],
affine=upsampled_clip.rio.transform(),
stats="sum",
nodata=np.nan,
all_touched=False,
)
gdf["population"] = [s["sum"] for s in stats]
if gdf["population"].isna().sum() > 0:
print("At least one nan value for population. Setting na to zero")
gdf.loc[gdf["population"].isna(), "population"] = 0
pop_dens = gdf["population"] / (gdf["voronoi"].area * 1e-6)
vor_area = gdf["voronoi"].area * 1e-6
gdf["population_density"] = pop_dens
gdf["voronoi_area"] = vor_area
gdf = gdf.to_crs(org_crs)
return gdf
def population_from_raster_to_graph(
G,
raster_path="/Users/jonas/Documents/PHD/osm-population/data/GHS/GHS_POP_P2030_GLOBE_R2022A_54009_100_V1_0.tif",
return_voronoi=False,
):
nodes = ox.graph_to_gdfs(G, edges=False)
voronoi_nodes = compoute_voronoi_polys_of_nodes(nodes)
voronoi_nodes = population_from_raster_to_gdf(voronoi_nodes, raster_path)
nx.set_node_attributes(G, voronoi_nodes["population"], "population")
nx.set_node_attributes(G, voronoi_nodes["population_density"], "population_density")
nx.set_node_attributes(G, voronoi_nodes["voronoi"], "voronoi")
if return_voronoi:
return G, voronoi_nodes
return G
# %%
# G = ox.graph_from_place('Cologne Germany', network_type='drive')
# G, vor_nodes = population_from_raster_to_graph(G, return_voronoi=True)
src/old/RoadNetwork.py deleted 100644 → 0
+ 0
− 570
View file @ b0885796
# %%
import sys
import importlib
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
import os
import pickle
import numpy as np
from pathlib import Path
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import colors
from shapely.geometry import Polygon, LineString, Point
# import src.GermanMobiltyPanel as gmp
import src.EmergencyModule as em
import src.GravityModule2 as gm
import src.DaganzoModule as dm
import src.FloodModule as fm
import src.Isochrones as im
import src.EffectiveVelocities as ev
import src.ParallelLoad as pl
import src.PopulationFromRaster as pfr
# %%
class RoadNetwork:
def __init__(self, north, east, south, west, date, roads, flooded=False, *args):
self.bbox = [north, east, south, west]
self.date = date
self.roads = roads
self.flooded = flooded
self.graph = graph_from_osmfile(north, east, south, west, date, roads)
self.graph = dm.set_number_of_lanes(self.graph)
nodes, edges = ox.graph_to_gdfs(self.graph)
self.nodes = nodes
self.edges = edges
self.population = nodes["population"].sum()
self.area = (
self.nodes.set_geometry("voronoi", crs=self.nodes.crs)
.to_crs("ESRI:54009")
.area.sum()
* 1e-6
)
if flooded:
remove_flooded_edges(self)
update_nodes_edges(self)
def loads(self, path_attr, cache=True, cutoff="default"):
self.graph = gm.commuter_loads(
self.graph, path_attr=path_attr, cache=cache, cutoff=cutoff
)
update_nodes_edges(self)
def parallel_load(
self, path_attr, cpu_cores=5, jobs_per_cpu=5, cache=True, cutoff="default"
):
self.graph = pl.parallel_load_with_cache(
self.graph,
cpu_cores=cpu_cores,
cache=cache,
jobs_per_cpu=jobs_per_cpu,
cutoff=cutoff,
weight=path_attr,
)
update_nodes_edges(self)
def add_emergencies(self, key, tag):
north, east, south, west = self.bbox
date = self.date
em.add_emergencies(self.graph, north, east, south, west, date, key, tag)
update_nodes_edges(self)
def effective_spillover_velocities(self, N_road):
self.graph = ev.effective_spillover_velocities(self.graph, N_road)
update_nodes_edges(self)
def effective_velocities(self, N_road):
self.graph = ev.effective_velocities(self.graph, N_road)
update_nodes_edges(self)
def daganzo_method(self, alpha, kph_min):
if not "load" in self.edges:
raise NameError("Run loads() method before daganzo method.")
dm.set_daganzo_velocities(self.graph, alpha, kph_min)
edges = ox.graph_to_gdfs(self.graph, nodes=False)
# edges = self.edges
# convert distance meters to km, and speed km per hour to km per second
distance_km = edges["length"] / 1000
speed_km_h = edges["effective_velocity"]
# calculate edge travel time in seconds
travel_time = (distance_km / speed_km_h) * 60 * 60
# add travel time attribute to graph edges
edges["effective_travel_time"] = travel_time.round(1).values
nx.set_edge_attributes(
self.graph,
values=edges["effective_travel_time"],
name="effective_travel_time",
)
update_nodes_edges(self)
def shortest_path_from_emergency(self, weight, key, tag, cutoff=None):
G = self.graph
nodes = self.nodes
emergency_nodes = nodes[nodes[key] == tag]
# list(nodes[nodes[emergency]].index),
path_length_emerg = nx.multi_source_dijkstra_path_length(
G, list(emergency_nodes.index), weight=weight, cutoff=cutoff
)
# path_length_emerg = { k : v for k, v in path_length_emerg.items() }
self.nodes[f"{weight}_from_{key}_{tag}"] = path_length_emerg
nx.set_node_attributes(
self.graph, path_length_emerg, f"{weight}_from_{key}_{tag}"
)
def shortest_path_to_emergency(self, weight, key, tag, cutoff=None):
G_reversed = self.graph.reverse()
nodes = self.nodes
emergency_nodes = nodes[nodes[key] == tag]
path_length_emerg = nx.multi_source_dijkstra_path_length(
G_reversed, list(emergency_nodes.index), weight=weight, cutoff=cutoff
)
# path_length_emerg = { k : v for k, v in path_length_emerg.items() }
self.nodes[f"{weight}_to_{key}_{tag}"] = path_length_emerg
nx.set_node_attributes(
self.graph, path_length_emerg, f"{weight}_to_{key}_{tag}"
)
def traffic_delay(self):
if "effective_travel_time" in self.edges:
distribution = (
self.edges["effective_travel_time"] / self.edges["travel_time"]
).to_numpy()
distribution = distribution[np.isfinite(distribution)]
weights = self.edges["length"].to_numpy()
weights = weights[np.isfinite(weights)]
numerator = np.sum(
[distribution[i] * weights[i] for i in range(len(distribution))]
)
denominator = np.sum(weights)
return numerator / denominator
else:
raise NameError("Use self.daganzo_methdod() first.")
def to_undirected(self):
self.graph = ox.utils_graph.get_undirected(self.graph)
update_nodes_edges(self)
def isochrones_from_shortest_path(
self,
centernodes,
time_arr=np.array([0, 2, 4, 8, 16]) * 60,
weight="travel_time",
alpha=100,
buffer=1e-3,
):
iso_gpd = im.isochrones_from_shortest_path(
self.graph, centernodes, time_arr, weight=weight, alpha=alpha
)
return iso_gpd
def plot(
self,
arg,
pad=0.2,
ax=False,
label=None,
vmin=False,
vmax=False,
center=None,
emergency_key=None,
emergency_tag=None,
xlim=None,
ylim=None,
cmap=mpl.colormaps["cividis"],
norm=None,
color_flooded_edges=False,
ascending=True,
legend=True,
):
if not vmin:
vmin = np.nanmin(arg[(arg != np.inf) & (arg != -np.inf)])
if not vmax:
vmax = np.nanmax(arg[(arg != np.inf) & (arg != -np.inf)])
if center == None:
center = (vmin + vmax) / 2
edges = self.edges.loc[arg.index]
if not ax:
fig, ax = plt.subplots()
if emergency_key != None:
emerg_nodes = self.nodes[self.nodes[emergency_key] == emergency_tag]
emerg_nodes.plot(marker="v", ax=ax, color="red", markersize=100, zorder=4)
cax = None
if legend:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=pad)
edges["cscheme"] = arg
if norm == None:
norm = colors.TwoSlopeNorm(vmin=vmin, vcenter=center, vmax=vmax)
sorted_edges = edges.sort_values(by="cscheme", key=abs, ascending=ascending)
prim = [
"motorway",
"trunk",
"primary",
"motorway_link",
"trunk_link",
"primary_link",
]
prim_edg = sorted_edges[sorted_edges["highway"].apply(lambda x: x in prim)]
sec = ["secondary", "secondary_link"]
sec_edg = sorted_edges[sorted_edges["highway"].apply(lambda x: x in sec)]
tert = ["tertiary", "tertiary_link"]
tert_edg = sorted_edges[sorted_edges["highway"].apply(lambda x: x in tert)]
if self.flooded:
f = gpd.GeoSeries(fm.flood_footprint())
f.plot(ax=ax, color="grey", alpha=1, zorder=5, aspect=None, hatch="\\\\")
if color_flooded_edges != False:
ft_edg = tert_edg[tert_edg["length"] == np.inf]
ft_edg.plot(ax=ax, edgecolor=color_flooded_edges, linewidth=2, zorder=4)
fs_edg = sec_edg[sec_edg["length"] == np.inf]
fs_edg.plot(ax=ax, edgecolor=color_flooded_edges, linewidth=4, zorder=4)
fp_edg = prim_edg[prim_edg["length"] == np.inf]
fp_edg.plot(ax=ax, edgecolor=color_flooded_edges, linewidth=7, zorder=4)
sorted_edges.plot(
column=arg,
ax=ax,
legend=legend,
vmin=vmin,
vmax=vmax,
cax=cax,
legend_kwds={"label": label},
linewidth=1,
cmap=cmap,
norm=norm,
)
tert_edg.plot(
column="cscheme",
ax=ax,
linewidth=2,
vmin=vmin,
vmax=vmax,
cmap=cmap,
norm=norm,
)
sec_edg.plot(
column="cscheme",
ax=ax,
linewidth=4,
vmin=vmin,
vmax=vmax,
cmap=cmap,
norm=norm,
)
prim_edg.plot(
column="cscheme",
ax=ax,
linewidth=7,
vmin=vmin,
vmax=vmax,
cmap=cmap,
norm=norm,
)
north, east, south, west = self.bbox
if xlim == None:
ax.set_xlim([west, east])
else:
ax.set_xlim(xlim)
if ylim == None:
ax.set_ylim([south, north])
else:
ax.set_ylim(ylim)
# return ax
def remove_edge(self, e):
graph = self.graph
u, v, k = e
graph[u][v][k]["speed_kph"] = 0
graph[u][v][k]["maxspeed"] = 0
graph[u][v][k]["travel_time"] = np.inf
graph[u][v][k]["length"] = np.inf
graph[u][v][k]["removed"] = True
update_nodes_edges(self)
def update_nodes_edges(self):
nodes, edges = ox.graph_to_gdfs(self.graph)
self.nodes = nodes
self.edges = edges
def remove_flooded_edges(self):
flood_poly = fm.flood_footprint()
nodes, edges = ox.graph_to_gdfs(self.graph)
flooded_edges = pd.Series(
{
e: fm.poly_great_circle_vec(flood_poly, e_poly) == 0
for e, e_poly in edges["geometry"].items()
},
name="flooded",
)
flooded_nodes = pd.Series(
{
n: fm.poly_great_circle_vec(flood_poly, n_poly) == 0
for n, n_poly in nodes["geometry"].items()
},
name="flooded",
)
nodes["flooded"] = flooded_nodes
edges["flooded"] = flooded_edges
nx.set_node_attributes(self.graph, nodes["flooded"], "flooded")
nx.set_edge_attributes(self.graph, edges["flooded"], "flooded")
edges.loc[flooded_edges, "speed_kph"] = 0
edges.loc[flooded_edges, "travel_time"] = np.inf
edges.loc[flooded_edges, "length"] = np.inf
nx.set_edge_attributes(self.graph, edges["speed_kph"], "speed_kph")
nx.set_edge_attributes(self.graph, edges["travel_time"], "travel_time")
nx.set_edge_attributes(self.graph, edges["length"], "length")
update_nodes_edges(self)
# if 'travel_time' in edges.columns:
# edges[flooded_edges]['travel_time'] = np.inf
# if 'effective_travel_time' in edges.column:
# edges[flooded_edges]['effective_travel_time'] = np.inf
#
# ox.graph_from_gdfs(nodes, edges)
#
# update_nodes_edges(self)
def graph_from_osmfile(north, east, south, west, date, roads):
"""Returns a graph of a region bounded by north, west, south, east.
Parameters
----------
date : Date YYYY-MM-DD of the osm-file
north : northern limit
east : eastern limit
south : southern limit
west : western limit
roads : Which road types are included in the graph. Can be 'primary', 'primary_secondary',
'primary_secondary_tertiary', 'all', or specified roads.
"""
if roads == "primary":
highw_names = "motorway,trunk,primary,motorway_link,trunk_link,primary_link"
elif roads == "primary_secondary":
highw_names = "motorway,trunk,primary,motorway_link,trunk_link,primary_link,secondary,secondary_link"
elif roads == "primary_secondary_tertiary":
highw_names = "motorway,trunk,primary,motorway_link,trunk_link,primary_link,secondary,secondary_link,tertiary,tertiary_link"
elif roads == "all":
highw_names = "motorway,trunk,primary,secondary,tertiary,motorway_link,trunk_link,primary_link,secondary_link,tertiary_link,unclassified,residential,living_street"
else:
highw_names = roads
if Path(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/graph-{highw_names}.gpickle"
).is_file():
with open(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/graph-{highw_names}.gpickle",
"rb",
) as pickle_file:
G = pickle.load(pickle_file)
# G = nx.read_gpickle(f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/graph-{highw_names}.gpickle")
print("Nodes:", len(G.nodes()))
print("Edges:", len(G.edges()))
return G
if not Path(f"data/cache/osmfiles/{date}").is_dir():
os.system(f"mkdir data/cache/osmfiles/{date}")
if not Path(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf"
).is_file():
os.system(f"mkdir data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}")
os.system(
f"osmium extract -b {west},{south},{east},{north} data/historic-data.nosync/{date}/germany.osm.pbf\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf --overwrite"
)
if not Path(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm"
).is_file():
os.system(
f"osmium tags-filter data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf\
w/highway={highw_names}\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.pbf --overwrite"
)
os.system(
f"osmium cat data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.pbf\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.bz2 --overwrite"
)
os.system(
f"bzip2 -d data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.bz2 -f"
)
G = ox.graph_from_xml(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm",
retain_all=False,
simplify=True,
)
# Assert geometry to all edges
nodes, edges = ox.graph_to_gdfs(G, fill_edge_geometry=True)
G = ox.graph_from_gdfs(nodes, edges, graph_attrs=G.graph)
for u, v, data in G.edges(keys=False, data=True):
assert "geometry" in data
# impute speed on all edges missing data
G = ox.add_edge_speeds(G)
# calculate travel time (seconds) for all edges
G = ox.add_edge_travel_times(G)
# add population estimate
G = pfr.population_from_raster_to_graph(G)
# R = residential_graph(date, 'residential,living_street', north, south, west, east)
# remove edges not in GGC
# G = ox.utils_graph.get_largest_component(G, strongly=True)
with open(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/graph-{highw_names}.gpickle",
"wb",
) as pickle_file:
pickle.dump(G, pickle_file)
print("Nodes:", len(G.nodes()))
print("Edges:", len(G.edges()))
return G
def residential_graph(date, highw_names, north, south, west, east):
"""Returns (unconnected) graph with only residential roads. Used for population estimate."""
if not Path(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm"
).is_file():
os.system(
f"osmium tags-filter data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/nofilter.osm.pbf\
w/highway={highw_names}\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.pbf --overwrite"
)
os.system(
f"osmium cat data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.pbf\
-o data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.bz2 --overwrite"
)
os.system(
f"bzip2 -d data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm.bz2 -f"
)
R = ox.graph_from_xml(
f"data/cache/osmfiles/{date}/{north}-{south}-{west}-{east}/filter-{highw_names}.osm",
retain_all=True,
simplify=True,
)
return R
# %%
def remove_edge(graph, edgs):
graph_copy = graph.copy()
for e in edgs:
u, v, k = e
graph_copy[u][v][k]["speed_kph"] = 0
graph_copy[u][v][k]["maxspeed"] = 0
graph_copy[u][v][k]["travel_time"] = np.inf
graph_copy[u][v][k]["length"] = np.inf
graph_copy[u][v][k]["removed"] = True
return graph_copy
def finite_sum(list):
return sum(list[np.isfinite(list)])
def finite_mean(list):
return np.mean(list[np.isfinite(list)])
def finite_max(list):
return max(list[np.isfinite(list)])
def finite_min(list):
return min(list[np.isfinite(list)])
def LODFs(graph, graph_r, kwd="load"):
edge = list(nx.get_edge_attributes(graph_r, "removed").keys())
if len(edge) == 0:
print("Graph does not contain any removed edges. Returning nothing.")
return
loads = np.array(list(nx.get_edge_attributes(graph, kwd).values()))
loads_r = np.array(list(nx.get_edge_attributes(graph_r, kwd).values()))
if len(loads) == 0:
print("Graph does not contain kwd. Returning nothing.")
return
loads_diff = loads_r - loads
lodfs = loads_diff / loads
nx.set_edge_attributes(graph_r, dict(zip(graph_r.edges, lodfs)), "LODF")
return graph_r
# %%
west, east = 7, 7.41
south, north = 50.0, 50.2
A = RoadNetwork(north, east, south, west, "2021-01-02", "primary_secondary_tertiary")
A.loads("travel_time")
load0 = nx.get_edge_attributes(A.graph, "load")
# %%
A.parallel_load("travel_time", cache=False)
load1 = nx.get_edge_attributes(A.graph, "load")
# %%
# %%
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment