Repository: martinfleis/momepy (python)

Issue:
ENH: handle NaNs in limit_range
`limit_range` returns an empty array if there's NaN. It should be ignored (maybe add a keyword to control?).

Pull request:
ENH: handle NaNs in limit_range
Closes #227

Patch:
diff --git a/ci/envs/dev.yaml b/ci/envs/dev.yaml
index 369f086..ebb1120 100644
--- a/ci/envs/dev.yaml
+++ b/ci/envs/dev.yaml
@@ -10,7 +10,6 @@ dependencies:
   - pytest
   - pytest-cov
   - codecov
-  - mapclassify
   - pip
   - osmnx
   - black=19.10b0
@@ -18,6 +17,8 @@ dependencies:
   - inequality
   - pygeos
   - dask
+  # for mapclassify
+  - scikit-learn
   - pip:
       - git+https://github.com/geopandas/geopandas.git
       - git+https://github.com/pysal/libpysal.git
diff --git a/momepy/intensity.py b/momepy/intensity.py
index 3cf79c9..a5ee69f 100644
--- a/momepy/intensity.py
+++ b/momepy/intensity.py
@@ -278,7 +278,8 @@ class BlocksCount:
     """
     Calculates the weighted number of blocks
 
-    Number of blocks within neighbours defined in ``spatial_weights``.
+    Number of blocks within neighbours defined in ``spatial_weights`` divided by the area
+    you have covered by the neighbours.
 
     .. math::
 
diff --git a/momepy/preprocessing.py b/momepy/preprocessing.py
index 65b0d3b..1380fa7 100644
--- a/momepy/preprocessing.py
+++ b/momepy/preprocessing.py
@@ -192,7 +192,7 @@ def remove_false_nodes(gdf):
     Returns
     -------
     gdf : GeoDataFrame, GeoSeries
-    
+
     See also
     --------
     momepy.extend_lines
@@ -343,7 +343,7 @@ class CheckTessellationInput:
     have to be fixed prior Tessellation. Features which will split will cause issues
     only sometimes, so
     should be checked and fixed if necessary. Features which will collapse could
-    be ignored, but they will have to excluded from next steps of 
+    be ignored, but they will have to excluded from next steps of
     tessellation-based analysis.
 
     Parameters
@@ -570,7 +570,7 @@ def close_gaps(gdf, tolerance):
     Returns
     -------
     GeoSeries
-    
+
     See also
     --------
     momepy.extend_lines
@@ -609,7 +609,7 @@ def close_gaps(gdf, tolerance):
 
 
 def extend_lines(gdf, tolerance, target=None, barrier=None, extension=0):
-    """ Extends lines from gdf to istelf or target within a set tolerance
+    """Extends lines from gdf to istelf or target within a set tolerance
 
     Extends unjoined ends of LineString segments to join with other segments or
     target. If ``target`` is passed, extend lines to target. Otherwise extend
@@ -619,7 +619,7 @@ def extend_lines(gdf, tolerance, target=None, barrier=None, extension=0):
     with ``barrier``. If they intersect, extended line is not returned. This
     can be useful if you don't want to extend street network segments through
     buildings.
-    
+
     Parameters
     ----------
     gdf : GeoDataFrame
@@ -628,19 +628,19 @@ def extend_lines(gdf, tolerance, target=None, barrier=None, extension=0):
         tolerance in snapping (by how much could be each segment
         extended).
     target : GeoDataFrame, GeoSeries
-        target geometry to which ``gdf`` gets extended. Has to be 
+        target geometry to which ``gdf`` gets extended. Has to be
         (Multi)LineString geometry.
     barrier : GeoDataFrame, GeoSeries
         extended line is not used if it intersects barrier
     extension : float
         by how much to extend line beyond the snapped geometry. Useful
         when creating enclosures to avoid floating point imprecision.
-    
+
     Returns
     -------
     GeoDataFrame
         GeoDataFrame of with extended geometry
-    
+
     See also
     --------
     momepy.close_gaps
diff --git a/momepy/shape.py b/momepy/shape.py
index 43876a1..2e1c9fd 100644
--- a/momepy/shape.py
+++ b/momepy/shape.py
@@ -414,7 +414,7 @@ class CircularCompactness:
     gdf : GeoDataFrame
         GeoDataFrame containing objects
     areas : str, list, np.array, pd.Series (default None)
-        the name of the dataframe column, ``np.array``, or ``pd.Series`` where is 
+        the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
         stored area value. If set to ``None``, function will calculate areas
         during the process without saving them separately.
 
@@ -727,7 +727,7 @@ class Corners:
     """
     Calculates number of corners of each object in given GeoDataFrame.
 
