Repository: martinfleis/momepy
Pull request #186 (2020-07-17 20:13:47)

REF: StreetProfile performance refactor
Refactor of `StreetProfile` to get better performance. There's still space for more work.

Benchmark:
```
       before           after         ratio
     [05ca8623]       [fece9f90]
     <master>         <profile_perf~3>
-      3.17±0.01s       1.18±0.05s     0.37  bench_dimension.TimeDimension.time_StreetProfile
```

I've also fixed a bug I did not know about which sometimes picked incorrect height.


Patch:
diff --git a/ci/travis/37-pysal.yaml b/ci/travis/37-pysal.yaml
index f45e62c..723160e 100644
--- a/ci/travis/37-pysal.yaml
+++ b/ci/travis/37-pysal.yaml
@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python=3.7
-  - geopandas
+  - geopandas=0.7
   - libpysal
   - networkx
   - tqdm
diff --git a/momepy/__init__.py b/momepy/__init__.py
index 3eb2d08..b4856b5 100644
--- a/momepy/__init__.py
+++ b/momepy/__init__.py
@@ -1,13 +1,14 @@
+import momepy.datasets
+
 from .dimension import *
-from .intensity import *
-from .utils import *
-from .elements import *
-from .diversity import *
 from .distribution import *
+from .diversity import *
+from .elements import *
 from .graph import *
+from .intensity import *
 from .shape import *
+from .utils import *
 from .weights import *
-import momepy.datasets
 
 __author__ = "Martin Fleischmann"
 __author_email__ = "martin@martinfleischmann.net"
diff --git a/momepy/datasets/__init__.py b/momepy/datasets/__init__.py
index eb45339..4e7a5b9 100644
--- a/momepy/datasets/__init__.py
+++ b/momepy/datasets/__init__.py
@@ -1,6 +1,5 @@
 import os
 
-
 __all__ = ["available", "get_path"]
 
 _module_path = os.path.dirname(__file__)
diff --git a/momepy/dimension.py b/momepy/dimension.py
index 44b5541..8aaf0f0 100644
--- a/momepy/dimension.py
+++ b/momepy/dimension.py
@@ -5,7 +5,9 @@
 # definitions of dimension characters
 
 import math
+from distutils.version import LooseVersion
 
+import geopandas as gpd
 import numpy as np
 import pandas as pd
 import scipy as sp
@@ -14,6 +16,8 @@ from tqdm import tqdm
 
 from .shape import _make_circle
 
+GPD_08 = str(gpd.__version__) >= LooseVersion("0.8.0")
+
 __all__ = [
     "Area",
     "Perimeter",
@@ -538,25 +542,25 @@ class StreetProfile:
         heights_deviations_list = []
         openness_list = []
 
-        for idx, row in tqdm(left.iterrows(), total=left.shape[0]):
+        for shapely_line in tqdm(left.geometry, total=left.shape[0]):
             # list to hold all the point coords
             list_points = []
             # set the current distance to place the point
             current_dist = distance
-            # make shapely MultiLineString object
-            shapely_line = row.geometry
             # get the total length of the line
             line_length = shapely_line.length
             # append the starting coordinate to the list
-            list_points.append(Point(list(shapely_line.coords)[0]))
+            list_points.append(list(shapely_line.coords)[0])
             # https://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/
             # while the current cumulative distance is less than the total length of the line
             while current_dist < line_length:
                 # use interpolate and increase the current distance
-                list_points.append(shapely_line.interpolate(current_dist))
+                list_points.append(
+                    list(shapely_line.interpolate(current_dist).coords)[0]
+                )
                 current_dist += distance
             # append end coordinate to the list
-            list_points.append(Point(list(shapely_line.coords)[-1]))
+            list_points.append(list(shapely_line.coords)[-1])
 
             ticks = []
             for num, pt in enumerate(list_points, 1):
@@ -566,8 +570,8 @@ class StreetProfile:
                     line_end_1 = self._getPoint1(pt, angle, tick_length / 2)
                     angle = self._getAngle(line_end_1, pt)
                     line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
-                    tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
-                    tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
+                    tick1 = LineString([line_end_1, pt])
+                    tick2 = LineString([line_end_2, pt])
                     ticks.append([tick1, tick2])
 
                 # everything in between
@@ -578,18 +582,8 @@ class StreetProfile:
                     )
                     angle = self._getAngle(line_end_1, list_points[num])
                     line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
