Repository: martinfleis/momepy
Pull request #22 (2019-07-10 14:23:47)

ENH: node_density() based on nodes and length
Another (hopefully final) version of `node_density`. It no longer uses area, but length of reached network as a normalisation.

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

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

```diff
@@            Coverage Diff             @@
##           master      #22      +/-   ##
==========================================
+ Coverage   88.86%   89.53%   +0.66%     
==========================================
  Files           9        9              
  Lines        2327     2323       -4     
==========================================
+ Hits         2068     2080      +12     
+ Misses        259      243      -16
```


| [Impacted Files](https://codecov.io/gh/martinfleis/momepy/pull/22?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [momepy/intensity.py](https://codecov.io/gh/martinfleis/momepy/pull/22/diff?src=pr&el=tree#diff-bW9tZXB5L2ludGVuc2l0eS5weQ==) | `92.35% <90%> (+9.01%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/martinfleis/momepy/pull/22?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/22?src=pr&el=footer). Last update [6e3e729...c982185](https://codecov.io/gh/martinfleis/momepy/pull/22?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/docs/api.rst b/docs/api.rst
index 06306a3..6961d82 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -88,6 +88,7 @@ intensity
    elements_count
    floor_area_ratio
    gross_density
+   node_density
    reached
 
 connectivity
diff --git a/momepy/elements.py b/momepy/elements.py
index b160749..520b325 100644
--- a/momepy/elements.py
+++ b/momepy/elements.py
@@ -809,6 +809,9 @@ def get_network_id(elements, streets, unique_id, network_id, min_size=100):
 
     print('Generating centroids...')
     buildings_c = elements.copy()
+    if network_id in buildings_c.columns:
+        buildings_c = buildings_c.drop([network_id], axis=1)
+
     buildings_c['geometry'] = buildings_c.centroid  # make centroids
 
     print('Generating list of points...')
@@ -841,7 +844,11 @@ def get_network_id(elements, streets, unique_id, network_id, min_size=100):
     cleaned = joined[[unique_id, network_id]]
 
     print('Merging with elements...')
-    elements_m = elements.merge(cleaned, on=unique_id)
+    if network_id in elements.columns:
+        elements_copy = elements.copy().drop([network_id], axis=1)
+        elements_m = elements_copy.merge(cleaned, on=unique_id)
+    else:
+        elements_m = elements.merge(cleaned, on=unique_id)
 
     if elements_m[network_id].isnull().any():
         import warnings
diff --git a/momepy/intensity.py b/momepy/intensity.py
index ab0ccaf..e20bf22 100644
--- a/momepy/intensity.py
+++ b/momepy/intensity.py
@@ -474,27 +474,32 @@ def reached(streets, elements, unique_id, spatial_weights=None, mode='count', va
     return series
 
 
-def node_density(objects, nodes, spatial_weights=None, order=9, node_id='nodeID'):
+def node_density(left, right, spatial_weights, weighted=False, node_degree=None, node_start='node_start', node_end='node_end'):
     """
-    Calculate the density of nodes within topological steps of morphological tessellation.
-    Node is marked as reached if reached cell con
+    Calculate the density of nodes within topological steps defined in spatial_weights.
+
+    Calculated as number of nodes within k steps / cummulative length of street network within k steps.
+    node_start and node_end is standard output of :py:func:`momepy.nx_to_gdf` and is compulsory.
 
     .. math::
 
 
     Parameters
     ----------
-    objects : GeoDataFrame
-        GeoDataFrame containing tessellation objects to analyse
-    nodes : GeoDataFrame
-        GeoDataFrame containing nodes
+    left : GeoDataFrame
+        GeoDataFrame containing nodes of street network
+    right : GeoDataFrame
+        GeoDataFrame containing edges of street network
     spatial_weights : libpysal.weights, optional
-        spatial weights matrix - If None, Queen contiguity matrix of selected order will be calculated
-        based on objects.
-    order : int
-        order of Queen contiguity
-    node_id : str
-        name of the column of objects gdf with node id.
+        spatial weights matrix capturing relationship between nodes within set topological distance
+    weighted : bool
+        if True density will take into account node degree as k-1
+    node_degree : str
+        name of the column of left gdf containing node degree
+    node_start : str
+        name of the column of right gdf containing id of starting node
+    node_end : str
+        name of the column of right gdf containing id of ending node
 
     Returns
     -------
@@ -512,25 +517,27 @@ def node_density(objects, nodes, spatial_weights=None, order=9, node_id='nodeID'
     # define empty list for results
     results_list = []
 
-    print('Calculating gross density...')
-
-    if not all(objects.index == range(len(objects))):
-        raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
-
-    if spatial_weights is None:
-        print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
-        from momepy import Queen_higher
-        # matrix to define area of analysis (more steps)
-        spatial_weights = Queen_higher(k=order, geodataframe=objects)
+    print('Calculating node density...')
 
     # iterating over rows one by one
-    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
-        neighbours = spatial_weights.neighbors[index]
+    for index, row in tqdm(left.iterrows(), total=left.shape[0]):
+
+        neighbours = list(spatial_weights.neighbors[index])
         neighbours.append(index)
-        sub = objects.iloc[neighbours]
-        subnodes = set(sub['nodeID'])
-        results_list.append((len(subnodes) / sum(sub.geometry.area)) * 10000)
+        if weighted:
+            neighbour_nodes = left.iloc[neighbours]
+            number_nodes = sum(neighbour_nodes[node_degree] - 1)
+        else:
+            number_nodes = len(neighbours)
+
+        edg = right.loc[right['node_start'].isin(neighbours)].loc[right['node_end'].isin(neighbours)]
+        length = sum(edg.geometry.length)
+
+        if length > 0:
+            results_list.append(number_nodes / length)
+        else:
+            results_list.append(0)
 
     series = pd.Series(results_list)
-    print('Gross density calculated.')
+    print('Node density calculated.')
     return series


Test patch:
diff --git a/tests/test_elements.py b/tests/test_elements.py
index b908262..01221ea 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -34,6 +34,11 @@ class TestElements:
         buildings_id = mm.get_network_id(self.df_buildings, self.df_streets, 'uID', 'nID')
         assert not buildings_id.isna().any()
 
+    def test_get_network_id_duplicate(self):
+        self.df_buildings['nID'] = range(len(self.df_buildings))
+        buildings_id = mm.get_network_id(self.df_buildings, self.df_streets, 'uID', 'nID')
+        assert not buildings_id.isna().any()
+
     def test_get_node_id(self):
         nx = mm.gdf_to_nx(self.df_streets)
         nodes, edges = mm.nx_to_gdf(nx)
diff --git a/tests/test_intensity.py b/tests/test_intensity.py
index 06cc154..a09cf8c 100644
--- a/tests/test_intensity.py
+++ b/tests/test_intensity.py
@@ -82,3 +82,13 @@ class TestIntensity:
         assert max(area_v) == 79169.31385861784
         assert max(mean_v) == 7916.931385861784
         assert max(std_v) == 8995.18003493457
+
+    def test_node_density(self):
+        nx = mm.gdf_to_nx(self.df_streets)
+        nx = mm.node_degree(nx)
+        nodes, edges, W = mm.nx_to_gdf(nx, spatial_weights=True)
+        sw = mm.Queen_higher(k=3, weights=W)
+        density = mm.node_density(nodes, edges, sw)
+        weighted = mm.node_density(nodes, edges, sw, weighted=True, node_degree='degree')
+        assert density.mean() == 0.012690163074599968
+        assert weighted.mean() == 0.023207675994368446