Repository: martinfleis/momepy
Pull request #72 (2019-09-08 17:34:36)

ENH: Dual graph, minor graph changes
`gdf_to_nx` is now able to generate both primal and dual graph based on the street network. `nx_to_gdf` can turn both back to geodataframe, while for dual only lines are possible. 

+ some minor changes to graph algorithms

codecov[bot]: # [Codecov](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=h1) Report
> Merging [#72](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=desc) into [master](https://codecov.io/gh/martinfleis/momepy/commit/1f18b47e9318021121ed90bc425c8b81c469d237?src=pr&el=desc) will **increase** coverage by `0.2%`.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/martinfleis/momepy/pull/72/graphs/tree.svg?width=650&token=i04999UxkK&height=150&src=pr)](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##           master      #72     +/-   ##
=========================================
+ Coverage   93.12%   93.33%   +0.2%     
=========================================
  Files           9        9             
  Lines        2256     2311     +55     
=========================================
+ Hits         2101     2157     +56     
+ Misses        155      154      -1
```


| [Impacted Files](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [momepy/graph.py](https://codecov.io/gh/martinfleis/momepy/pull/72/diff?src=pr&el=tree#diff-bW9tZXB5L2dyYXBoLnB5) | `98.19% <100%> (+0.45%)` | :arrow_up: |
| [momepy/utils.py](https://codecov.io/gh/martinfleis/momepy/pull/72/diff?src=pr&el=tree#diff-bW9tZXB5L3V0aWxzLnB5) | `81.34% <100%> (+2.95%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=footer). Last update [1f18b47...90d5b1d](https://codecov.io/gh/martinfleis/momepy/pull/72?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/momepy/graph.py b/momepy/graph.py
index 00c51c6..d89bf1a 100644
--- a/momepy/graph.py
+++ b/momepy/graph.py
@@ -660,7 +660,7 @@ def _closeness_centrality(G, u=None, length=None, wf_improved=True, len_graph=No
 
 
 def local_closeness_centrality(
-    graph, radius=5, name="closeness", distance=None, length=None
+    graph, radius=5, name="closeness", distance=None, weight=None
 ):
     """
     Calculates local closeness for each node based on the defined distance.
@@ -684,8 +684,8 @@ def local_closeness_centrality(
     distance : str, optional
         Use specified edge data key as distance.
         For example, setting distance=’weight’ will use the edge weight to
-        measure the distance from the node n.
-    length : str, optional
+        measure the distance from the node n during ego_graph generation.
+    weight : str, optional
       Use the specified edge attribute as the edge distance in shortest
       path calculations in closeness centrality algorithm
 
@@ -711,13 +711,13 @@ def local_closeness_centrality(
             netx, n, radius=radius, undirected=True, distance=distance
         )  # define subgraph of steps=radius
         netx.nodes[n][name] = _closeness_centrality(
-            sub, n, length=length, len_graph=lengraph
+            sub, n, length=weight, len_graph=lengraph
         )
 
     return netx
 
 
-def global_closeness_centrality(graph, name="closeness", length="mm_len", **kwargs):
+def global_closeness_centrality(graph, name="closeness", weight="mm_len", **kwargs):
     """
     Calculates the closeness centrality for nodes.
 
@@ -733,8 +733,8 @@ def global_closeness_centrality(graph, name="closeness", length="mm_len", **kwar
         Ideally genereated from GeoDataFrame using :py:func:`momepy.gdf_to_nx`
     name : str, optional
         calculated attribute name
-    length : str (default 'mm_len')
-        attribute holding length of edge
+    weight : str (default 'mm_len')
+        attribute holding the weight of edge (e.g. length, angle)
     **kwargs
         kwargs for ``networkx.closeness_centrality``
 
@@ -749,14 +749,14 @@ def global_closeness_centrality(graph, name="closeness", length="mm_len", **kwar
     """
     netx = graph.copy()
 
-    vals = nx.closeness_centrality(netx, distance=length, **kwargs)
+    vals = nx.closeness_centrality(netx, distance=weight, **kwargs)
     nx.set_node_attributes(netx, vals, name)
 
     return netx
 
 
 def betweenness_centrality(
-    graph, name="betweenness", mode="nodes", length="mm_len", **kwargs
+    graph, name="betweenness", mode="nodes", weight="mm_len", endpoints=True, **kwargs
 ):
     """
     Calculates the shortest-path betweenness centrality for nodes.
@@ -773,8 +773,8 @@ def betweenness_centrality(
         calculated attribute name
     mode : str, default 'nodes'
         mode of betweenness calculation. 'node' for node-based, 'edges' for edge-based
-    length : str (default 'mm_len')
-        attribute holding length of edge
+    weight : str (default 'mm_len')
+        attribute holding the weight of edge (e.g. length, angle)
     **kwargs
         kwargs for ``networkx.betweenness_centrality`` or ``networkx.edge_betweenness_centrality``
 
@@ -798,16 +798,18 @@ def betweenness_centrality(
     G = nx.Graph()
     for u, v, k, data in netx.edges(data=True, keys=True):
         if G.has_edge(u, v):
-            if G[u][v][length] > netx[u][v][k][length]:
+            if G[u][v][weight] > netx[u][v][k][weight]:
                 nx.set_edge_attributes(G, {(u, v): data})
         else:
             G.add_edge(u, v, **data)
 
     if mode == "nodes":
-        vals = nx.betweenness_centrality(G, **kwargs)
+        vals = nx.betweenness_centrality(
+            G, weight=weight, endpoints=endpoints, **kwargs
+        )
         nx.set_node_attributes(netx, vals, name)
     elif mode == "edges":
-        vals = nx.edge_betweenness_centrality(G, **kwargs)
+        vals = nx.edge_betweenness_centrality(G, weight=weight, **kwargs)
         for u, v, k in netx.edges(keys=True):
             try:
                 val = vals[u, v]
@@ -827,7 +829,7 @@ def _euclidean(n, m):
     return math.sqrt((n[0] - m[0]) ** 2 + (n[1] - m[1]) ** 2)
 
 
-def _straightness_centrality(G, distance, normalized=True):
+def _straightness_centrality(G, weight, normalized=True):
     """
     Calculates straightness centrality.
     """
@@ -835,7 +837,7 @@ def _straightness_centrality(G, distance, normalized=True):
 
     for n in tqdm(G.nodes(), total=G.number_of_nodes()):
         straightness = 0
-        sp = nx.single_source_dijkstra_path_length(G, n, weight=distance)
+        sp = nx.single_source_dijkstra_path_length(G, n, weight=weight)
 
         if len(sp) > 0 and len(G) > 1:
             for target in sp:
@@ -854,7 +856,7 @@ def _straightness_centrality(G, distance, normalized=True):
 
 
 def straightness_centrality(
-    graph, length="mm_len", normalized=True, name="straightness"
+    graph, weight="mm_len", normalized=True, name="straightness"
 ):
     """
     Calculates the straightness centrality for nodes.
@@ -867,7 +869,7 @@ def straightness_centrality(
     graph : networkx.Graph
         Graph representing street network.
         Ideally genereated from GeoDataFrame using :py:func:`momepy.gdf_to_nx`
-    length : str (default 'mm_len')
+    weight : str (default 'mm_len')
         attribute holding length of edge
     normalized : bool
         normalize to number of nodes-1 in connected part
@@ -890,7 +892,7 @@ def straightness_centrality(
     """
     netx = graph.copy()
 
-    vals = _straightness_centrality(netx, distance=length, normalized=normalized)
+    vals = _straightness_centrality(netx, weight=weight, normalized=normalized)
     nx.set_node_attributes(netx, vals, name)
 
     return netx
@@ -911,7 +913,7 @@ def subgraph(
     edge_node_ratio=True,
     gamma=True,
     local_closeness=True,
-    closeness_length=None,
+    closeness_weight=None,
 ):
     """
     Calculates all subgraph-based characters.
@@ -954,7 +956,7 @@ def subgraph(
         Calculate gamma index (True/False)
     local_closeness : bool, default True
         Calculate local closeness centrality (True/False)
-    closeness_length : str, optional
+    closeness_weight : str, optional
       Use the specified edge attribute as the edge distance in shortest
       path calculations in closeness centrality algorithm
 
@@ -1012,7 +1014,7 @@ def subgraph(
         if local_closeness:
             lengraph = len(netx)
             netx.nodes[n]["local_closeness"] = _closeness_centrality(
-                sub, n, length=closeness_length, len_graph=lengraph
+                sub, n, length=closeness_weight, len_graph=lengraph
             )
 
     return netx
diff --git a/momepy/utils.py b/momepy/utils.py
index de56107..33a81a6 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -97,7 +97,71 @@ def sw_high(k, gdf=None, weights=None, ids=None, contiguity="queen", silent=True
     return joined
 
 
-def gdf_to_nx(gdf_network, length="mm_len"):
+def _angle(a, b, c):
+    """
+    Measure angle between a-b, b-c. In radians.
+    Helper for gdf_to_nx.
+    """
+    ba = [aa - bb for aa, bb in zip(a, b)]
+    bc = [cc - bb for cc, bb in zip(c, b)]
+    nba = math.sqrt(sum((x ** 2.0 for x in ba)))
+    ba = [x / nba for x in ba]
+    nbc = math.sqrt(sum((x ** 2.0 for x in bc)))
+    bc = [x / nbc for x in bc]
+    scal = sum((aa * bb for aa, bb in zip(ba, bc)))
+    angle = math.acos(round(scal, 10))
+    return angle
+
+
+def _generate_primal(G, gdf_network, fields):
+    """
+    Generate primal graph.
+    Helper for gdf_to_nx.
+    """
+    G.graph["approach"] = "primal"
+    key = 0
+    for index, row in gdf_network.iterrows():
+        first = row.geometry.coords[0]
+        last = row.geometry.coords[-1]
+
+        data = [row[f] for f in fields]
+        attributes = dict(zip(fields, data))
+        G.add_edge(first, last, key=key, **attributes)
+        key += 1
+
+
+def _generate_dual(G, gdf_network, fields):
+    """
+    Generate dual graph
+    Helper for gdf_to_nx.
+    """
+    G.graph["approach"] = "dual"
+    sw = libpysal.weights.Queen.from_dataframe(gdf_network)
+    gdf_network["mm_cent"] = gdf_network.geometry.centroid
+
+    for index, row in gdf_network.iterrows():
+        centroid = (row.mm_cent.x, row.mm_cent.y)
+        data = [row[f] for f in fields]
+        attributes = dict(zip(fields, data))
+        G.add_node(centroid, **attributes)
+
+        if sw.cardinalities[index] > 0:
+            for n in sw.neighbors[index]:
+                start = centroid
+                end = list(gdf_network.loc[n, "mm_cent"].coords)[0]
+                p0 = row.geometry.coords[0]
+                p1 = row.geometry.coords[-1]
+                p2 = gdf_network.iloc[n]["geometry"].coords[0]
+                p3 = gdf_network.iloc[n]["geometry"].coords[-1]
+                points = [p0, p1, p2, p3]
+                shared = [x for x in points if points.count(x) > 1]
+                if shared:  # fix for non-planar graph
+                    remaining = [e for e in points if e not in [shared[0]]]
+                    angle = _angle(remaining[0], shared[0], remaining[1])
+                    G.add_edge(start, end, angle=angle)
+
+
+def gdf_to_nx(gdf_network, approach="primal", length="mm_len"):
     """
     Convert LineString GeoDataFrame to networkx.Graph
 
@@ -105,7 +169,9 @@ def gdf_to_nx(gdf_network, length="mm_len"):
     ----------
     gdf_network : GeoDataFrame
         GeoDataFrame containing objects to convert
-    length : str, optional
+    approach : str, default 'primal'
+        Decide wheter genereate 'primal' or 'dual' graph.
+    length : str, default mm_len
         name of attribute of segment length (geographical) which will be saved to graph
 
     Returns
@@ -114,26 +180,102 @@ def gdf_to_nx(gdf_network, length="mm_len"):
         Graph
 
     """
+    gdf_network = gdf_network.copy()
     # generate graph from GeoDataFrame of LineStrings
     net = nx.MultiGraph()
     net.graph["crs"] = gdf_network.crs
     gdf_network[length] = gdf_network.geometry.length
     fields = list(gdf_network.columns)
 
-    key = 0
-    for index, row in gdf_network.iterrows():
-        first = row.geometry.coords[0]
-        last = row.geometry.coords[-1]
+    if approach == "primal":
+        _generate_primal(net, gdf_network, fields)
 
-        data = [row[f] for f in fields]
-        attributes = dict(zip(fields, data))
-        net.add_edge(first, last, key=key, **attributes)
-        key += 1
+    elif approach == "dual":
+        _generate_dual(net, gdf_network, fields)
+
+    else:
+        raise ValueError(
+            "Approach {} is not supported. Use 'primal' or 'dual'.".format(approach)
+        )
 
     return net
 
 
-def nx_to_gdf(net, nodes=True, edges=True, spatial_weights=False, nodeID="nodeID"):
+def _points_to_gdf(net, spatial_weights, nodeID):
+    """
+    Generate point gdf from nodes.
+    Helper for nx_to_gdf.
+    """
+    nid = 1
+    for n in net:
+        net.nodes[n][nodeID] = nid
+        nid += 1
+    node_xy, node_data = zip(*net.nodes(data=True))
+    gdf_nodes = gpd.GeoDataFrame(
+        list(node_data), geometry=[Point(i, j) for i, j in node_xy]
+    )
+    gdf_nodes.crs = net.graph["crs"]
+    return gdf_nodes
+
+
+def _lines_to_gdf(net, lines, points, nodeID):
+    """
+    Generate linestring gdf from edges.
+    Helper for nx_to_gdf.
+    """
+    starts, ends, edge_data = zip(*net.edges(data=True))
+    if lines is True:
+        node_start = []
+        node_end = []
+        for s in starts:
+            node_start.append(net.node[s][nodeID])
+        for e in ends:
+            node_end.append(net.node[e][nodeID])
+    gdf_edges = gpd.GeoDataFrame(list(edge_data))
+    if points is True:
+        gdf_edges["node_start"] = node_start
+        gdf_edges["node_end"] = node_end
+    gdf_edges.crs = net.graph["crs"]
+    return gdf_edges
+
+
+def _primal_to_gdf(net, points, lines, spatial_weights, nodeID):
+    """
+    Generate gdf(s) from primal network.
+    Helper for nx_to_gdf.
+    """
+    if points is True:
+        gdf_nodes = _points_to_gdf(net, spatial_weights, nodeID)
+
+        if spatial_weights is True:
+            W = libpysal.weights.W.from_networkx(net)
+
+    if lines is True:
+        gdf_edges = _lines_to_gdf(net, lines, points, nodeID)
+
+    if points is True and lines is True:
+        if spatial_weights is True:
+            return gdf_nodes, gdf_edges, W
+        return gdf_nodes, gdf_edges
+    if points is True and lines is False:
+        if spatial_weights is True:
+            return gdf_nodes, W
+        return gdf_nodes
+    return gdf_edges
+
+
+def _dual_to_gdf(net):
+    """
+    Generate linestring gdf from dual network.
+    Helper for nx_to_gdf.
+    """
+    starts, edge_data = zip(*net.nodes(data=True))
+    gdf_edges = gpd.GeoDataFrame(list(edge_data))
+    gdf_edges.crs = net.graph["crs"]
+    return gdf_edges
+
+
+def nx_to_gdf(net, points=True, lines=True, spatial_weights=False, nodeID="nodeID"):
     """
     Convert networkx.Graph to LineString GeoDataFrame and Point GeoDataFrame
 
@@ -157,43 +299,27 @@ def nx_to_gdf(net, nodes=True, edges=True, spatial_weights=False, nodeID="nodeID
 
     """
     # generate nodes and edges geodataframes from graph
-    if nodes is True:
-        nid = 1
-        for n in net:
-            net.nodes[n][nodeID] = nid
-            nid += 1
-        node_xy, node_data = zip(*net.nodes(data=True))
-        gdf_nodes = gpd.GeoDataFrame(
-            list(node_data), geometry=[Point(i, j) for i, j in node_xy]
+    try:
+        if net.graph["approach"] == "primal":
+            return _primal_to_gdf(
+                net,
+                points=points,
+                lines=lines,
+                spatial_weights=spatial_weights,
+                nodeID=nodeID,
+            )
+        if net.graph["approach"] == "dual":
+            return _dual_to_gdf(net)
+        raise ValueError(
+            "Approach {} is not supported. Use 'primal' or 'dual'.".format(
+                net.graph["approach"]
+            )
+        )
+    except KeyError:
+        raise KeyError(
+            "Approach is not set. Graph needs an attribute 'approach' to be set"
+            "either to 'primal' or to 'dual'.".format(net.graph["approach"])
         )
-        gdf_nodes.crs = net.graph["crs"]
-        if spatial_weights is True:
-            W = libpysal.weights.W.from_networkx(net)
-
-    if edges is True:
-        starts, ends, edge_data = zip(*net.edges(data=True))
-        if nodes is True:
-            node_start = []
-            node_end = []
-            for s in starts:
-                node_start.append(net.node[s][nodeID])
-            for e in ends:
-                node_end.append(net.node[e][nodeID])
-        gdf_edges = gpd.GeoDataFrame(list(edge_data))
-        if nodes is True:
-            gdf_edges["node_start"] = node_start
-            gdf_edges["node_end"] = node_end
-        gdf_edges.crs = net.graph["crs"]
-
-    if nodes is True and edges is True:
-        if spatial_weights is True:
-            return gdf_nodes, gdf_edges, W
-        return gdf_nodes, gdf_edges
-    if nodes is True and edges is False:
-        if spatial_weights is True:
-            return gdf_nodes, W
-        return gdf_nodes
-    return gdf_edges
 
 
 def limit_range(vals, rng):


Test patch:
diff --git a/tests/test_graph.py b/tests/test_graph.py
index bb45e20..0fd87a9 100644
--- a/tests/test_graph.py
+++ b/tests/test_graph.py
@@ -10,6 +10,7 @@ class TestGraph:
         self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.network = mm.gdf_to_nx(self.df_streets)
         self.network = mm.node_degree(self.network)
+        self.dual = mm.gdf_to_nx(self.df_streets, approach="dual")
 
     def test_node_degree(self):
         deg1 = mm.node_degree(graph=self.network)
@@ -71,26 +72,28 @@ class TestGraph:
         check = 0.4444444444444444
         assert net.nodes[(1603650.450422848, 6464368.600601688)]["gamma"] == check
 
-    def test_local_closeness(self):
+    def test_local_closeness_centrality(self):
         net = mm.local_closeness_centrality(self.network)
         check = 0.27557319223985893
         assert net.nodes[(1603650.450422848, 6464368.600601688)]["closeness"] == check
-        net2 = mm.local_closeness_centrality(self.network, length="mm_len")
+        net2 = mm.local_closeness_centrality(self.network, weight="mm_len")
         check2 = 0.0015544070362478774
         assert net2.nodes[(1603650.450422848, 6464368.600601688)]["closeness"] == check2
 
     def test_global_closeness_centrality(self):
-        net = mm.global_closeness_centrality(self.network, length="mm_len")
+        net = mm.global_closeness_centrality(self.network, weight="mm_len")
         check = 0.0016066095164175716
         assert net.nodes[(1603650.450422848, 6464368.600601688)]["closeness"] == check
 
     def test_betweenness_centrality(self):
         net = mm.betweenness_centrality(self.network)
         net2 = mm.betweenness_centrality(self.network, mode="edges")
+        angular = mm.betweenness_centrality(self.dual, weight="angle")
         with pytest.raises(ValueError):
             mm.betweenness_centrality(self.network, mode="nonexistent")
-        node = 0.15806878306878305
-        edge = 0.20320197044334976
+        node = 0.2413793103448276
+        edge = 0.16995073891625617
+        ang_b = 0.16470588235294117
         assert net.nodes[(1603650.450422848, 6464368.600601688)]["betweenness"] == node
         assert (
             net2.edges[
@@ -100,6 +103,10 @@ class TestGraph:
             ]["betweenness"]
             == edge
         )
+        assert (
+            angular.nodes[(1603315.3564306537, 6464044.376339891)]["betweenness"]
+            == ang_b
+        )
 
     def test_straightness_centrality(self):
         net = mm.straightness_centrality(self.network)
@@ -133,7 +140,7 @@ class TestGraph:
 
     def test_subgraph(self):
         net = mm.subgraph(self.network)
-        nodes = mm.nx_to_gdf(net, edges=False)
+        nodes = mm.nx_to_gdf(net, lines=False)
         cols = [
             "meshedness",
             "cds_length",
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 5f67e7c..eeb699e 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -2,6 +2,7 @@ import momepy as mm
 import geopandas as gpd
 import numpy as np
 import libpysal
+import networkx
 
 import pytest
 
@@ -39,6 +40,11 @@ class TestUtils:
         nx = mm.gdf_to_nx(self.df_streets)
         assert nx.number_of_nodes() == 29
         assert nx.number_of_edges() == 35
+        dual = mm.gdf_to_nx(self.df_streets, approach="dual")
+        assert dual.number_of_nodes() == 35
+        assert dual.number_of_edges() == 148
+        with pytest.raises(ValueError):
+            mm.gdf_to_nx(self.df_streets, approach="nonexistent")
 
     def test_nx_to_gdf(self):
         nx = mm.gdf_to_nx(self.df_streets)
@@ -49,13 +55,22 @@ class TestUtils:
         nodes, edges = mm.nx_to_gdf(nx)
         assert len(nodes) == 29
         assert len(edges) == 35
-        edges = mm.nx_to_gdf(nx, nodes=False)
+        edges = mm.nx_to_gdf(nx, points=False)
         assert len(edges) == 35
-        nodes, W = mm.nx_to_gdf(nx, edges=False, spatial_weights=True)
+        nodes, W = mm.nx_to_gdf(nx, lines=False, spatial_weights=True)
         assert len(nodes) == 29
         assert W.n == 29
-        nodes = mm.nx_to_gdf(nx, edges=False, spatial_weights=False)
+        nodes = mm.nx_to_gdf(nx, lines=False, spatial_weights=False)
         assert len(nodes) == 29
+        dual = mm.gdf_to_nx(self.df_streets, approach="dual")
+        edges = mm.nx_to_gdf(dual)
+        assert len(edges) == 35
+        dual.graph["approach"] = "nonexistent"
+        with pytest.raises(ValueError):
+            mm.nx_to_gdf(dual)
+        G = networkx.Graph()
+        with pytest.raises(KeyError):
+            mm.nx_to_gdf(G)
 
     def test_limit_range(self):
         assert mm.limit_range(range(10), rng=(25, 75)) == [2, 3, 4, 5, 6, 7]