-                    tick1 = LineString(
-                        [
-                            (line_end_1.x, line_end_1.y),
-                            (list_points[num].x, list_points[num].y),
-                        ]
-                    )
-                    tick2 = LineString(
-                        [
-                            (line_end_2.x, line_end_2.y),
-                            (list_points[num].x, list_points[num].y),
-                        ]
-                    )
+                    tick1 = LineString([line_end_1, list_points[num]])
+                    tick2 = LineString([line_end_2, list_points[num]])
                     ticks.append([tick1, tick2])
 
                 # end chainage
@@ -598,61 +592,35 @@ class StreetProfile:
                     line_end_1 = self._getPoint1(pt, angle, tick_length / 2)
                     angle = self._getAngle(line_end_1, pt)
                     line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
-                    tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
-                    tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
+                    tick1 = LineString([line_end_1, pt])
+                    tick2 = LineString([line_end_2, pt])
                     ticks.append([tick1, tick2])
             # widths = []
             m_heights = []
             lefts = []
             rights = []
             for duo in ticks:
-
                 for ix, tick in enumerate(duo):
-                    possible_intersections_index = list(
-                        sindex.intersection(tick.bounds)
-                    )
-                    possible_intersections = right.iloc[possible_intersections_index]
-                    real_intersections = possible_intersections.intersects(tick)
-                    get_height = right.loc[list(real_intersections.index)]
-                    possible_int = get_height.exterior.intersection(tick)
-
-                    if not possible_int.is_empty.all():
-                        true_int = []
-                        for one in list(possible_int.index):
-                            if possible_int[one].type == "Point":
-                                true_int.append(possible_int[one])
-                            elif possible_int[one].type == "MultiPoint":
-                                for p in possible_int[one]:
-                                    true_int.append(p)
-
-                        if len(true_int) > 1:
-                            distances = []
-                            ix = 0
-                            for p in true_int:
-                                dist = p.distance(Point(tick.coords[-1]))
-                                distances.append(dist)
-                                ix = ix + 1
-                            minimal = min(distances)
-                            if ix == 0:
-                                lefts.append(minimal)
-                            else:
-                                rights.append(minimal)
+                    if GPD_08:
+                        int_blg = right.iloc[sindex.query(tick, predicate="intersects")]
+                    else:
+                        possible_intersections_index = list(
+                            sindex.intersection(tick.bounds)
+                        )
+                        possible_intersections = right.iloc[
+                            possible_intersections_index
+                        ]
+                        real_intersections = possible_intersections.intersects(tick)
+                        int_blg = possible_intersections[real_intersections]
+                    if not int_blg.empty:
+                        true_int = int_blg.intersection(tick)
+                        dist = true_int.distance(Point(tick.coords[-1]))
+                        if ix == 0:
+                            lefts.append(dist.min())
                         else:
-                            if ix == 0:
-                                lefts.append(
-                                    true_int[0].distance(Point(tick.coords[-1]))
-                                )
-                            else:
-                                rights.append(
-                                    true_int[0].distance(Point(tick.coords[-1]))
-                                )
+                            rights.append(dist.min())
                         if heights is not None:
-                            indices = {}
-                            for idx, row in get_height.iterrows():
-                                dist = row.geometry.distance(Point(tick.coords[-1]))
-                                indices[idx] = dist
-                            minim = min(indices, key=indices.get)
-                            m_heights.append(right.loc[minim][heights])
+                            m_heights.append(int_blg.loc[dist.idxmin()][heights])
 
             openness = (len(lefts) + len(rights)) / len(ticks * 2)
             openness_list.append(1 - openness)
@@ -690,25 +658,34 @@ class StreetProfile:
     # https://glenbambrick.com/tag/perpendicular/
     # angle between two points
     def _getAngle(self, pt1, pt2):
-        x_diff = pt2.x - pt1.x
-        y_diff = pt2.y - pt1.y
+        """
+        pt1, pt2 : tuple
+        """
+        x_diff = pt2[0] - pt1[0]
+        y_diff = pt2[1] - pt1[1]
         return math.degrees(math.atan2(y_diff, x_diff))
 
     # start and end points of chainage tick
     # get the first end point of a tick
     def _getPoint1(self, pt, bearing, dist):
+        """
+        pt : tuple
+        """
         angle = bearing + 90
         bearing = math.radians(angle)
-        x = pt.x + dist * math.cos(bearing)
-        y = pt.y + dist * math.sin(bearing)
-        return Point(x, y)
+        x = pt[0] + dist * math.cos(bearing)
+        y = pt[1] + dist * math.sin(bearing)
+        return (x, y)
 
     # get the second end point of a tick
     def _getPoint2(self, pt, bearing, dist):
