Repository: martinfleis/momepy
Pull request #37 (2019-07-29 11:17:15)

BUG: use explode() in the preprocess to avoid multigeometry
In case of MultiPolygon geometry, `preprocess` may fail. This ensures that there is only Polygon geometry before running actual cleaning process.


Patch:
diff --git a/momepy/elements.py b/momepy/elements.py
index 5a7478f..c5aa41d 100644
--- a/momepy/elements.py
+++ b/momepy/elements.py
@@ -19,82 +19,63 @@ import operator
 from libpysal.weights import Queen
 
 
-def tessellation(gdf, unique_id='uID', cut_buffer=50, queen_corners=False, minimum=2):
+def buffered_limit(gdf, buffer=100):
     """
-    Generate morphological tessellation around given buildings.
+    Define limit for tessellation as a buffer around buildings.
 
     Parameters
     ----------
     gdf : GeoDataFrame
         GeoDataFrame containing building footprints
-    unique_id : str
-        name of the column with unique id. If there is none, it could be generated by unique_id().
-    cut_buffer : float
+    buffer : float
         buffer around buildings limiting the extend of tessellation
 
     Returns
     -------
-    GeoDataFrame
-        GeoDataFrame of morphological tessellation with the unique id based on original buildings.
+    MultiPolygon
+        MultiPolygon or Polygon defining the study area
 
-    Notes
-    -------
-    Fix saw-like geometry.
     """
+    study_area = gdf.copy()
+    study_area['geometry'] = study_area.buffer(100)
+    study_area['diss'] = 1
+    built_up_df = study_area.dissolve(by='diss')
+    built_up = built_up_df.geometry[1]
+    return built_up
 
-    tqdm.pandas()
-    objects = gdf.copy()
-    # move dataframe close to 0 to eliminate imprecision of Qhull
-
-    def _get_centre(gdf):
-        bounds = gdf['geometry'].bounds
-        centre_x = (bounds['maxx'].max() + bounds['minx'].min()) / 2
-        centre_y = (bounds['maxy'].max() + bounds['miny'].min()) / 2
-        return centre_x, centre_y
-
-    centre = _get_centre(gdf)
-    objects['geometry'] = objects['geometry'].translate(xoff=-centre[0], yoff=-centre[1])
-
-    # buffer geometry to resolve shared walls
-    print('Bufferring geometry...')
-    objects['geometry'] = objects.geometry.apply(lambda g: g.buffer(-0.5, cap_style=2, join_style=2))
-
-    # simplify geometry before Voronoi
-    print('Simplifying geometry...')
-    objects['geometry'] = objects.simplify(0.25, preserve_topology=True)
-    obj_simple = objects.copy()
-
-    print('Preparing buffer zone for edge resolving (buffering)...')
-    obj_simple['geometry'] = obj_simple.buffer(cut_buffer)
-
-    print('Preparing buffer zone for edge resolving (dissolving)...')
-    obj_simple['diss'] = 0
-    built_up_df = obj_simple.dissolve(by='diss')
-    built_up = built_up_df.geometry[0]
-
-    # resolve multipart polygons, singlepart are needed for densification
-    print('Converting multipart geometry to singlepart...')
-    objects = objects.explode()
-    objects.reset_index(inplace=True, drop=True)
 
-    # densify geometry before Voronoi tesselation
-    def _densify(geom):
-        poly = geom
-        wkt = geom.wkt  # shapely Polygon to wkt
-        geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
-        geom.Segmentize(2)  # densify geometry by 2 metres
-        geom.CloseRings()  # fix for GDAL 2.4.1 bug
-        wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
-        try:
-            new = loads(wkt2)  # wkt to shapely Polygon
-            return new
-        except:
-            return poly
+def _get_centre(gdf):
+    """
+    Returns centre coords of gdf.
+    """
+    bounds = gdf['geometry'].bounds
+    centre_x = (bounds['maxx'].max() + bounds['minx'].min()) / 2
+    centre_y = (bounds['maxy'].max() + bounds['miny'].min()) / 2
+    return centre_x, centre_y
 
-    print('Densifying geometry...')
-    objects['geometry'] = objects['geometry'].progress_map(_densify)
 
-    print('Generating input point array...')
+# densify geometry before Voronoi tesselation
+def _densify(geom, segment):
+    """
+    Returns densified geoemtry with segments no longer than `segment`.
+    """
+    poly = geom
+    wkt = geom.wkt  # shapely Polygon to wkt
+    geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
+    geom.Segmentize(segment)  # densify geometry by 2 metres
+    geom.CloseRings()  # fix for GDAL 2.4.1 bug
+    wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
+    try:
+        new = loads(wkt2)  # wkt to shapely Polygon
+        return new
+    except:
+        return poly
+
+
+def _point_array(objects, unique_id):
+    """
+    Returns lists of points and ids based on geometry and unique_id.
+    """
     points = []
     ids = []
     for idx, row in tqdm(objects.iterrows(), total=objects.shape[0]):