-    Uses only external shape (``shapely.geometry.exterior``), courtyards are not 
+    Uses only external shape (``shapely.geometry.exterior``), courtyards are not
     included.
 
     .. math::
@@ -781,33 +781,65 @@ class Corners:
 
         # fill new column with the value of area, iterating over rows one by one
         for geom in tqdm(gdf.geometry, total=gdf.shape[0], disable=not verbose):
-            corners = 0  # define empty variables
-            points = list(geom.exterior.coords)  # get points of a shape
-            stop = len(points) - 1  # define where to stop
-            for i in np.arange(
-                len(points)
-            ):  # for every point, calculate angle and add 1 if True angle
-                if i == 0:
-                    continue
-                elif i == stop:
-                    a = np.asarray(points[i - 1])
-                    b = np.asarray(points[i])
-                    c = np.asarray(points[1])
-
-                    if _true_angle(a, b, c) is True:
-                        corners = corners + 1
-                    else:
+            if geom.type == "Polygon":
+                corners = 0  # define empty variables
+                points = list(geom.exterior.coords)  # get points of a shape
+                stop = len(points) - 1  # define where to stop
+                for i in np.arange(
+                    len(points)
+                ):  # for every point, calculate angle and add 1 if True angle
+                    if i == 0:
                         continue
+                    elif i == stop:
+                        a = np.asarray(points[i - 1])
+                        b = np.asarray(points[i])
+                        c = np.asarray(points[1])
 
-                else:
-                    a = np.asarray(points[i - 1])
-                    b = np.asarray(points[i])
-                    c = np.asarray(points[i + 1])
+                        if _true_angle(a, b, c) is True:
+                            corners = corners + 1
+                        else:
+                            continue
 
-                    if _true_angle(a, b, c) is True:
-                        corners = corners + 1
                     else:
-                        continue
+                        a = np.asarray(points[i - 1])
+                        b = np.asarray(points[i])
+                        c = np.asarray(points[i + 1])
+
+                        if _true_angle(a, b, c) is True:
+                            corners = corners + 1
+                        else:
+                            continue
+            elif geom.type == "MultiPolygon":
+                corners = 0  # define empty variables
+                for g in geom.geoms:
+                    points = list(g.exterior.coords)  # get points of a shape
+                    stop = len(points) - 1  # define where to stop
+                    for i in np.arange(
+                        len(points)
+                    ):  # for every point, calculate angle and add 1 if True angle
+                        if i == 0:
+                            continue
+                        elif i == stop:
+                            a = np.asarray(points[i - 1])
+                            b = np.asarray(points[i])
+                            c = np.asarray(points[1])
+
+                            if _true_angle(a, b, c) is True:
+                                corners = corners + 1
+                            else:
+                                continue
+
+                        else:
+                            a = np.asarray(points[i - 1])
+                            b = np.asarray(points[i])
+                            c = np.asarray(points[i + 1])
+
+                            if _true_angle(a, b, c) is True:
+                                corners = corners + 1
+                            else:
+                                continue
+            else:
+                corners = np.nan
 
             results_list.append(corners)
 
@@ -827,6 +859,8 @@ class Squareness:
 
     Adapted from :cite:`dibble2017`.
 
+    Returns ``np.nan`` for MultiPolygons.
+
     Parameters
     ----------
     gdf : GeoDataFrame
@@ -1037,6 +1071,8 @@ class CentroidCorners:
 
     Adapted from :cite:`schirmer2015` and :cite:`cimburova2017`.
 
+    Returns ``np.nan`` for MultiPolygons.
+
     Parameters
     ----------
     gdf : GeoDataFrame
@@ -1087,52 +1123,57 @@ class CentroidCorners:
 
         # iterating over rows one by one
         for geom in tqdm(gdf.geometry, total=gdf.shape[0], disable=not verbose):
