import os
import numpy as np
import geopandas as gpd
from ..utils.geo import get_subgraph
import shapely
from shapely.geometry import Point, LineString
import networkx as nx
import rasterio as rio
import fiona
import pickle
from multiprocessing import Pool
class Node(object):
"""An object to hold node attributes.
Attributes
----------
idx : int
The numerical index of the node. Used as a unique identifier
when the nodes are added to the graph.
x : `int` or `float`
Numeric x location of the node, in either a geographic CRS or in pixel
coordinates.
y : `int` or `float`
Numeric y location of the node, in either a geographic CRS or in pixel
coordinates.
"""
def __init__(self, idx, x, y):
self.idx = idx
self.x = x
self.y = y
def __repr__(self):
return 'Node {} at ({}, {})'.format(self.idx, self.x, self.y)
class Edge(object):
"""An object to hold edge attributes.
Attributes
----------
nodes : 2-`tuple` of :class:`Node` s
:class:`Node` instances connected by the edge.
weight : int or float
The weight of the edge.
"""
def __init__(self, nodes, edge_weight=None):
self.nodes = nodes
self.weight = edge_weight
def __repr__(self):
return 'Edge between {} and {} with weight {}'.format(self.nodes[0],
self.nodes[1],
self.weight)
def set_edge_weight(self, normalize_factor=None, inverse=False):
"""Get the edge weight based on Euclidean distance between nodes.
Note
----
This method does not account for spherical deformation (i.e. does not
use the Haversine equation). It is a simple linear distance.
Arguments
---------
normalize_factor : `int` or `float`, optional
a number to multiply (or divide, if
``inverse=True``) the Euclidean distance by. Defaults to ``None``
(no normalization)
inverse : bool, optional
if ``True``, the Euclidean distance weight will be divided by
``normalize_factor`` instead of multiplied by it.
"""
weight = np.linalg.norm(
np.array((self.nodes[0].x, self.nodes[0].y)) -
np.array((self.nodes[1].x, self.nodes[1].y)))
if normalize_factor is not None:
if inverse:
weight = weight/normalize_factor
else:
weight = weight*normalize_factor
self.weight = weight
def get_node_idxs(self):
"""Return the Node.idx for the nodes in the edge."""
return (self.nodes[0].idx, self.nodes[1].idx)
class Path(object):
"""An object to hold :class:`Edge` s with common properties.
Attributes
----------
edges : `list` of :class:`Edge` s
A `list` of :class:`Edge` s
properties : dict
A dictionary of property: value pairs that provide relevant metadata
about edges along the path (e.g. road type, speed limit, etc.)
"""
def __init__(self, edges=None, properties=None):
self.edges = edges
if properties is None:
properties = {}
self.properties = properties
def __repr__(self):
return 'Path including {}'.format([e for e in self.edges])
def add_edge(self, edge):
"""Add an edge to the path."""
self.edges.append(edge)
def set_edge_weights(self, data_key=None, inverse=False, overwrite=True):
"""Calculate edge weights for all edges in the Path."""
for edge in self.edges:
if not overwrite and edge.weight is not None:
continue
if data_key is not None:
edge.set_edge_weight(
normalize_factor=self.properties[data_key],
inverse=inverse)
else:
edge.set_edge_weight()
def add_data(self, property, value):
"""Add a property: value pair to the Path.properties attribute."""
self.properties[property] = value
def __iter__(self):
"""Iterate through edges in the path."""
yield from self.edges
def geojson_to_graph(geojson, graph_name=None, retain_all=True,
valid_road_types=None, road_type_field='type', edge_idx=0,
first_node_idx=0, weight_norm_field=None, inverse=False,
workers=1, verbose=False, output_path=None):
"""Convert a geojson of path strings to a network graph.
Arguments
---------
geojson : str
Path to a geojson file (or any other OGR-compatible vector file) to
load network edges and nodes from.
graph_name : str, optional
Name of the graph. If not provided, graph will be named ``'unnamed'`` .
retain_all : bool, optional
If ``True`` , the entire graph will be returned even if some parts are
not connected. Defaults to ``True``.
valid_road_types : :class:`list` of :class:`int` s, optional
The road types to permit in the graph. If not provided, it's assumed
that all road types are permitted. The possible values are integers
``1``-``7``, which map as follows::
1: Motorway
2: Primary
3: Secondary
4: Tertiary
5: Residential
6: Unclassified
7: Cart track
road_type_field : str, optional
The name of the property in the vector data that delineates road type.
Defaults to ``'type'`` .
edge_idx : int, optional
The first index to use for an edge. This can be set to a higher value
so that a graph's edge indices don't overlap with existing values in
another graph.
first_node_idx : int, optional
The first index to use for a node. This can be set to a higher value
so that a graph's node indices don't overlap with existing values in
another graph.
weight_norm_field : str, optional
The name of a field in `geojson` to pass to argument ``data_key`` in
:func:`Path.set_edge_weights`. Defaults to ``None``, in which case
no weighting is performed (weights calculated solely using Euclidean
distance.)
workers : int, optional
Number of parallel processes to run for parallelization. Defaults to 1.
Should not be greater than the number of CPUs available.
verbose : bool, optional
Verbose print output. Defaults to ``False`` .
output_path : str, optional
Path to a pickle file to save the output graph to. Nothing will be
saved to disk if not provided.
Returns
-------
G : :class:`networkx.MultiDiGraph`
A :class:`networkx.MultiDiGraph` containing all of the nodes and edges
from the geojson (or only the largest connected component if
`retain_all` = ``False``). Edge lengths are weighted based on
geographic distance.
"""
with fiona.open(geojson, 'r') as f:
crs = f.crs
f.close()
# due to an annoying feature of loading these graphs, the numeric road
# type identifiers are presented as string versions. we therefore reformat
# the valid_road_types list as strings.
if valid_road_types is not None:
valid_road_types = [str(i) for i in valid_road_types]
# create the graph as a MultiGraph and set the original CRS to EPSG 4326
# extract nodes and paths
nodes, paths = get_nodes_paths(geojson,
valid_road_types=valid_road_types,
first_node_idx=first_node_idx,
road_type_field=road_type_field,
workers=workers, verbose=verbose)
# nodes is a dict of node_idx: node_params (e.g. location, metadata)
# pairs.
# paths is a dict of path dicts. the path key is the path_idx.
# each path dict has a list of node_idxs as well as properties metadata.
# initialize the graph object
G = nx.MultiDiGraph(name=graph_name, crs=crs)
if not nodes: # if there are no nodes in the graph
return G
if verbose:
print("nodes:", nodes)
print("paths:", paths)
# add each osm node to the graph
for node in nodes:
G.add_node(node.idx, **{'x': node.x, 'y': node.y})
# add each path to the graph
for path in paths:
# calculate edge length using euclidean distance and a weighting term
path.set_edge_weights(data_key=weight_norm_field, inverse=inverse)
edges = [(*[node.idx for node in edge.nodes],
edge.weight) for edge in path]
if verbose:
print(edges)
G.add_weighted_edges_from(edges)
if not retain_all:
# keep only largest connected component of graph unless retain_all
# code modified from osmnx.core.get_largest_component & induce_subgraph
largest_cc = max(nx.weakly_connected_components(G), key=len)
G = get_subgraph(G, largest_cc)
if output_path:
with open(output_path, 'wb') as f:
pickle.dump(G, f)
f.close()
return G
def get_nodes_paths(vector_file, first_node_idx=0, node_gdf=gpd.GeoDataFrame(),
valid_road_types=None, road_type_field='type', workers=1,
verbose=False):
"""
Extract nodes and paths from a vector file.
Arguments
---------
vector_file : str
Path to an OGR-compatible vector file containing line segments (e.g.,
JSON response from from the Overpass API, or a SpaceNet GeoJSON).
first_path_idx : int, optional
The first index to use for a path. This can be set to a higher value
so that a graph's path indices don't overlap with existing values in
another graph.
first_node_idx : int, optional
The first index to use for a node. This can be set to a higher value
so that a graph's node indices don't overlap with existing values in
another graph.
node_gdf : :class:`geopandas.GeoDataFrame` , optional
A :class:`geopandas.GeoDataFrame` containing nodes to add to the graph.
New nodes will be added to this object incrementally during the
function call.
valid_road_types : :class:`list` of :class:`int` s, optional
The road types to permit in the graph. If not provided, it's assumed
that all road types are permitted. The possible values are integers
``1``-``7``, which map as follows::
1: Motorway
2: Primary
3: Secondary
4: Tertiary
5: Residential
6: Unclassified
7: Cart track
road_type_field : str, optional
The name of the attribute containing road type information in
`vector_file`. Defaults to ``'type'``.
workers : int, optional
Number of worker processes to use for parallelization. Defaults to 1.
Should not exceed the number of CPUs available.
verbose : bool, optional
Verbose print output. Defaults to ``False``.
Returns
-------
nodes, paths : `tuple` of `dict` s
nodes : list
A `list` of :class:`Node` s to be added to the graph.
paths : list
A list of :class:`Path` s containing the :class:`Edge` s and
:class:`Node` s to be added to the graph.
"""
if valid_road_types is None:
valid_road_types = ['1', '2', '3', '4', '5', '6', '7']
with fiona.open(vector_file, 'r') as source:
with Pool(processes=workers) as pool:
node_list = pool.map(_get_all_nodes, source,
chunksize=10)
pool.close()
source.close()
# convert to geoseries and drop duplicates (have to flatten first)
node_series = gpd.GeoSeries([i for sublist in node_list for i in sublist])
# NOTE: It is ESSENTIAL to use keep='last' in the line below; otherwise, it
# misses a duplicate if it includes the first element of the series.
node_series = node_series.drop_duplicates(keep='last')
node_series = node_series.reset_index(drop=True)
node_series.name = 'geometry'
node_series.index.name = 'node_idx'
node_gdf = gpd.GeoDataFrame(node_series.reset_index())
node_gdf['node'] = node_gdf.apply(
lambda p: Node(p['node_idx'], p['geometry'].x, p['geometry'].y),
axis=1)
# create another parallelized operation to iterate through edges
# _init_worker passes the node_series to every process in the pool
with fiona.open(vector_file, 'r') as source:
with Pool(
processes=workers, initializer=_init_worker,
initargs=(node_gdf, valid_road_types, road_type_field)
) as pool:
zipped_edges_properties = pool.map(parallel_linestring_to_path,
source, chunksize=10)
pool.close()
source.close()
nodes = node_gdf['node'].tolist()
paths = []
# it would've been better to do this within the multiprocessing pool but
# it's REALLY hard to share objects in memory across processes without
# copies being made (and therefore nodes being duplicated)
for edges, properties in zipped_edges_properties:
path = Path(
edges=[Edge((nodes[edge[0]], nodes[edge[1]])) for edge in edges],
properties=properties
)
paths.append(path)
return nodes, paths
def parallel_linestring_to_path(feature):
"""Read in a feature line from a fiona-opened shapefile and get the edges.
Arguments
---------
feature : dict
An item from a :class:`fiona.open` iterable with the key ``'geometry'``
containing :class:`shapely.geometry.line.LineString` s or
:class:`shapely.geometry.line.MultiLineString` s.
Returns
-------
A list of :class:`Path` s containing all edges in the LineString or
MultiLineString.
Notes
-----
This function depends on ``node_series`` and ``valid_road_types``, which
are passed by an initializer as items in ``var_dict``.
"""
properties = feature['properties']
# TODO: create more adjustable filter
if var_dict['road_type_field'] in properties:
road_type = properties[var_dict['road_type_field']]
elif 'highway' in properties:
road_type = properties['highway']
elif 'road_type' in properties:
road_type = properties['road_type']
else:
road_type = 'None'
geom = feature['geometry']
if geom['type'] == 'LineString' or \
geom['type'] == 'MultiLineString':
if road_type not in var_dict['valid_road_types'] or \
'LINESTRING EMPTY' in properties.values():
return
if geom['type'] == 'LineString':
linestring = shapely.geometry.shape(geom)
edges = linestring_to_edges(linestring, var_dict['node_gdf'])
elif geom['type'] == 'MultiLineString':
# do the same thing as above, but do it for each piece
edges = []
for linestring in shapely.geometry.shape(geom):
edge_set, node_idx, node_gdf = linestring_to_edges(
linestring, var_dict['node_gdf'])
edges.extend(edge_set)
return edges, properties
def linestring_to_edges(linestring, node_gdf):
"""Collect nodes in a linestring and add them to an edge.
Arguments
---------
linestring : :class:`shapely.geometry.LineString`
A :class:`shapely.geometry.LineString` object to extract nodes and
edges from.
node_series : :class:`geopandas.GeoSeries`
A :class:`geopandas.GeoSeries` containing a
:class:`shapely.geometry.point.Point` for every node to be added to the
graph.
Returns
-------
edges : list
A list of :class:`Edge` s from ``linestring``.
"""
edges = []
nodes = []
for point in linestring.coords:
point_shp = shapely.geometry.shape(Point(point))
nodes.append(
node_gdf.node_idx[node_gdf.distance(point_shp) == 0.0].values[0]
)
if len(nodes) > 1:
edges.append(nodes[-2:])
return edges
def graph_to_geojson(G, output_path, encoding='utf-8', overwrite=False,
verbose=False):
"""
Save graph to two geojsons: one containing nodes, the other edges.
Arguments
---------
G : :class:`networkx.MultiDiGraph`
A graph object to save to geojson files.
output_path : str
Path to save the geojsons to. ``'_nodes.geojson'`` and
``'_edges.geojson'`` will be appended to ``output_path`` (after
stripping the extension).
encoding : str, optional
The character encoding for the saved files.
overwrite : bool, optional
Should files at ``output_path`` be overwritten? Defaults to no
(``False``).
verbose : bool, optional
Switch to print relevant values. Defaults to no (``False``).
Notes
-----
This function is based on ``osmnx.save_load.save_graph_shapefile``, with
tweaks to make it work with our graph objects. It will save two geojsons:
a file containing all of the nodes and a file containing all of the edges.
When writing to geojson, must convert the coordinate reference system
(crs) to string if it's a dict, otherwise no crs will be appended to the
geojson.
Returns
-------
None
"""
# convert directed graph G to an undirected graph for saving as a shapefile
G_to_save = G.copy().to_undirected()
# create GeoDataFrame containing all of the nodes
nodes, data = zip(*G_to_save.nodes(data=True))
gdf_nodes = gpd.GeoDataFrame(list(data), index=nodes)
# get coordinate reference system
g_crs = G_to_save.graph['crs']
if type(g_crs) == dict:
# convert from dict
g_crs = rio.crs.CRS.from_dict(g_crs)
gdf_nodes.crs = g_crs
if verbose:
print("crs:", g_crs)
gdf_nodes['geometry'] = gdf_nodes.apply(
lambda row: Point(row['x'], row['y']), axis=1
)
gdf_nodes = gdf_nodes.drop(['x', 'y'], axis=1)
# gdf_nodes['node_idx'] = gdf_nodes['node_idx'].astype(np.int32)
# # make everything but geometry column a string
# for col in [c for c in gdf_nodes.columns if not c == 'geometry']:
# gdf_nodes[col] = gdf_nodes[col].fillna('').map(make_str)
# create GeoDataFrame containing all of the edges
edges = []
for u, v, key, data in G_to_save.edges(keys=True, data=True):
edge = {'key': key}
for attr_key in data:
edge[attr_key] = data[attr_key]
if 'geometry' not in data:
point_u = Point((G_to_save.nodes[u]['x'], G_to_save.nodes[u]['y']))
point_v = Point((G_to_save.nodes[v]['x'], G_to_save.nodes[v]['y']))
edge['geometry'] = LineString([point_u, point_v])
edges.append(edge)
gdf_edges = gpd.GeoDataFrame(edges)
gdf_edges.crs = g_crs
for col in [c for c in gdf_nodes.columns if c != 'geometry']:
gdf_nodes[col] = gdf_nodes[col].fillna('').apply(str)
for col in [c for c in gdf_edges.columns if c != 'geometry']:
gdf_edges[col] = gdf_edges[col].fillna('').apply(str)
# make directory structure
if not os.path.exists(os.path.split(output_path)[0]):
os.makedirs(os.path.split(output_path)[0])
edges_path = os.path.splitext(output_path)[0] + '_edges.geojson'
nodes_path = os.path.splitext(output_path)[0] + '_nodes.geojson'
if overwrite:
if os.path.exists(edges_path):
os.remove(edges_path)
if os.path.exists(nodes_path):
os.remove(nodes_path)
gdf_edges.to_file(edges_path, encoding=encoding, driver='GeoJSON')
gdf_nodes.to_file(nodes_path, encoding=encoding, driver='GeoJSON')
def _get_all_nodes(feature):
"""Create a list of node geometries from a geojson of (multi)linestrings.
Note
----
This function is intended to be used with pool.imap_unordered for
parallelization.
Returns
-------
A list of :class:`shapely.geometry.Point` instances. DUPLICATES CAN EXIST.
"""
points = []
geom = feature['geometry']
if geom['type'] == 'LineString':
linestring = shapely.geometry.shape(geom)
points.extend(_get_linestring_points(linestring))
elif geom['type'] == 'MultiLineString':
for linestring in shapely.geometry.shape(geom):
points.extend(_get_linestring_points(linestring))
return points
def _get_linestring_points(linestring):
points = []
for point in linestring.coords:
points.append(shapely.geometry.shape(Point(point)))
return points
def _init_worker(node_gdf, valid_road_types, road_type_field):
the_dict = {'node_gdf': node_gdf,
'valid_road_types': valid_road_types,
'road_type_field': road_type_field}
global var_dict
var_dict = the_dict