@@ -115,21 +96,13 @@ def tessellation(gdf, unique_id='uID', cut_buffer=50, queen_corners=False, minim
                     ids.append(row[unique_id])
             else:
                 raise Exception('Boundary type is {}'.format(poly_ext.type))
+    return points, ids
 
-    # add convex hull buffered by cut distance to eliminate infinity issues
-    print('Generating convex hull...')
-    hull = built_up.convex_hull.buffer(cut_buffer)
-    hull_array = np.array(hull.boundary.coords).tolist()
-    for i in range(len(hull_array)):
-        points.append(hull_array[i])
-        ids.append(-1)
-
-    voronoi_points = np.array(points)  # array of points as an input of Voronoi
 
-    print('Generating Voronoi diagram...')
-    voronoi_diagram = Voronoi(voronoi_points)
-
-    print('Generating GeoDataFrame...')
+def _regions(voronoi_diagram, unique_id, ids, crs):
+    """
+    Generate GeoDataFrame of Voronoi regions from scipy.spatial.Voronoi.
+    """
     # generate DataFrame of results
     regions = pd.DataFrame()
     regions[unique_id] = ids  # add unique id
@@ -155,24 +128,29 @@ def tessellation(gdf, unique_id='uID', cut_buffer=50, queen_corners=False, minim
     regions_gdf = gpd.GeoDataFrame(regions.dropna(), geometry='geometry')
     regions_gdf = regions_gdf.loc[regions_gdf['geometry'].length < 1000000]  # delete errors
     regions_gdf = regions_gdf.loc[regions_gdf[unique_id] != -1]  # delete hull-based cells
-    regions_gdf.crs = gdf.crs
+    regions_gdf.crs = crs
+    return regions_gdf
 
-    print('Dissolving Voronoi polygons...')
-    morphological_tessellation = regions_gdf[[unique_id, 'geometry']].dissolve(by=unique_id, as_index=False)
 
+def _cut(tessellation, limit, unique_id):
+    """
+    Cut tessellation by the limit (Multi)Polygon.
+
+    ADD: add option to delete everything outside of limit. Now it keeps it.
+    """
     # cut infinity of voronoi by set buffer (thanks for script to Geoff Boeing)
     print('Preparing buffer zone for edge resolving (quadrat cut)...')
-    geometry = built_up.boundary
+    geometry = limit.boundary
     geometry_cut = ox.quadrat_cut_geometry(geometry, quadrat_width=100)
 
     print('Building R-tree...')
-    sindex = morphological_tessellation.sindex
+    sindex = tessellation.sindex
     # find the points that intersect with each subpolygon and add them to points_within_geometry
     to_cut = pd.DataFrame()
     for poly in geometry_cut:
         # find approximate matches with r-tree, then precise matches from those approximate ones
         possible_matches_index = list(sindex.intersection(poly.bounds))
-        possible_matches = morphological_tessellation.iloc[possible_matches_index]
+        possible_matches = tessellation.iloc[possible_matches_index]
         precise_matches = possible_matches[possible_matches.intersects(poly)]
         to_cut = to_cut.append(precise_matches)
 
@@ -181,27 +159,33 @@ def tessellation(gdf, unique_id='uID', cut_buffer=50, queen_corners=False, minim
     subselection = list(to_cut.index)
 
     print('Cutting...')
-    for idx, row in tqdm(morphological_tessellation.loc[subselection].iterrows(), total=morphological_tessellation.loc[subselection].shape[0]):
-        intersection = row.geometry.intersection(built_up)
+    for idx, row in tqdm(tessellation.loc[subselection].iterrows(), total=tessellation.loc[subselection].shape[0]):
+        intersection = row.geometry.intersection(limit)
         if intersection.type == 'MultiPolygon':
             areas = {}
             for p in range(len(intersection)):
                 area = intersection[p].area
                 areas[p] = area
             maximal = max(areas.items(), key=operator.itemgetter(1))[0]
-            morphological_tessellation.loc[idx, 'geometry'] = intersection[maximal]
+            tessellation.loc[idx, 'geometry'] = intersection[maximal]
         elif intersection.type == 'GeometryCollection':
             for geom in list(intersection.geoms):
                 if geom.type != 'Polygon':
                     pass
                 else:
-                    morphological_tessellation.loc[idx, 'geometry'] = geom
+                    tessellation.loc[idx, 'geometry'] = geom
         else:
-            morphological_tessellation.loc[idx, 'geometry'] = intersection
+            tessellation.loc[idx, 'geometry'] = intersection
+    return tessellation, sindex
+
 
+def _check_result(tesselation, orig_gdf, unique_id):
+    """
+    Check whether result of tessellation matches buildings and contains only Polygons.
+    """
     # check against input layer
-    ids_original = list(gdf[unique_id])
-    ids_generated = list(morphological_tessellation[unique_id])
+    ids_original = list(orig_gdf[unique_id])
+    ids_generated = list(tesselation[unique_id])
     if len(ids_original) != len(ids_generated):
         import warnings
         diff = set(ids_original).difference(ids_generated)
@@ -209,112 +193,165 @@ def tessellation(gdf, unique_id='uID', cut_buffer=50, queen_corners=False, minim
                       "during generation - unique_id: {i}".format(len=len(diff), i=diff))
 
     # check MultiPolygons - usually caused by error in input geometry
-    uids = morphological_tessellation[morphological_tessellation.geometry.type == 'MultiPolygon'][unique_id]
+    uids = tesselation[tesselation.geometry.type == 'MultiPolygon'][unique_id]
     if len(uids) > 0:
         import warnings
         warnings.warn('Tessellation contains MultiPolygon elements. Initial objects should be edited. '
                       'unique_id of affected elements: {}'.format(list(uids)))
 
-    # resolve unprecise corners
-    if queen_corners is True:
-        print('Generating queen corners...')
-        print(' Generating spatial index...')
-        changes = {}
-        qid = 0
-        # detect points which should be changed and calculate new coordinates
-        print(' Detecting points of change...')
-        for ix, row in tqdm(morphological_tessellation.iterrows(), total=morphological_tessellation.shape[0]):
-            corners = []
-            change = []
-
-            cell = row.geometry
-            coords = cell.exterior.coords
-            for i in coords:
-                point = Point(i)
-                possible_matches_index = list(sindex.intersection(point.bounds))
-                possible_matches = morphological_tessellation.iloc[possible_matches_index]
-                precise_matches = sum(possible_matches.intersects(point))
-                if precise_matches > 2:
-                    corners.append(point)
-
-            if len(corners) > 2:
-                for c in range(len(corners)):
-                    next_c = c + 1
-                    if c == (len(corners) - 1):
-                        next_c = 0
-                    if corners[c].distance(corners[next_c]) < minimum:
-                        change.append([corners[c], corners[next_c]])
-            elif len(corners) == 2:
-                if corners[0].distance(corners[1]) > 0:
-                    if corners[0].distance(corners[1]) < minimum:
-                        change.append([corners[0], corners[1]])
-
-            if change:
-                for points in change:
-                    x_new = np.mean([points[0].x, points[1].x])
-                    y_new = np.mean([points[0].y, points[1].y])
-                    new = [(x_new, y_new), id]
-                    changes[(points[0].x, points[0].y)] = new
-                    changes[(points[1].x, points[1].y)] = new
-                    qid = qid + 1
-
-        print(' Generating new geometry...')
-        for ix, row in tqdm(morphological_tessellation.iterrows(), total=morphological_tessellation.shape[0]):
-            cell = row.geometry
-            coords = list(cell.exterior.coords)
-
-            moves = {}
-            for x in coords:
-                if x in changes.keys():
-                    moves[coords.index(x)] = changes[x]
-            keys = list(moves.keys())
-            delete_points = []
-            for move in range(len(keys)):
-                if move < len(keys) - 1:
-                    if moves[keys[move]][1] == moves[keys[move + 1]][1] and keys[move + 1] - keys[move] < 5:
-                        delete_points = delete_points + (coords[keys[move]:keys[move + 1]])
-                        # change the code above to have if based on distance not number
-
-            newcoords = [changes[x][0] if x in changes.keys() else x for x in coords]
-            for coord in newcoords:
-                if coord in delete_points:
-                    newcoords.remove(coord)
-            if coords != newcoords:
-                if not cell.interiors:
-                    # newgeom = Polygon(newcoords).buffer(0)
-                    be = Polygon(newcoords).exterior
-                    mls = be.intersection(be)
-                    if len(list(shapely.ops.polygonize(mls))) > 1:
-                        newgeom = MultiPolygon(shapely.ops.polygonize(mls))
-                        geoms = []
-                        for g in range(len(newgeom)):
-                            geoms.append(newgeom[g].area)
-                        newgeom = newgeom[geoms.index(max(geoms))]
-                    else:
-                        newgeom = list(shapely.ops.polygonize(mls))[0]
+
+def _queen_corners(tessellation, sensitivity, sindex):
+    """
+    Experimental: Fix unprecise corners.
+    """
+    changes = {}
+    qid = 0
+
+    for ix, row in tqdm(tessellation.iterrows(), total=tessellation.shape[0]):
+        corners = []
+        change = []
+
+        cell = row.geometry
+        coords = cell.exterior.coords
+        for i in coords:
+            point = Point(i)
+            possible_matches_index = list(sindex.intersection(point.bounds))
+            possible_matches = tessellation.iloc[possible_matches_index]
+            precise_matches = sum(possible_matches.intersects(point))
+            if precise_matches > 2:
+                corners.append(point)
+
+        if len(corners) > 2:
+            for c in range(len(corners)):
+                next_c = c + 1
+                if c == (len(corners) - 1):
+                    next_c = 0
+                if corners[c].distance(corners[next_c]) < sensitivity:
+                    change.append([corners[c], corners[next_c]])
+        elif len(corners) == 2:
+            if corners[0].distance(corners[1]) > 0:
+                if corners[0].distance(corners[1]) < sensitivity:
+                    change.append([corners[0], corners[1]])
+
+        if change:
+            for points in change:
+                x_new = np.mean([points[0].x, points[1].x])
+                y_new = np.mean([points[0].y, points[1].y])
+                new = [(x_new, y_new), id]
+                changes[(points[0].x, points[0].y)] = new
+                changes[(points[1].x, points[1].y)] = new
+                qid = qid + 1
+
+    for ix, row in tqdm(tessellation.iterrows(), total=tessellation.shape[0]):
+        cell = row.geometry
+        coords = list(cell.exterior.coords)
+
+        moves = {}
+        for x in coords:
+            if x in changes.keys():
+                moves[coords.index(x)] = changes[x]
+        keys = list(moves.keys())
+        delete_points = []
+        for move in range(len(keys)):
+            if move < len(keys) - 1:
+                if moves[keys[move]][1] == moves[keys[move + 1]][1] and keys[move + 1] - keys[move] < 5:
+                    delete_points = delete_points + (coords[keys[move]:keys[move + 1]])
+                    # change the code above to have if based on distance not number
+
+        newcoords = [changes[x][0] if x in changes.keys() else x for x in coords]
+        for coord in newcoords:
+            if coord in delete_points:
+                newcoords.remove(coord)
+        if coords != newcoords:
+            if not cell.interiors:
+                # newgeom = Polygon(newcoords).buffer(0)
+                be = Polygon(newcoords).exterior
+                mls = be.intersection(be)
+                if len(list(shapely.ops.polygonize(mls))) > 1:
+                    newgeom = MultiPolygon(shapely.ops.polygonize(mls))
+                    geoms = []
+                    for g in range(len(newgeom)):
+                        geoms.append(newgeom[g].area)
+                    newgeom = newgeom[geoms.index(max(geoms))]
                 else:
-                    newgeom = Polygon(newcoords, holes=cell.interiors)
-                morphological_tessellation.loc[ix, 'geometry'] = newgeom
-
-        # check against input layer
-        ids_original = list(gdf[unique_id])
-        ids_generated = list(morphological_tessellation[unique_id])
-        if len(ids_original) != len(ids_generated):
-            import warnings
-            diff = set(ids_original).difference(ids_generated)
-            warnings.warn("Tessellation does not fully match buildings. {len} element(s) collapsed "
-                          "during generation - unique_id: {i}".format(len=len(diff), i=diff))
-
-        # check MultiPolygons - usually caused by error in input geometry
-        uids = morphological_tessellation[morphological_tessellation.geometry.type == 'MultiPolygon'][unique_id]
-        if len(uids) > 0:
-            import warnings
-            warnings.warn('Tessellation contains MultiPolygon elements. Initial objects should be edited. '
-                          'unique_id of affected elements: {}'.format(list(uids)))
-
-    # translate back to true position
+                    newgeom = list(shapely.ops.polygonize(mls))[0]
+            else:
+                newgeom = Polygon(newcoords, holes=cell.interiors)
+            tessellation.loc[ix, 'geometry'] = newgeom
+    return tessellation
+
+
+def tessellation(gdf, unique_id, limit, shrink=0.4, segment=0.5, queen_corners=False, sensitivity=2):
+    """
+    Generate morphological tessellation around given buildings.
+
+    Parameters
+    ----------
+    gdf : GeoDataFrame
+        GeoDataFrame containing building footprints
+    unique_id : str
+        name of the column with unique id
+    limit : MultiPolygon or Polygon
+        MultiPolygon or Polygon defining the study area limiting tessellation (otherwise it could go to infinity).
+    shrink : float (default 0.4)
+        distance for negative buffer to generate space between adjacent buildings.
+    segment : float (default 0.5)
+        maximum distance between points on Polygon after discretisation
+
+    Returns
+    -------
+    GeoDataFrame
+        GeoDataFrame of morphological tessellation with the unique id based on original buildings.
+
+    Notes
+    -------
+    queen_corners and sensitivity are currently experimental only and can cause errors.
+    """
+    objects = gdf.copy()
+
+    centre = _get_centre(objects)
+    objects['geometry'] = objects['geometry'].translate(xoff=-centre[0], yoff=-centre[1])
+
+    print('Bufferring geometry...')
+    objects['geometry'] = objects.geometry.apply(lambda g: g.buffer(-shrink, cap_style=2, join_style=2))
+
+    print('Converting multipart geometry to singlepart...')
+    objects = objects.explode()
+    objects.reset_index(inplace=True, drop=True)
+
+    print('Densifying geometry...')
+    objects['geometry'] = objects['geometry'].apply(_densify, segment=segment)
+
+    print('Generating input point array...')
+    points, ids = _point_array(objects, unique_id)
+
+    # add convex hull buffered large distance to eliminate infinity issues
+    series = gpd.GeoSeries(limit, crs=gdf.crs).translate(xoff=-centre[0], yoff=-centre[1])
+    hull = series.geometry[0].convex_hull.buffer(300)
+    hull = _densify(hull, 20)
+    hull_array = np.array(hull.boundary.coords).tolist()
+    for i in range(len(hull_array)):
+        points.append(hull_array[i])
+        ids.append(-1)
+
+    print('Generating Voronoi diagram...')
+    voronoi_diagram = Voronoi(np.array(points))
+
+    print('Generating GeoDataFrame...')
+    regions_gdf = _regions(voronoi_diagram, unique_id, ids, crs=gdf.crs)
+
+    print('Dissolving Voronoi polygons...')
+    morphological_tessellation = regions_gdf[[unique_id, 'geometry']].dissolve(by=unique_id, as_index=False)
+
     morphological_tessellation['geometry'] = morphological_tessellation['geometry'].translate(xoff=centre[0], yoff=centre[1])
-    print('Tessellation finished.')
+
+    morphological_tessellation, sindex = _cut(morphological_tessellation, limit, unique_id)
+
+    if queen_corners is True:
+        morphological_tessellation = _queen_corners(morphological_tessellation, sensitivity, sindex)
+
+    _check_result(morphological_tessellation, gdf, unique_id=unique_id)
+
     return morphological_tessellation
 
 
diff --git a/momepy/utils.py b/momepy/utils.py
index 6825ed8..22e9a3e 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -241,6 +241,8 @@ def preprocess(buildings, size=30, compactness=True, islands=True):
         GeoDataFrame containing preprocessed geometry
     """
     blg = buildings.copy()
+    blg = blg.explode()
+    blg.reset_index(drop=True, inplace=True)
     for l in range(0, 2):
         print('Loop', l + 1, 'out of 2.')
         blg.reset_index(inplace=True, drop=True)


Test patch:
diff --git a/tests/test_elements.py b/tests/test_elements.py
index 01221ea..a391699 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -12,13 +12,14 @@ class TestElements:
         self.df_tessellation = gpd.read_file(test_file_path, layer='tessellation')
         self.df_streets = gpd.read_file(test_file_path, layer='streets')
         self.df_streets['nID'] = range(len(self.df_streets))
+        self.limit = mm.buffered_limit(self.df_buildings, 50)
 
     def test_tessellation(self):
-        tessellation = mm.tessellation(self.df_buildings)
+        tessellation = mm.tessellation(self.df_buildings, 'uID', self.limit, segment=2)
         assert len(tessellation) == len(self.df_tessellation)
-        queen_corners = mm.tessellation(self.df_buildings, queen_corners=True)
+        queen_corners = mm.tessellation(self.df_buildings, 'uID', self.limit, segment=2, queen_corners=True)
         w = libpysal.weights.Queen.from_dataframe(queen_corners)
-        assert w.neighbors[14] == [35, 36, 13, 15, 26, 27, 28, 30]
+        assert w.neighbors[14] == [35, 36, 13, 15, 26, 27, 28, 30, 31]
 
     def test_snap_street_network_edge(self):
         snapped = mm.snap_street_network_edge(self.df_streets, self.df_buildings, self.df_tessellation, 20, 70)