Repository: martinfleis/momepy
Pull request #5 (2019-04-25 15:38:57)

ENH: Effective mesh modes
Allow interquartile and interdecile modes for effective mesh. During that, I've found out issues with spatial weights, which will be fixed soon (hopefully).


Patch:
diff --git a/momepy/dimension.py b/momepy/dimension.py
index 92c4dfb..a9ad480 100644
--- a/momepy/dimension.py
+++ b/momepy/dimension.py
@@ -299,7 +299,7 @@ def longest_axis_length(objects):
     return series
 
 
-def effective_mesh(objects, spatial_weights=None, areas=None, order=3):
+def effective_mesh(objects, spatial_weights=None, areas=None, order=3, mode=None):
     """
     Calculates the effective mesh size of morphological tessellation
 
@@ -320,6 +320,8 @@ def effective_mesh(objects, spatial_weights=None, areas=None, order=3):
         during the process without saving them separately.
     order : int
         order of Queen contiguity
+    mode : str
+        mode of calculation. None will use all values, `iq` will use values within interquartile range, `id` will use values within interdecile range
 
     Returns
     -------
@@ -364,6 +366,9 @@ def effective_mesh(objects, spatial_weights=None, areas=None, order=3):
         from momepy import Queen_higher
         # matrix to define area of analysis (more steps)
         spatial_weights = Queen_higher(k=order, geodataframe=objects)
+    else:
+        if not all(objects.index == range(len(objects))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
     if areas is not None:
         if not isinstance(areas, str):
@@ -373,10 +378,15 @@ def effective_mesh(objects, spatial_weights=None, areas=None, order=3):
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
         neighbours = spatial_weights.neighbors[index]
         if areas is not None:
-            total_area = sum(objects.iloc[neighbours][areas]) + row[areas]
+            area_list = objects.iloc[neighbours][areas].tolist()
+            area_list.append(row[areas])
         else:
-            total_area = sum(objects.iloc[neighbours].geometry.area) + row.geometry.area
-        results_list.append(total_area / (len(neighbours) + 1))
+            area_list = objects.iloc[neighbours].geometry.area.tolist()
+            area_list.append(row.geometry.area)
+        if mode:
+            from momepy import limit_range
+            area_list = limit_range(area_list, mode=mode)
+        results_list.append(sum(area_list) / len(area_list))
 
     series = pd.Series(results_list)
 
diff --git a/momepy/elements.py b/momepy/elements.py
index d61153e..6eda887 100644
--- a/momepy/elements.py
+++ b/momepy/elements.py
@@ -875,7 +875,7 @@ def get_node_id(objects, nodes, node_id, max_distance=200):
         GeoDataFrame containing objects
     nodes : GeoDataFrame
         GeoDataFrame containing street nodes with unique node ID.
-        If there is none, it could be generated by :py:func:`momepy.elements.unique_id`.
+        If there is none, it could be generated by :py:func:`momepy.unique_id`.
     node_id : str, list, np.array, pd.Series (default None)
         the name of the nodes dataframe column, np.array, or pd.Series with unique id
     max_distance : int (default=200)
diff --git a/momepy/utils.py b/momepy/utils.py
index 8afe89d..fa13351 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -6,6 +6,7 @@ import libpysal
 from shapely.geometry import Point
 import networkx as nx
 import pandas as pd
+import numpy as np
 
 
 def unique_id(objects):
@@ -38,7 +39,8 @@ def Queen_higher(k, geodataframe=None, weights=None):
     k : int
         order of contiguity
     geodataframe : GeoDataFrame
-        GeoDataFrame containing objects to analyse
+        GeoDataFrame containing objects to analyse. Index has to be consecutive range 0:x.
+        Otherwise, spatial weights will not match objects.
     weights : libpysal.weights
         libpysal.weights of order 1
 
@@ -60,6 +62,8 @@ def Queen_higher(k, geodataframe=None, weights=None):
     if weights is not None:
         first_order = weights
     elif geodataframe is not None:
+        if not all(geodataframe.index == range(len(geodataframe))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
         first_order = libpysal.weights.Queen.from_dataframe(geodataframe)
     else:
         raise Warning('GeoDataFrame of spatial weights must be given.')
@@ -173,3 +177,34 @@ def multi2single(gpdf):
 
     gpdf_singlepoly.reset_index(inplace=True, drop=True)
     return gpdf_singlepoly
+
+
+def limit_range(vals, mode):
+    """
+    Extract values within selected mode ('id', 'iq')
+
+    Parameters
+    ----------
+    vals : array
+
+    mode : str
+        'iq' for interquartile range, 'id' for interdecile range
+
+    Returns
+    -------
+    array
+        limited array
+    """
+    limited = []
+    if mode == 'iq':
+        mval = 25
+        xval = 75
+    elif mode == 'id':
+        mval = 10
+        xval = 90
+    lower = np.percentile(vals, mval)
+    higher = np.percentile(vals, xval)
+    for x in vals:
+        if x >= lower and x <= higher:
+            limited.append(x)
+    return limited


Test patch:
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 7c7ad2e..a0f03a4 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -101,36 +101,25 @@ class TestDimensions:
     def test_effective_mesh(self):
         self.df_tessellation['mesh'] = mm.effective_mesh(self.df_tessellation)
         spatial_weights = Queen_higher(k=3, geodataframe=self.df_tessellation)
+        self.df_tessellation['mesh_sw'] = mm.effective_mesh(self.df_tessellation, spatial_weights)
+        self.df_tessellation['area'] = area = self.df_tessellation.geometry.area
+        self.df_tessellation['mesh_ar'] = mm.effective_mesh(self.df_tessellation, spatial_weights, areas='area')
+        self.df_tessellation['mesh_array'] = mm.effective_mesh(self.df_tessellation, spatial_weights, areas=area)
+        self.df_tessellation['mesh_id'] = mm.effective_mesh(self.df_tessellation, spatial_weights, areas='area', mode='id')
+        self.df_tessellation['mesh_iq'] = mm.effective_mesh(self.df_tessellation, spatial_weights, areas='area', mode='iq')
         neighbours = spatial_weights.neighbors[38]
         total_area = sum(self.df_tessellation.iloc[neighbours].geometry.area) + self.df_tessellation.geometry.area[38]
         check = total_area / (len(neighbours) + 1)
         assert self.df_tessellation['mesh'][38] == check
-
-    def test_effective_mesh_sw(self):
-        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
-        self.df_tessellation['mesh'] = mm.effective_mesh(self.df_tessellation, sw)
-        neighbours = sw.neighbors[38]
-        total_area = sum(self.df_tessellation.iloc[neighbours].geometry.area) + self.df_tessellation.geometry.area[38]
-        check = total_area / (len(neighbours) + 1)
-        assert self.df_tessellation['mesh'][38] == check
-
-    def test_effective_mesh_area(self):
-        self.df_tessellation['area'] = self.df_tessellation.geometry.area
-        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
-        self.df_tessellation['mesh'] = mm.effective_mesh(self.df_tessellation, sw, areas='area')
-        neighbours = sw.neighbors[38]
-        total_area = sum(self.df_tessellation.iloc[neighbours].geometry.area) + self.df_tessellation.geometry.area[38]
-        check = total_area / (len(neighbours) + 1)
-        assert self.df_tessellation['mesh'][38] == check
-
-    def test_effective_mesh_array(self):
-        area = self.df_tessellation.geometry.area
-        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
-        self.df_tessellation['mesh'] = mm.effective_mesh(self.df_tessellation, sw, areas=area)
-        neighbours = sw.neighbors[38]
-        total_area = sum(self.df_tessellation.iloc[neighbours].geometry.area) + self.df_tessellation.geometry.area[38]
-        check = total_area / (len(neighbours) + 1)
-        assert self.df_tessellation['mesh'][38] == check
+        assert self.df_tessellation['mesh_sw'][38] == check
+        assert self.df_tessellation['mesh_ar'][38] == check
+        assert self.df_tessellation['mesh_array'][38] == check
+        assert self.df_tessellation['mesh_id'][38] == 2206.611646069303
+        assert self.df_tessellation['mesh_iq'][38] == 2118.609142733066
+        gdf = self.df_tessellation
+        gdf.index = gdf.index + 20
+        with pytest.raises(ValueError):
+            self.df_tessellation['meshE'] = mm.effective_mesh(gdf)
 
     def test_street_profile(self):
         widths, heights, profile = mm.street_profile(self.df_streets, self.df_buildings, heights='height')
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 6281559..6f993e9 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -33,6 +33,11 @@ class TestUtils:
         with pytest.raises(Warning):
             mm.Queen_higher(2, geodataframe=None, weights=None)
 
+        with pytest.raises(ValueError):
+            gdf = self.df_tessellation
+            gdf.index = gdf.index + 20
+            mm.Queen_higher(2, geodataframe=gdf, weights=None)
+
     def test_multi2single(self):
         polygon = Polygon([(0, 0), (1, 1), (1, 0)])
         polygon2 = Polygon([(2, 3), (1, 2), (2, 4)])
@@ -62,3 +67,7 @@ class TestUtils:
         assert W.n == 32
         nodes = mm.nx_to_gdf(nx, edges=False, spatial_weights=False)
         assert len(nodes) == 32
+
+    def test_limit_range(self):
+        assert mm.limit_range(range(10), mode='iq') == [3, 4, 5, 6]
+        assert mm.limit_range(range(10), mode='id') == [1, 2, 3, 4, 5, 6, 7, 8]