-            distances = []  # set empty list of distances
-            centroid = geom.centroid  # define centroid
-            points = list(geom.exterior.coords)  # get points of a shape
-            stop = len(points) - 1  # define where to stop
-            for i in np.arange(
-                len(points)
-            ):  # for every point, calculate angle and add 1 if True angle
-                if i == 0:
-                    continue
-                elif i == stop:
-                    a = np.asarray(points[i - 1])
-                    b = np.asarray(points[i])
-                    c = np.asarray(points[1])
-                    p = Point(points[i])
-
-                    if true_angle(a, b, c) is True:
-                        distance = centroid.distance(
-                            p
-                        )  # calculate distance point - centroid
-                        distances.append(distance)  # add distance to the list
-                    else:
+            if geom.type == "Polygon":
+                distances = []  # set empty list of distances
+                centroid = geom.centroid  # define centroid
+                points = list(geom.exterior.coords)  # get points of a shape
+                stop = len(points) - 1  # define where to stop
+                for i in np.arange(
+                    len(points)
+                ):  # for every point, calculate angle and add 1 if True angle
+                    if i == 0:
                         continue
+                    elif i == stop:
+                        a = np.asarray(points[i - 1])
+                        b = np.asarray(points[i])
+                        c = np.asarray(points[1])
+                        p = Point(points[i])
+
+                        if true_angle(a, b, c) is True:
+                            distance = centroid.distance(
+                                p
+                            )  # calculate distance point - centroid
+                            distances.append(distance)  # add distance to the list
+                        else:
+                            continue
 
-                else:
-                    a = np.asarray(points[i - 1])
-                    b = np.asarray(points[i])
-                    c = np.asarray(points[i + 1])
-                    p = Point(points[i])
-
-                    if true_angle(a, b, c) is True:
-                        distance = centroid.distance(p)
-                        distances.append(distance)
                     else:
-                        continue
-            if not distances:  # circular buildings
-                if geom.has_z:
-                    coords = [
-                        (coo[0], coo[1]) for coo in geom.convex_hull.exterior.coords
-                    ]
+                        a = np.asarray(points[i - 1])
+                        b = np.asarray(points[i])
+                        c = np.asarray(points[i + 1])
+                        p = Point(points[i])
+
+                        if true_angle(a, b, c) is True:
+                            distance = centroid.distance(p)
+                            distances.append(distance)
+                        else:
+                            continue
+                if not distances:  # circular buildings
+                    if geom.has_z:
+                        coords = [
+                            (coo[0], coo[1]) for coo in geom.convex_hull.exterior.coords
+                        ]
+                    else:
+                        coords = geom.convex_hull.exterior.coords
+                    results_list.append(_circle_radius(coords))
+                    results_list_sd.append(0)
                 else:
-                    coords = geom.convex_hull.exterior.coords
-                results_list.append(_circle_radius(coords))
-                results_list_sd.append(0)
+                    results_list.append(np.mean(distances))  # calculate mean
+                    results_list_sd.append(np.std(distances))  # calculate st.dev
             else:
-                results_list.append(np.mean(distances))  # calculate mean
-                results_list_sd.append(np.std(distances))  # calculate st.dev
+                results_list.append(np.nan)
+                results_list_sd.append(np.nan)
+
         self.mean = pd.Series(results_list, index=gdf.index)
         self.std = pd.Series(results_list_sd, index=gdf.index)
 
diff --git a/momepy/utils.py b/momepy/utils.py
index d9b2a90..c4a4528 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -287,12 +287,16 @@ def limit_range(vals, rng):
     array
         limited array
     """
+    vals = np.asarray(vals)
     if len(vals) > 2:
         rng = sorted(rng)
-        lower = np.percentile(vals, rng[0], interpolation="nearest")
-        higher = np.percentile(vals, rng[1], interpolation="nearest")
-        limited = [x for x in vals if x >= lower and x <= higher]
-        return np.array(limited)
+        if np.isnan(vals).any():
+            lower = np.nanpercentile(vals, rng[0], interpolation="nearest")
+            higher = np.nanpercentile(vals, rng[1], interpolation="nearest")
+        else:
+            lower = np.percentile(vals, rng[0], interpolation="nearest")
+            higher = np.percentile(vals, rng[1], interpolation="nearest")
+        return vals[(lower <= vals) & (vals <= higher)]
     return vals
 
 


Test patch:
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 94f0427..84c9fee 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -95,3 +95,6 @@ class TestUtils:
         assert list(mm.limit_range(range(10), rng=(25, 75))) == [2, 3, 4, 5, 6, 7]
         assert list(mm.limit_range(range(10), rng=(10, 90))) == [1, 2, 3, 4, 5, 6, 7, 8]
         assert list(mm.limit_range([0, 1], rng=(25, 75))) == [0, 1]
+        assert list(
+            mm.limit_range(np.array([0, 1, 2, 3, 4, np.nan]), rng=(25, 75))
+        ) == [1, 2, 3]