+        """
+        pt : tuple
+        """
         bearing = math.radians(bearing)
-        x = pt.x + dist * math.cos(bearing)
-        y = pt.y + dist * math.sin(bearing)
-        return Point(x, y)
+        x = pt[0] + dist * math.cos(bearing)
+        y = pt[1] + dist * math.sin(bearing)
+        return (x, y)
 
 
 class WeightedCharacter:
diff --git a/momepy/distribution.py b/momepy/distribution.py
index 71b6b0d..0891d0a 100644
--- a/momepy/distribution.py
+++ b/momepy/distribution.py
@@ -13,7 +13,6 @@ from tqdm import tqdm  # progress bar
 
 from .utils import _azimuth
 
-
 __all__ = [
     "Orientation",
     "SharedWallsRatio",
@@ -570,9 +569,6 @@ class MeanInterbuildingDistance:
     --------
     >>> buildings_df['mean_interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings_df, sw, 'uID').series
     Generating weights matrix (Queen) of 3 topological steps...
-    Generating adjacency matrix based on weights matrix...
-    Computing interbuilding distances...
-    100%|██████████| 746/746 [00:03<00:00, 200.14it/s]
     Computing mean interbuilding distances...
     100%|██████████| 144/144 [00:00<00:00, 317.42it/s]
     >>> buildings_df['mean_interbuilding_distance'][0]
@@ -604,26 +600,13 @@ class MeanInterbuildingDistance:
         # define empty list for results
         results_list = []
 
-        print("Generating adjacency matrix based on weights matrix...")
         # define adjacency list from lipysal
         adj_list = spatial_weights.to_adjlist()
-        adj_list["distance"] = -1
-
-        print("Computing interbuilding distances...")
-        # measure each interbuilding distance of neighbours and save them to adjacency list
-        for row in tqdm(adj_list.itertuples(), total=adj_list.shape[0]):
-            inverted = adj_list[(adj_list.focal == row.neighbor)][
-                (adj_list.neighbor == row.focal)
-            ].iloc[0]["distance"]
-            if inverted == -1:
-                building_object = data.loc[row.focal]
-
-                building_neighbour = data.loc[row.neighbor]
-                adj_list.loc[row.Index, "distance"] = building_neighbour.distance(
-                    building_object
-                )
-            else:
-                adj_list.at[row.Index, "distance"] = inverted
+        adj_list["distance"] = (
+            data.loc[adj_list.focal]
+            .reset_index()
+            .distance(data.loc[adj_list.neighbor].reset_index())
+        )
 
         print("Computing mean interbuilding distances...")
         # iterate over objects to get the final values
diff --git a/momepy/elements.py b/momepy/elements.py
index 6f525ea..50d3fea 100644
--- a/momepy/elements.py
+++ b/momepy/elements.py
@@ -6,16 +6,15 @@
 import operator
 
 import geopandas as gpd
+import libpysal
 import numpy as np
 import pandas as pd
 import shapely
-import libpysal
 from scipy.spatial import Voronoi
 from shapely.geometry import MultiPolygon, Point, Polygon
 from shapely.wkt import loads
 from tqdm import tqdm
 
-
 __all__ = ["buffered_limit", "Tessellation", "Blocks", "get_network_id", "get_node_id"]
 
 
diff --git a/momepy/weights.py b/momepy/weights.py
index 920eeb4..b7274bd 100644
--- a/momepy/weights.py
+++ b/momepy/weights.py
@@ -4,7 +4,6 @@
 import libpysal
 import numpy as np
 
-
 __all__ = ["DistanceBand", "sw_high"]
 
 


Test patch:
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 0dc708b..7d22af8 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -205,10 +205,10 @@ class TestDimensions:
         )
         assert results2.w[0] == 67.563796073073
         assert results2.wd[0] == 8.791875291865827
-        assert results2.h[0] == 23.756643356643362
-        assert results2.p[0] == 0.3516179483306353
+        assert results2.h[0] == 23.74545454545455
+        assert results2.p[0] == 0.3514523446813568
         assert results2.o[0] == 0.5535714285714286
-        assert results2.hd[0] == 5.526848034418866
+        assert results2.hd[0] == 5.5188633970914065
 
         results3 = mm.StreetProfile(self.df_streets, self.df_buildings)
         assert results3.w[0] == 46.7290758769204