#!/usr/bin/env python import math import operator import warnings from copy import deepcopy import geopandas as gpd import libpysal import networkx as nx import numpy as np import pandas as pd import shapely from packaging.version import Version from scipy.signal import find_peaks from scipy.stats import gaussian_kde from shapely.geometry import LineString, Point from shapely.ops import linemerge, polygonize, split from tqdm.auto import tqdm from .coins import COINS from .graph import node_degree from .shape import CircularCompactness from .utils import gdf_to_nx, nx_to_gdf __all__ = [ "preprocess", "remove_false_nodes", "CheckTessellationInput", "close_gaps", "extend_lines", "roundabout_simplification", "consolidate_intersections", "FaceArtifacts", ] GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") def preprocess( buildings, size=30, compactness=0.2, islands=True, loops=2, verbose=True ): """ Preprocesses building geometry to eliminate additional structures being single features. Certain data providers (e.g. Ordnance Survey in GB) do not provide building geometry as one feature, but divided into different features depending their level (if they are on ground floor or not - passages, overhangs). Ideally, these features should share one building ID on which they could be dissolved. If this is not the case, series of steps needs to be done to minimize errors in morphological analysis. This script attempts to preprocess such geometry based on several condidions: If feature area is smaller than set size it will be a) deleted if it does not touch any other feature; b) will be joined to feature with which it shares the longest boundary. If feature is fully within other feature, these will be joined. If feature's circular compactness (:py:class:`momepy.CircularCompactness`) is < 0.2, it will be joined to feature with which it shares the longest boundary. Function does multiple loops through. Parameters ---------- buildings : geopandas.GeoDataFrame geopandas.GeoDataFrame containing building layer size : float (default 30) maximum area of feature to be considered as additional structure. Set to None if not wanted. compactness : float (default .2) if set, function will resolve additional structures identified based on their circular compactness. islands : bool (default True) if True, function will resolve additional structures which are fully within other structures (share 100% of exterior boundary). loops : int (default 2) number of loops verbose : bool (default True) if True, shows progress bars in loops and indication of steps Returns ------- GeoDataFrame GeoDataFrame containing preprocessed geometry """ blg = buildings.copy() blg = blg.explode(ignore_index=True) for loop in range(0, loops): print("Loop", loop + 1, f"out of {loops}.") if verbose else None blg.reset_index(inplace=True, drop=True) blg["mm_uid"] = range(len(blg)) sw = libpysal.weights.contiguity.Rook.from_dataframe( blg, silence_warnings=True, use_index=False ) blg["neighbors"] = sw.neighbors.values() blg["n_count"] = blg.apply(lambda row: len(row.neighbors), axis=1) blg["circu"] = CircularCompactness(blg).series # idetify those smaller than x with only one neighbor and attaches it to it. join = {} delete = [] for row in tqdm( blg.itertuples(), total=blg.shape[0], desc="Identifying changes", disable=not verbose, ): if size and row.geometry.area < size: if row.n_count == 1: uid = blg.iloc[row.neighbors[0]].mm_uid join.setdefault(uid, []).append(row.mm_uid) elif row.n_count > 1: shares = {} for n in row.neighbors: shares[n] = row.geometry.intersection( blg.at[n, blg.geometry.name] ).length maximal = max(shares.items(), key=operator.itemgetter(1))[0] uid = blg.loc[maximal].mm_uid join.setdefault(uid, []).append(row.mm_uid) else: delete.append(row.Index) if compactness and row.circu < compactness: if row.n_count == 1: uid = blg.iloc[row.neighbors[0]].mm_uid join.setdefault(uid, []).append(row.mm_uid) elif row.n_count > 1: shares = {} for n in row.neighbors: shares[n] = row.geometry.intersection( blg.at[n, blg.geometry.name] ).length maximal = max(shares.items(), key=operator.itemgetter(1))[0] uid = blg.loc[maximal].mm_uid join.setdefault(uid, []).append(row.mm_uid) if islands and row.n_count == 1: shared = row.geometry.intersection( blg.at[row.neighbors[0], blg.geometry.name] ).length if shared == row.geometry.exterior.length: uid = blg.iloc[row.neighbors[0]].mm_uid join.setdefault(uid, []).append(row.mm_uid) for key in tqdm( join, total=len(join), desc="Changing geometry", disable=not verbose ): selection = blg[blg["mm_uid"] == key] if not selection.empty: geoms = [selection.iloc[0].geometry] for j in join[key]: subset = blg[blg["mm_uid"] == j] if not subset.empty: geoms.append(blg[blg["mm_uid"] == j].iloc[0].geometry) blg.drop(blg[blg["mm_uid"] == j].index[0], inplace=True) new_geom = shapely.union_all(geoms) blg.loc[blg.loc[blg["mm_uid"] == key].index[0], blg.geometry.name] = ( new_geom ) blg.drop(delete, inplace=True) return blg[buildings.columns] def remove_false_nodes(gdf): """ Clean topology of existing LineString geometry by removal of nodes of degree 2. Returns the original gdf if there's no node of degree 2. Parameters ---------- gdf : GeoDataFrame, GeoSeries, array of shapely geometries (Multi)LineString data of street network Returns ------- gdf : GeoDataFrame, GeoSeries See also -------- momepy.extend_lines momepy.close_gaps """ if isinstance(gdf, gpd.GeoDataFrame | gpd.GeoSeries): # explode to avoid MultiLineStrings # reset index due to the bug in GeoPandas explode df = gdf.reset_index(drop=True).explode(ignore_index=True) # get underlying shapely geometry geom = df.geometry.array else: geom = gdf df = gpd.GeoSeries(gdf) # extract array of coordinates and number per geometry start_points = shapely.get_point(geom, 0) end_points = shapely.get_point(geom, -1) points = shapely.points( np.unique( shapely.get_coordinates(np.concatenate([start_points, end_points])), axis=0 ) ) # query LineString geometry to identify points intersecting 2 geometries tree = shapely.STRtree(geom) inp, res = tree.query(points, predicate="intersects") unique, counts = np.unique(inp, return_counts=True) mask = np.isin(inp, unique[counts == 2]) merge_res = res[mask] merge_inp = inp[mask] if len(merge_res): g = nx.Graph(list(zip(merge_inp * -1, merge_res, strict=True))) new_geoms = [] for c in nx.connected_components(g): valid = [ix for ix in c if ix > -1] new_geoms.append(shapely.line_merge(shapely.union_all(geom[valid]))) df = df.drop(merge_res) final = gpd.GeoSeries(new_geoms, crs=df.crs).explode(ignore_index=True) if isinstance(gdf, gpd.GeoDataFrame): combined = pd.concat( [ df, gpd.GeoDataFrame( {df.geometry.name: final}, geometry=df.geometry.name, crs=df.crs ), ], ignore_index=True, ) else: combined = pd.concat([df, final], ignore_index=True) # re-order closed loops fixed_loops = [] fixed_index = [] nodes = nx_to_gdf( node_degree( gdf_to_nx( combined if isinstance(combined, gpd.GeoDataFrame) else combined.to_frame("geometry") ) ), lines=False, ) loops = combined[combined.is_ring] node_ix, loop_ix = loops.sindex.query(nodes.geometry, predicate="intersects") for ix in np.unique(loop_ix): loop_geom = loops.geometry.iloc[ix] target_nodes = nodes.geometry.iloc[node_ix[loop_ix == ix]] if len(target_nodes) == 2: node_coords = shapely.get_coordinates(target_nodes) coords = np.array(loop_geom.coords) new_start = ( node_coords[0] if (node_coords[0] != coords[0]).all() else node_coords[1] ) new_start_idx = np.where(coords == new_start)[0][0] rolled_coords = np.roll(coords[:-1], -new_start_idx, axis=0) new_sequence = np.append(rolled_coords, rolled_coords[[0]], axis=0) fixed_loops.append(shapely.LineString(new_sequence)) fixed_index.append(ix) fixed_loops = gpd.GeoSeries(fixed_loops, crs=df.crs).explode(ignore_index=True) if isinstance(gdf, gpd.GeoDataFrame): return pd.concat( [ combined.drop(loops.iloc[fixed_index].index), gpd.GeoDataFrame( {df.geometry.name: fixed_loops}, geometry=df.geometry.name, crs=df.crs, ), ], ignore_index=True, ) else: return pd.concat( [combined.drop(loops.iloc[fixed_index].index), fixed_loops], ignore_index=True, ) # if there's nothing to fix, return the original dataframe return gdf class CheckTessellationInput: """ Check input data for :class:`Tessellation` for potential errors. :class:`Tessellation` requires data of relatively high level of precision and there are three particular patterns causing issues.\n 1. Features will collapse into empty polygon - these do not have tessellation cell in the end.\n 2. Features will split into MultiPolygon - at some cases, features with narrow links between parts split into two during 'shrinking'. In most cases that is not an issue and resulting tessellation is correct anyway, but sometimes this result in a cell being MultiPolygon, which is not correct.\n 3. Overlapping features - features which overlap even after 'shrinking' cause invalid tessellation geometry.\n :class:`CheckTessellationInput` will check for all of these. Overlapping features 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 tessellation-based analysis. Parameters ---------- gdf : GeoDataFrame or GeoSeries GeoDataFrame containing objects to be used as ``gdf`` in :class:`Tessellation` shrink : float (default 0.4) distance for negative buffer collapse : bool (default True) check for features which would collapse to empty polygon split : bool (default True) check for features which would split into Multi-type overlap : bool (default True) check for overlapping features (after negative buffer) Attributes ---------- collapse : GeoDataFrame or GeoSeries features which would collapse to empty polygon split : GeoDataFrame or GeoSeries features which would split into Multi-type overlap : GeoDataFrame or GeoSeries overlapping features (after negative buffer) Examples -------- >>> check = CheckTessellationData(df) Collapsed features : 3157 Split features : 519 Overlapping features: 22 """ warnings.filterwarnings("ignore", "GeoSeries.isna", UserWarning) def __init__(self, gdf, shrink=0.4, collapse=True, split=True, overlap=True): data = gdf[~gdf.is_empty] if split: types = data.geom_type shrink = data.buffer(-shrink) if shrink != 0 else data if collapse: emptycheck = shrink.is_empty self.collapse = gdf[emptycheck] collapsed = len(self.collapse) else: collapsed = "NA" if split: type_check = shrink.geom_type != types self.split = gdf[type_check] split_count = len(self.split) else: split_count = "NA" if overlap: shrink = shrink.reset_index(drop=True) shrink = shrink[~(shrink.is_empty | shrink.geometry.isna())] sindex = shrink.sindex hits = shrink.bounds.apply( lambda row: list(sindex.intersection(row)), axis=1 ) od_matrix = pd.DataFrame( { "origin": np.repeat(hits.index, hits.apply(len)), "dest": np.concatenate(hits.values), } ) od_matrix = od_matrix[od_matrix.origin != od_matrix.dest] duplicated = pd.DataFrame(np.sort(od_matrix, axis=1)).duplicated() od_matrix = od_matrix.reset_index(drop=True)[~duplicated] od_matrix = od_matrix.join( shrink.geometry.rename("o_geom"), on="origin" ).join(shrink.geometry.rename("d_geom"), on="dest") intersection = od_matrix.o_geom.values.intersection(od_matrix.d_geom.values) type_filter = gpd.GeoSeries(intersection).geom_type == "Polygon" empty_filter = intersection.is_empty overlapping = od_matrix.reset_index(drop=True)[empty_filter ^ type_filter] over_rows = sorted( pd.concat([overlapping.origin, overlapping.dest]).unique() ) self.overlap = gdf.iloc[over_rows] overlapping_c = len(self.overlap) else: overlapping_c = "NA" print( f"Collapsed features : {collapsed}\n" f"Split features : {split_count}\n" f"Overlapping features: {overlapping_c}" ) def close_gaps(gdf, tolerance): """Close gaps in LineString geometry where it should be contiguous. Snaps both lines to a centroid of a gap in between. Parameters ---------- gdf : GeoDataFrame, GeoSeries GeoDataFrame or GeoSeries containing LineString representation of a network. tolerance : float nodes within a tolerance will be snapped together Returns ------- GeoSeries See also -------- momepy.extend_lines momepy.remove_false_nodes """ geom = gdf.geometry.array coords = shapely.get_coordinates(geom) indices = shapely.get_num_coordinates(geom) # generate a list of start and end coordinates and create point geometries edges = [0] i = 0 for ind in indices: ix = i + ind edges.append(ix - 1) edges.append(ix) i = ix edges = edges[:-1] points = shapely.points(np.unique(coords[edges], axis=0)) buffered = shapely.buffer(points, tolerance / 2) dissolved = shapely.union_all(buffered) exploded = [ shapely.get_geometry(dissolved, i) for i in range(shapely.get_num_geometries(dissolved)) ] centroids = shapely.centroid(exploded) snapped = shapely.snap(geom, shapely.union_all(centroids), tolerance) return gpd.GeoSeries(snapped, crs=gdf.crs) def extend_lines(gdf, tolerance, target=None, barrier=None, extension=0): """Extends lines from gdf to itself 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 lines to itself. If ``barrier`` is passed, each extended line is checked for intersection 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 GeoDataFrame containing LineString geometry tolerance : float tolerance in snapping (by how much could be each segment extended). target : GeoDataFrame, GeoSeries 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 momepy.remove_false_nodes """ # explode to avoid MultiLineStrings # reset index due to the bug in GeoPandas explode df = gdf.reset_index(drop=True).explode(ignore_index=True) if target is None: target = df itself = True else: itself = False # get underlying shapely geometry geom = df.geometry.array # extract array of coordinates and number per geometry coords = shapely.get_coordinates(geom) indices = shapely.get_num_coordinates(geom) # generate a list of start and end coordinates and create point geometries edges = [0] i = 0 for ind in indices: ix = i + ind edges.append(ix - 1) edges.append(ix) i = ix edges = edges[:-1] points = shapely.points(np.unique(coords[edges], axis=0)) # query LineString geometry to identify points intersecting 2 geometries tree = shapely.STRtree(geom) inp, res = tree.query(points, predicate="intersects") unique, counts = np.unique(inp, return_counts=True) ends = np.unique(res[np.isin(inp, unique[counts == 1])]) new_geoms = [] # iterate over cul-de-sac-like segments and attempt to snap them to street network for line in ends: l_coords = shapely.get_coordinates(geom[line]) start = shapely.points(l_coords[0]) end = shapely.points(l_coords[-1]) first = list(tree.query(start, predicate="intersects")) second = list(tree.query(end, predicate="intersects")) first.remove(line) second.remove(line) t = target if not itself else target.drop(line) if first and not second: snapped = _extend_line(l_coords, t, tolerance) if ( barrier is not None and barrier.sindex.query( shapely.linestrings(snapped), predicate="intersects" ).size > 0 ): new_geoms.append(geom[line]) else: if extension == 0: new_geoms.append(shapely.linestrings(snapped)) else: new_geoms.append( shapely.linestrings( _extend_line(snapped, t, extension, snap=False) ) ) elif not first and second: snapped = _extend_line(np.flip(l_coords, axis=0), t, tolerance) if ( barrier is not None and barrier.sindex.query( shapely.linestrings(snapped), predicate="intersects" ).size > 0 ): new_geoms.append(geom[line]) else: if extension == 0: new_geoms.append(shapely.linestrings(snapped)) else: new_geoms.append( shapely.linestrings( _extend_line(snapped, t, extension, snap=False) ) ) elif not first and not second: one_side = _extend_line(l_coords, t, tolerance) one_side_e = _extend_line(one_side, t, extension, snap=False) snapped = _extend_line(np.flip(one_side_e, axis=0), t, tolerance) if ( barrier is not None and barrier.sindex.query( shapely.linestrings(snapped), predicate="intersects" ).size > 0 ): new_geoms.append(geom[line]) else: if extension == 0: new_geoms.append(shapely.linestrings(snapped)) else: new_geoms.append( shapely.linestrings( _extend_line(snapped, t, extension, snap=False) ) ) df.iloc[ends, df.columns.get_loc(df.geometry.name)] = new_geoms return df def _extend_line(coords, target, tolerance, snap=True): """ Extends a line geometry to snap on the target within a tolerance. """ if snap: extrapolation = _get_extrapolated_line( coords[-4:] if len(coords.shape) == 1 else coords[-2:].flatten(), tolerance, ) int_idx = target.sindex.query(extrapolation, predicate="intersects") intersection = shapely.intersection( target.iloc[int_idx].geometry.array, extrapolation ) if intersection.size > 0: if len(intersection) > 1: distances = {} ix = 0 for p in intersection: distance = shapely.distance(p, shapely.points(coords[-1])) distances[ix] = distance ix = ix + 1 minimal = min(distances.items(), key=operator.itemgetter(1))[0] new_point_coords = shapely.get_coordinates(intersection[minimal]) else: new_point_coords = shapely.get_coordinates(intersection[0]) coo = np.append(coords, new_point_coords) new = np.reshape(coo, (len(coo) // 2, 2)) return new return coords extrapolation = _get_extrapolated_line( coords[-4:] if len(coords.shape) == 1 else coords[-2:].flatten(), tolerance, point=True, ) return np.vstack([coords, extrapolation]) def _get_extrapolated_line(coords, tolerance, point=False): """ Creates a shapely line extrapoled in p1->p2 direction. """ p1 = coords[:2] p2 = coords[2:] a = p2 # defining new point based on the vector between existing points if p1[0] >= p2[0] and p1[1] >= p2[1]: b = ( p2[0] - tolerance * math.cos( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), p2[1] - tolerance * math.sin( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), ) elif p1[0] <= p2[0] and p1[1] >= p2[1]: b = ( p2[0] + tolerance * math.cos( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), p2[1] - tolerance * math.sin( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), ) elif p1[0] <= p2[0] and p1[1] <= p2[1]: b = ( p2[0] + tolerance * math.cos( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), p2[1] + tolerance * math.sin( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), ) else: b = ( p2[0] - tolerance * math.cos( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), p2[1] + tolerance * math.sin( math.atan( math.fabs(p1[1] - p2[1] + 0.000001) / math.fabs(p1[0] - p2[0] + 0.000001) ) ), ) if point: return b return shapely.linestrings([a, b]) def _polygonize_ifnone(edges, polys): if polys is None: pre_polys = polygonize(edges.geometry) polys = gpd.GeoDataFrame(geometry=list(pre_polys), crs=edges.crs) return polys def _selecting_rabs_from_poly( gdf, area_col="area", circom_threshold=0.7, area_threshold=0.85, include_adjacent=True, diameter_factor=1.5, ): """ From a GeoDataFrame of polygons, returns a GDF of polygons that are above the Circular Compactness threshold. Return ________ GeoDataFrames : roundabouts and adjacent polygons """ # calculate parameters if area_col == "area": gdf.loc[:, area_col] = gdf.geometry.area circom_serie = CircularCompactness(gdf, area_col).series # selecting roundabout polygons based on compactness mask = circom_serie > circom_threshold rab = gdf[mask] # exclude those above the area threshold area_threshold_val = gdf.area.quantile(area_threshold) rab = rab[rab[area_col] < area_threshold_val] if include_adjacent is True: bounds = rab.geometry.bounds rab = pd.concat([rab, bounds], axis=1) rab["deltax"] = rab.maxx - rab.minx rab["deltay"] = rab.maxy - rab.miny rab["rab_diameter"] = rab[["deltax", "deltay"]].max(axis=1) # selecting the adjacent areas that are of smaller than itself rab_adj = gpd.sjoin(gdf, rab, predicate="intersects") area_right = area_col + "_right" area_left = area_col + "_left" area_mask = rab_adj[area_right] >= rab_adj[area_left] rab_adj = rab_adj[area_mask] rab_adj.index.name = "index" # adding a hausdorff_distance threshold rab_adj["hdist"] = 0.0 # TODO: (should be a way to vectorize) for i, group in rab_adj.groupby("index_right"): for g in group.itertuples(): hdist = g.geometry.hausdorff_distance(rab.loc[i].geometry) rab_adj.loc[g.Index, "hdist"] = hdist rab_plus = rab_adj[rab_adj.hdist < (rab_adj.rab_diameter * diameter_factor)] else: rab["index_right"] = rab.index rab_plus = rab # only keeping relevant fields geom_col = rab_plus.geometry.name rab_plus = rab_plus[[geom_col, "index_right"]] return rab_plus def _rabs_center_points(gdf, center_type="centroid"): """ From a selection of roundabouts, returns an aggregated GeoDataFrame per roundabout with extra column with center_type. """ # temporary DataFrame where geometry is the array of shapely geometries # Hack until shapely 2.0 is out. # TODO: replace shapely with shapely 2.0 tmp = pd.DataFrame(gdf.copy()) # creating a copy avoids warnings tmp["geometry"] = tmp.geometry.array shapely_geoms = ( tmp.groupby("index_right") .geometry.apply(shapely.multipolygons) .rename("geometry") ) shapely_geoms = shapely.make_valid(shapely_geoms) rab_multipolygons = gpd.GeoDataFrame(shapely_geoms, crs=gdf.crs) # make_valid is transforming the multipolygons into geometry collections because of # shared edges if center_type == "centroid": # geometry centroid of the actual circle rab_multipolygons["center_pt"] = gdf[ gdf.index == gdf.index_right ].geometry.centroid elif center_type == "mean": coords, idxs = shapely.get_coordinates(shapely_geoms, return_index=True) means = {} for i in np.unique(idxs): tmps = coords[idxs == i] target_idx = rab_multipolygons.index[i] means[target_idx] = Point(tmps.mean(axis=0)) rab_multipolygons["center_pt"] = gpd.GeoSeries(means, crs=gdf.crs) # centerpoint of minimum_bounding_circle # TODO # minimun_bounding_circle() should be available in Shapely 2.0. return rab_multipolygons def _coins_filtering_many_incoming(incoming_many, angle_threshold=0): """ Used only for the cases when more than one incoming line touches the roundabout. """ idx_out_many_incoming = [] # For each new connection, evaluate COINS and select the group from which the new # line belongs # TODO ideally use the groupby object on line_wkt used earlier for _g, x in incoming_many.groupby("line_wkt"): gs = gpd.GeoSeries(pd.concat([x.geometry, x.line]), crs=incoming_many.crs) gdf = gpd.GeoDataFrame(geometry=gs) gdf = gdf.drop_duplicates() coins = COINS(gdf, angle_threshold=angle_threshold) group_series = coins.stroke_attribute() gdf["coins_group"] = group_series # selecting the incoming and its extension coins_group_filter = gdf.groupby("coins_group").count() == 1 f = gdf.coins_group.map(coins_group_filter.geometry) idxs_remove = gdf[f].index idx_out_many_incoming.extend(idxs_remove) incoming_many_reduced = incoming_many.drop(idx_out_many_incoming, axis=0) return incoming_many_reduced def _selecting_incoming_lines(rab_multipolygons, edges, angle_threshold=0): """Selecting only the lines that are touching but not covered by the ``rab_plus``. If more than one LineString is incoming to ``rab_plus``, COINS algorithm is used to select the line to be extended further. """ # selecting the lines that are touching but not covered by touching = gpd.sjoin(edges, rab_multipolygons, predicate="touches") if GPD_GE_013: edges_idx, _ = rab_multipolygons.sindex.query( edges.geometry, predicate="covered_by" ) else: edges_idx, _ = rab_multipolygons.sindex.query_bulk( edges.geometry, predicate="covered_by" ) idx_drop = edges.index.take(edges_idx) touching_idx = touching.index ls = list(set(touching_idx) - set(idx_drop)) incoming = touching.loc[ls] # figuring out which ends of incoming edges need to be connected to the center_pt incoming["first_pt"] = incoming.geometry.apply(lambda x: Point(x.coords[0])) incoming["dist_first_pt"] = incoming.center_pt.distance(incoming.first_pt) incoming["last_pt"] = incoming.geometry.apply(lambda x: Point(x.coords[-1])) incoming["dist_last_pt"] = incoming.center_pt.distance(incoming.last_pt) lines = [] for _i, row in incoming.iterrows(): if row.dist_first_pt < row.dist_last_pt: lines.append(LineString([row.first_pt, row.center_pt])) else: lines.append(LineString([row.last_pt, row.center_pt])) incoming["line"] = gpd.GeoSeries(lines, index=incoming.index, crs=edges.crs) # checking if there are more than one incoming lines arriving to the same point # which would create several new lines incoming["line_wkt"] = incoming.line.to_wkt() grouped_lines = incoming.groupby(["line_wkt"])["line_wkt"] count_s = grouped_lines.count() # separating the incoming roads that come on their own to those that come in groups filter_count_one = pd.DataFrame(count_s[count_s == 1]) filter_count_many = pd.DataFrame(count_s[count_s > 1]) incoming_ones = pd.merge( incoming, filter_count_one, left_on="line_wkt", right_index=True, how="inner" ) incoming_many = pd.merge( incoming, filter_count_many, left_on="line_wkt", right_index=True, how="inner" ) incoming_many_reduced = _coins_filtering_many_incoming( incoming_many, angle_threshold=angle_threshold ) incoming_all = gpd.GeoDataFrame( pd.concat([incoming_ones, incoming_many_reduced]), crs=edges.crs ) return incoming_all, idx_drop def _ext_lines_to_center(edges, incoming_all, idx_out): """ Extends the LineStrings geometries to the centerpoint defined by _rabs_center_points. Also deletes the lines that originally defined the roundabout. Creates a new column labled with the 'rab' number. Returns ------- GeoDataFrame GeoDataFrame with updated geometry """ incoming_all["geometry"] = incoming_all.apply( lambda row: linemerge([row.geometry, row.line]), axis=1 ) new_edges = edges.drop(idx_out, axis=0) # creating a unique group label for returned gdf _, inv = np.unique(incoming_all.index_right, return_inverse=True) incoming_label = pd.Series(inv, index=incoming_all.index) incoming_label = incoming_label[~incoming_label.index.duplicated(keep="first")] # maintaining the same gdf shape as the original incoming_all = incoming_all[edges.columns] new_edges = pd.concat([new_edges, incoming_all]) # adding a new column to match new_edges["simplification_group"] = incoming_label.astype("Int64") return new_edges def roundabout_simplification( edges, polys=None, area_col="area", circom_threshold=0.7, area_threshold=0.85, include_adjacent=True, diameter_factor=1.5, center_type="centroid", angle_threshold=0, ): """ Selects the roundabouts from ``polys`` to create a center point to merge all incoming edges. If None is passed, the function will perform shapely polygonization. All ``edges`` attributes are preserved and roundabouts are deleted. Note that some attributes, like length, may no longer reflect the reality of newly constructed geometry. If ``include_adjacent`` is True, adjacent polygons to the actual roundabout are also selected for simplification if two conditions are met: - the area of adjacent polygons is less than the actual roundabout - adjacent polygons do not extend beyond a factor of the diameter of the actual roundabout. This uses hausdorff_distance algorithm. Parameters ---------- edges : GeoDataFrame GeoDataFrame containing LineString geometry of urban network polys : GeoDataFrame GeoDataFrame containing Polygon geometry derived from polygonyzing ``edges`` GeoDataFrame. area_col : string Column name containing area values if ``polys`` GeoDataFrame contains such information. Otherwise, it will circom_threshold : float (default 0.7) Circular compactness threshold to select roundabouts from ``polys`` GeoDataFrame. Polygons with a higher or equal threshold value will be considered for simplification. area_threshold : float (default 0.85) Percentile threshold value from the area of ``polys`` to leave as input geometry. Polygons with a higher or equal threshold will be considered as urban blocks not considered for simplification. include_adjacent : boolean (default True) Adjacent polygons to be considered also as part of the simplification. diameter_factor : float (default 1.5) The factor to be applied to the diameter of each roundabout that determines how far an adjacent polygon can stretch until it is no longer considered part of the overall roundabout group. Only applyies when include_adjacent = True. center_type : string (default 'centroid') Method to use for converging the incoming LineStrings. Current list of options available : 'centroid', 'mean'. - 'centroid': selects the centroid of the actual roundabout (ignoring adjacent geometries) - 'mean': calculates the mean coordinates from the points of polygons (including adjacent geometries) angle_threshold : int, float (default 0) The angle threshold for the COINS algorithm. Only used when multiple incoming LineStrings arrive at the same Point to the roundabout or to the adjacent polygons if set as True. eg. when two 'edges' touch the roundabout at the same point, COINS algorithm will evaluate which of those incoming lines should be extended according to their deflection angle. Segments will only be considered a part of the same street if the deflection angle is above the threshold. Returns ------- GeoDataFrame GeoDataFrame with an updated geometry and an additional column labeling modified edges. """ if len(edges[edges.geom_type != "LineString"]) > 0: raise TypeError( "Only LineString geometries are allowed. " "Try using the `explode()` method to explode MultiLineStrings." ) polys = _polygonize_ifnone(edges, polys) rab = _selecting_rabs_from_poly( polys, area_col=area_col, circom_threshold=circom_threshold, area_threshold=area_threshold, include_adjacent=include_adjacent, diameter_factor=diameter_factor, ) rab_multipolygons = _rabs_center_points(rab, center_type=center_type) incoming_all, idx_drop = _selecting_incoming_lines( rab_multipolygons, edges, angle_threshold=angle_threshold ) output = _ext_lines_to_center(edges, incoming_all, idx_drop) return output def consolidate_intersections( graph, tolerance=30, rebuild_graph=True, rebuild_edges_method="spider", x_att="x", y_att="y", edge_from_att="from", edge_to_att="to", ): """ Consolidate close street intersections into a single node, collapsing short edges. If rebuild_graph is True, new edges are drawn according to ``rebuild_edges_method`` which is one of: 1. Extension reconstruction: Edges are linearly extended from original endpoints until the new nodes. This method preserves most faithfully the network geometry but can result in overlapping geometry. 2. Spider-web reconstruction: Edges are cropped within a buffer of the new endpoints and linearly extended from there. This method improves upon linear reconstruction by mantaining, when possible, network planarity. 3. Euclidean reconstruction: Edges are ignored and new edges are built as straight lines between new origin and new destination. This method ignores geometry, but efficiently preserves adjacency. If ``rebuild_graph`` is False, graph is returned with consolidated nodes but without reconstructed edges i.e. graph is intentionally disconnected. Graph must be configured so that 1. All nodes have attributes determining their x and y coordinates; 2. All edges have attributes determining their origin, destination, and geometry. Parameters ---------- graph : Networkx.Graph, Networkx.DiGraph, Networkx.MultiGraph, or Networkx.MultiDiGraph tolerance : float, default 30 distance in network units below which nodes will be consolidated rebuild_graph : bool rebuild_edges_method : str 'extend' or 'spider' or 'euclidean', ignored if rebuild_graph is False x_att : str node attribute with the valid x-coordinate y_att : str node attribute with the valid y-coordinate edge_from_att : str edge attribute with the valid origin node id edge_to_att : str edge attribute with the valid destination node id Returns ------- Networkx.MultiGraph or Networkx.MultiDiGraph directionality inferred from input type """ # Collect nodes and their data: nodes, nodes_dict = zip(*graph.nodes(data=True), strict=False) nodes_df = pd.DataFrame(nodes_dict, index=nodes) graph_crs = graph.graph.get("crs") # Create a graph without the edges above a certain length and clean it # from isolated nodes (the unsimplifiable nodes): components_graph = deepcopy(graph) components_graph.remove_edges_from( [ edge for edge in graph.edges(keys=True, data=True) if edge[-1]["length"] > tolerance ] ) isolated_nodes_list = list(nx.isolates(components_graph)) components_graph.remove_nodes_from(isolated_nodes_list) # The connected components of this graph are node clusters we must individually # simplify. We collect them in a dataframe and retrieve node properties (x, y # coords mainly) from the original graph. components = nx.connected_components(components_graph) components_dict = dict(enumerate(components, start=max(nodes) + 1)) nodes_to_merge_dict = { node: cpt for cpt, nodes in components_dict.items() for node in nodes } new_nodes_df = pd.DataFrame.from_dict( nodes_to_merge_dict, orient="index", columns=["cluster"] ) nodes_to_merge_df = pd.concat( [new_nodes_df, nodes_df[[x_att, y_att]]], axis=1, join="inner" ) # The two node attributes we need for the clusters are the position of the cluster # centroids. Those are obtained by averaging the x and y columns. We also add # . attribtues referring to the original node ids in every cluster: cluster_centroids_df = nodes_to_merge_df.groupby("cluster").mean() cluster_centroids_df["simplified"] = True cluster_centroids_df["original_node_ids"] = cluster_centroids_df.index.map( components_dict ) cluster_geometries = gpd.points_from_xy( cluster_centroids_df[x_att], cluster_centroids_df[y_att] ) cluster_gdf = gpd.GeoDataFrame( cluster_centroids_df, crs=graph_crs, geometry=cluster_geometries ) cluster_nodes_list = list(cluster_gdf.to_dict("index").items()) # Create a simplified graph object: simplified_graph = graph.copy() # Rebuild edges if necessary: if rebuild_graph: rebuild_edges_method = rebuild_edges_method.lower() simplified_graph.graph["approach"] = "primal" edges_gdf = nx_to_gdf(simplified_graph, points=False, lines=True) simplified_edges = _get_rebuilt_edges( edges_gdf, nodes_to_merge_dict, cluster_gdf, method=rebuild_edges_method, buffer=1.5 * tolerance, edge_from_att=edge_from_att, edge_to_att=edge_to_att, ) # Replacing the collapsed nodes with centroids and adding edges: simplified_graph.remove_nodes_from(nodes_to_merge_df.index) simplified_graph.add_nodes_from(cluster_nodes_list) if rebuild_graph: simplified_graph.add_edges_from(simplified_edges) return simplified_graph def _get_rebuilt_edges( edges_gdf, nodes_dict, cluster_gdf, method="spider", buffer=45, edge_from_att="from", edge_to_att="to", ): """ Update origin and destination on network edges when original endpoints were replaced by a consolidated node cluster. New edges are drawn according to method which is one of: 1. Extension reconstruction: Edges are linearly extended from original endpoints until the new nodes. This method preserves most faithfully the network geometry. 2. Spider-web reconstruction: Edges are cropped within a buffer of the new endpoints and linearly extended from there. This method improves upon linear reconstruction by mantaining, when possible, network planarity. 3. Euclidean reconstruction: Edges are ignored and new edges are built as straightlines between new origin and new destination. This method ignores geometry, but efficiently preserves adjacency. Parameters ---------- edges_gdf : GeoDataFrame GeoDataFrame containing LineString geometry and columns determining origin and destination node ids nodes_dict : dict Dictionary whose keys are node ids and values are the corresponding consolidated node cluster ids. Only consolidated nodes are in the dictionary. cluster_gdf : GeoDataFrame GeoDataFrame containing consolidated node ids. method: string 'extension' or 'spider' or 'euclidean' buffer : float distance to buffer consolidated nodes in the Spider-web reconstruction edge_from_att : str edge attribute with the valid origin node id edge_to_att : str edge attribute with the valid destination node id Returns ---------- List list of edges that should be added to the network. Edges are in the format (origin_id, destination_id, data), where data is inferred from edges_gdf """ # Determine what endpoints were made into clusters: edges_gdf["origin_cluster"] = edges_gdf[edge_from_att].apply( lambda u: nodes_dict.get(u, -1) ) edges_gdf["destination_cluster"] = edges_gdf[edge_to_att].apply( lambda v: nodes_dict.get(v, -1) ) # Determine what edges need to be simplified (either between diff. # clusters or self-loops in a cluster): edges_tosimplify_gdf = edges_gdf.query( "origin_cluster != destination_cluster or " f"(('{edge_to_att}' == '{edge_from_att}') and origin_cluster >= 0)" ) # Determine the new point geometries (when exists): edges_tosimplify_gdf = edges_tosimplify_gdf.assign( new_origin_pt=edges_tosimplify_gdf.origin_cluster.map( cluster_gdf.geometry, None ) ) edges_tosimplify_gdf = edges_tosimplify_gdf.assign( new_destination_pt=edges_tosimplify_gdf.destination_cluster.map( cluster_gdf.geometry, None ) ) # Determine the new geometry according to the simplification method: if method == "extend": edges_simplified_geometries = edges_tosimplify_gdf.apply( lambda edge: _extension_simplification( edge.geometry, edge.new_origin_pt, edge.new_destination_pt ), axis=1, ) edges_simplified_gdf = edges_tosimplify_gdf.assign( new_geometry=edges_simplified_geometries ) elif method == "euclidean": edges_simplified_geometries = edges_tosimplify_gdf.apply( lambda edge: _euclidean_simplification( edge.geometry, edge.new_origin_pt, edge.new_destination_pt ), axis=1, ) edges_simplified_gdf = edges_tosimplify_gdf.assign( new_geometry=edges_simplified_geometries ) elif method == "spider": edges_simplified_geometries = edges_tosimplify_gdf.apply( lambda edge: _spider_simplification( edge.geometry, edge.new_origin_pt, edge.new_destination_pt, buffer ), axis=1, ) edges_simplified_gdf = edges_tosimplify_gdf.assign( new_geometry=edges_simplified_geometries ) else: msg = ( f"Simplification '{method}' not recognized. See documentation for options." ) raise ValueError(msg) # Rename and update the columns: cols_rename = { edge_from_att: "original_from", edge_to_att: "original_to", "origin_cluster": edge_from_att, "destination_cluster": edge_to_att, "geometry": "original_geometry", } new_edges_gdf = edges_simplified_gdf.rename(cols_rename, axis=1) cols_drop = ["new_origin_pt", "new_destination_pt"] new_edges_gdf = new_edges_gdf.drop(columns=cols_drop) new_edges_gdf = new_edges_gdf.set_geometry("new_geometry") new_edges_gdf.loc[:, "length"] = new_edges_gdf.length # Update the indices: new_edges_gdf.loc[:, edge_from_att] = new_edges_gdf[edge_from_att].where( new_edges_gdf[edge_from_att] >= 0, new_edges_gdf["original_from"] ) new_edges_gdf.loc[:, edge_to_att] = new_edges_gdf[edge_to_att].where( new_edges_gdf[edge_to_att] >= 0, new_edges_gdf["original_to"] ) # Get the edge list with (from, to, data): new_edges_list = list( zip( new_edges_gdf[edge_from_att], new_edges_gdf[edge_to_att], new_edges_gdf.iloc[:, 2:].to_dict("index").values(), strict=False, ) ) return new_edges_list def _extension_simplification(geometry, new_origin, new_destination): """ Extends edge geometry to new endpoints. If either new_origin or new_destination is None, maintains the respective current endpoint. Parameters ---------- geometry : shapely.LineString new_origin : shapely.Point or None new_destination: shapely.Point or None Returns ---------- shapely.LineString """ # If we are dealing with a self-loop the line has no endpoints: if new_origin == new_destination: current_node = Point(geometry.coords[0]) geometry = linemerge([LineString([new_origin, current_node]), geometry]) # Assuming the line is not closed, we can find its endpoints: else: current_origin, current_destination = geometry.boundary.geoms if new_origin is not None: geometry = linemerge([LineString([new_origin, current_origin]), geometry]) if new_destination is not None: geometry = linemerge( [geometry, LineString([current_destination, new_destination])] ) return geometry def _spider_simplification(geometry, new_origin, new_destination, buff=15): """ Extends edge geometry to new endpoints via a "spider-web" method. Breaks current geometry within a buffer of the new endpoint and then extends it linearly. Useful to maintain planarity. If either new_origin or new_destination is None, maintains the respective current endpoint. Parameters ---------- geometry : shapely.LineString new_origin : shapely.Point or None new_destination: shapely.Point or None buff : float distance from new endpoint to break current geometry Returns ---------- shapely.LineString """ # If we are dealing with a self-loop the line has no boundary # . and we just use the first coordinate: if new_origin == new_destination: current_node = Point(geometry.coords[0]) geometry = linemerge([LineString([new_origin, current_node]), geometry]) # Assuming the line is not closed, we can find its endpoints # via the boundary attribute: else: current_origin, current_destination = geometry.boundary.geoms if new_origin is not None: # Create a buffer around the new origin: new_origin_buffer = new_origin.buffer(buff) # Use shapely.ops.split to break the edge where it # intersects the buffer: geometry_split_by_buffer_list = list( split(geometry, new_origin_buffer).geoms ) # If only one geometry results, edge does not intersect # buffer and line should connect new origin to old origin if len(geometry_split_by_buffer_list) == 1: geometry_split_by_buffer = geometry_split_by_buffer_list[0] splitting_point = current_origin # If more than one geometry, merge all linestrings # but the first and get their origin else: geometry_split_by_buffer = linemerge(geometry_split_by_buffer_list[1:]) splitting_point = geometry_split_by_buffer.boundary.geoms[0] # Merge this into new geometry: additional_line = [LineString([new_origin, splitting_point])] # Consider MultiLineStrings separately: if geometry_split_by_buffer.geom_type == "MultiLineString": geometry = linemerge( additional_line + list(geometry_split_by_buffer.geoms) ) else: geometry = linemerge(additional_line + [geometry_split_by_buffer]) if new_destination is not None: # Create a buffer around the new destination: new_destination_buffer = new_destination.buffer(buff) # Use shapely.ops.split to break the edge where it # intersects the buffer: geometry_split_by_buffer_list = list( split(geometry, new_destination_buffer).geoms ) # If only one geometry results, edge does not intersect # . buffer and line should connect new destination to old destination if len(geometry_split_by_buffer_list) == 1: geometry_split_by_buffer = geometry_split_by_buffer_list[0] splitting_point = current_destination # If more than one geometry, merge all linestrings # but the last and get their destination else: geometry_split_by_buffer = linemerge(geometry_split_by_buffer_list[:-1]) splitting_point = geometry_split_by_buffer.boundary.geoms[1] # Merge this into new geometry: additional_line = [LineString([splitting_point, new_destination])] # Consider MultiLineStrings separately: if geometry_split_by_buffer.geom_type == "MultiLineString": geometry = linemerge( list(geometry_split_by_buffer.geoms) + additional_line ) else: geometry = linemerge([geometry_split_by_buffer] + additional_line) return geometry def _euclidean_simplification(geometry, new_origin, new_destination): """ Rebuilds edge geometry to new endpoints. Ignores current geometry and traces a straight line between new endpoints. If either new_origin or new_destination is None, maintains the respective current endpoint. Parameters ---------- geometry : shapely.LineString new_origin : shapely.Point or None new_destination : shapely.Point or None Returns ---------- shapely.LineString """ # If we are dealing with a self-loop, geometry will be null! if new_origin == new_destination: geometry = None # Assuming the line is not closed, we can find its endpoints: else: current_origin, current_destination = geometry.boundary.geoms if new_origin is not None: if new_destination is not None: geometry = LineString([new_origin, new_destination]) else: geometry = LineString([new_origin, current_destination]) else: if new_destination is not None: geometry = LineString([current_origin, new_destination]) return geometry class FaceArtifacts: """Identify face artifacts in street networks For a given street network composed of transportation-oriented geometry containing features representing things like roundabouts, dual carriegaways and complex intersections, identify areas enclosed by geometry that is considered a `face artifact` as per :cite:`fleischmann2023`. Face artifacts highlight areas with a high likelihood of being of non-morphological (e.g. transporation) origin and may require simplification prior morphological analysis. See :cite:`fleischmann2023` for more details. Parameters ---------- gdf : geopandas.GeoDataFrame GeoDataFrame containing street network represented as (Multi)LineString geometry index : str, optional A type of the shape compacntess index to be used. Available are ['circlular_compactness', 'isoperimetric_quotient', 'diameter_ratio'], by default "circular_compactness" height_mins : float, optional Required depth of valleys, by default -np.inf height_maxs : float, optional Required height of peaks, by default 0.008 prominence : float, optional Required prominence of peaks, by default 0.00075 Attributes ---------- threshold : float Identified threshold between face polygons and face artifacts face_artifacts : GeoDataFrame A GeoDataFrame of geometries identified as face artifacts polygons : GeoDataFrame All polygons resulting from polygonization of the input gdf with the face_artifact_index kde : scipy.stats._kde.gaussian_kde Representation of a kernel-density estimate using Gaussian kernels. pdf : numpy.ndarray Probability density function peaks : numpy.ndarray locations of peaks in pdf valleys : numpy.ndarray locations of valleys in pdf Examples -------- >>> fa = momepy.FaceArtifacts(street_network_prague) >>> fa.threshold 6.9634555986177045 >>> fa.face_artifacts.head() geometry face_artifact_index 6 POLYGON ((-744164.625 -1043922.362, -744167.39... 5.112844 9 POLYGON ((-744154.119 -1043804.734, -744152.07... 6.295660 10 POLYGON ((-744101.275 -1043738.053, -744103.80... 2.862871 12 POLYGON ((-744095.511 -1043623.478, -744095.35... 3.712403 17 POLYGON ((-744488.466 -1044533.317, -744489.33... 5.158554 """ def __init__( self, gdf, index="circular_compactness", height_mins=-np.inf, height_maxs=0.008, prominence=0.00075, ): try: from esda import shape except (ImportError, ModuleNotFoundError) as err: raise ImportError( "The `esda` package is required. You can install it using " "`pip install esda` or `conda install esda -c conda-forge`." ) from err # Polygonize street network polygons = gpd.GeoSeries( shapely.polygonize( # polygonize [gdf.dissolve().geometry.item()] ) ).explode(ignore_index=True) # Store geometries as a GeoDataFrame self.polygons = gpd.GeoDataFrame(geometry=polygons) if index == "circular_compactness": self.polygons["face_artifact_index"] = np.log( shape.minimum_bounding_circle_ratio(polygons) * polygons.area ) elif index == "isoperimetric_quotient": self.polygons["face_artifact_index"] = np.log( shape.isoperimetric_quotient(polygons) * polygons.area ) elif index == "diameter_ratio": self.polygons["face_artifact_index"] = np.log( shape.diameter_ratio(polygons) * polygons.area ) else: raise ValueError( f"'{index}' is not supported. Use one of ['circlular_compactness', " "'isoperimetric_quotient', 'diameter_ratio']" ) # parameters for peak/valley finding peak_parameters = { "height_mins": height_mins, "height_maxs": height_maxs, "prominence": prominence, } mylinspace = np.linspace( self.polygons["face_artifact_index"].min(), self.polygons["face_artifact_index"].max(), 1000, ) self.kde = gaussian_kde( self.polygons["face_artifact_index"], bw_method="silverman" ) self.pdf = self.kde.pdf(mylinspace) # find peaks self.peaks, self.d_peaks = find_peaks( x=self.pdf, height=peak_parameters["height_maxs"], threshold=None, distance=None, prominence=peak_parameters["prominence"], width=1, plateau_size=None, ) # find valleys self.valleys, self.d_valleys = find_peaks( x=-self.pdf + 1, height=peak_parameters["height_mins"], threshold=None, distance=None, prominence=peak_parameters["prominence"], width=1, plateau_size=None, ) # check if we have at least 2 peaks condition_2peaks = len(self.peaks) > 1 # check if we have at least 1 valley condition_1valley = len(self.valleys) > 0 conditions = [condition_2peaks, condition_1valley] # if both these conditions are true, we find the artifact index if all(conditions): # find list order of highest peak highest_peak_listindex = np.argmax(self.d_peaks["peak_heights"]) # find index (in linspace) of highest peak highest_peak_index = self.peaks[highest_peak_listindex] # define all possible peak ranges fitting our definition peak_bounds = list(zip(self.peaks[:-1], self.peaks[1:], strict=True)) peak_bounds_accepted = [b for b in peak_bounds if highest_peak_index in b] # find all valleys that lie between two peaks valleys_accepted = [ v_index for v_index in self.valleys if any(v_index in range(b[0], b[1]) for b in peak_bounds_accepted) ] # the value of the first of those valleys is our artifact index # get the order of the valley valley_index = valleys_accepted[0] # derive threshold value for given option from index/linspace self.threshold = float(mylinspace[valley_index]) self.face_artifacts = self.polygons[ self.polygons.face_artifact_index < self.threshold ] else: warnings.warn( "No threshold found. Either your dataset it too small or the " "distribution of the face artifact index does not follow the " "expected shape.", UserWarning, stacklevel=2, ) self.threshold = None self.face_artifacts = None