From a225ba8c5b586c713e602ed5fac26c0e169ffd8b Mon Sep 17 00:00:00 2001 From: Christian Heider Nielsen Date: Fri, 8 Nov 2024 10:54:19 +0100 Subject: [PATCH] desliver generalised --- jord/geometric_analysis/center_line.py | 50 ++--- .../serialisation/well_known_binary.py | 13 +- .../serialisation/well_known_text.py | 7 +- jord/shapely_utilities/__init__.py | 1 + jord/shapely_utilities/lines.py | 100 +++++++-- jord/shapely_utilities/morphology.py | 2 + jord/shapely_utilities/points.py | 7 +- jord/shapely_utilities/polygons.py | 2 +- jord/shapely_utilities/projection.py | 7 +- jord/shapely_utilities/rings.py | 96 ++++++++- jord/shapely_utilities/sliverin.py | 191 +++++++++++++++++- jord/shapely_utilities/subdivision.py | 76 +++++++ tests/shapely/test_polygon_utilities.py | 2 +- tests/shapely/test_sliverin.py | 47 ++++- 14 files changed, 534 insertions(+), 67 deletions(-) create mode 100644 jord/shapely_utilities/subdivision.py diff --git a/jord/geometric_analysis/center_line.py b/jord/geometric_analysis/center_line.py index 20014ca..79a09e2 100755 --- a/jord/geometric_analysis/center_line.py +++ b/jord/geometric_analysis/center_line.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -from typing import Union, Tuple, List, Iterable +from typing import Iterable, List, Tuple, Union import numpy from numpy import array @@ -40,20 +40,20 @@ def find_centerline( int(min(ext_xy[1])), ) - vertices, ridges = _get_voronoi_vertices_and_ridges( + vertices, ridges = get_voronoi_vertices_and_ridges( input_geometry, step_size, minx=_min_x, miny=_min_y ) lines = [] for ridge in ridges: - if _ridge_is_finite(ridge): - starting_point = _create_point_with_restored_coordinates( + if ridge_is_finite(ridge): + starting_point = create_point_with_restored_coordinates( x=vertices[ridge[0]][0], y=vertices[ridge[0]][1], _min_x=_min_x, _min_y=_min_y, ) - ending_point = _create_point_with_restored_coordinates( + ending_point = create_point_with_restored_coordinates( x=vertices[ridge[1]][0], y=vertices[ridge[1]][1], _min_x=_min_x, @@ -61,7 +61,7 @@ def find_centerline( ) linestring = LineString((starting_point, ending_point)) - if _linestring_is_within_input_geometry(linestring, input_geometry): + if linestring_is_within_input_geometry(linestring, input_geometry): lines.append(linestring) if len(lines) < 2: @@ -70,94 +70,94 @@ def find_centerline( return MultiLineString(lines=unary_union(lines)) -def _get_voronoi_vertices_and_ridges( +def get_voronoi_vertices_and_ridges( _input_geometry: BaseGeometry, _step_size: float, minx: float, miny: float, ) -> Tuple[numpy.ndarray, List[List[int]]]: - borders = _get_densified_borders(_input_geometry, _step_size, minx, miny) + borders = densify_border(_input_geometry, _step_size, minx, miny) voronoi_diagram = Voronoi(borders) return voronoi_diagram.vertices, voronoi_diagram.ridge_vertices -def _ridge_is_finite(ridge: Iterable) -> bool: +def ridge_is_finite(ridge: Iterable) -> bool: return -1 not in ridge -def _create_point_with_restored_coordinates( +def create_point_with_restored_coordinates( x: float, y: float, _min_x: float, _min_y: float ) -> Tuple[float, float]: return (x + _min_x, y + _min_y) -def _linestring_is_within_input_geometry( +def linestring_is_within_input_geometry( linestring: LineString, input_geometry: BaseGeometry ) -> bool: return linestring.within(input_geometry) and len(linestring.coords[0]) > 1 -def _get_densified_borders( +def densify_border( _input_geometry: BaseGeometry, _step_size, minx: float, miny: float ) -> numpy.ndarray: polygons = iter_polygons(_input_geometry) points = [] for polygon in polygons: - points += _get_interpolated_boundary(polygon.exterior, _step_size, minx, miny) + points += interpolated_boundary(polygon.exterior, _step_size, minx, miny) if polygon_has_interior_rings(polygon): for interior in polygon.interiors: - points += _get_interpolated_boundary( + points += interpolated_boundary( interior, _step_size, minx=minx, miny=miny ) return array(points) -def _get_interpolated_boundary( +def interpolated_boundary( boundary: BaseGeometry, _step_size: float, minx: float, miny: float ) -> List[Tuple[float, float]]: line = LineString(boundary) return ( - [_get_coordinates_of_first_point(line, minx, miny)] - + _get_coordinates_of_interpolated_points( + [get_coordinates_of_first_point(line, minx, miny)] + + get_coordinates_of_interpolated_points( line, _step_size, min_x=minx, min_y=miny ) - + [_get_coordinates_of_last_point(line, minx=minx, miny=miny)] + + [get_coordinates_of_last_point(line, minx=minx, miny=miny)] ) -def _create_point_with_reduced_coordinates( +def create_point_with_reduced_coordinates( x: float, y: float, _min_x: float, _min_y: float ) -> Tuple[float, float]: return (x - _min_x, y - _min_y) -def _get_coordinates_of_first_point( +def get_coordinates_of_first_point( linestring: LineString, minx: float, miny: float ) -> Tuple[float, float]: - return _create_point_with_reduced_coordinates( + return create_point_with_reduced_coordinates( x=linestring.xy[0][0], y=linestring.xy[1][0], _min_x=minx, _min_y=miny ) -def _get_coordinates_of_last_point( +def get_coordinates_of_last_point( linestring: LineString, minx: float, miny: float ) -> Tuple[float, float]: - return _create_point_with_reduced_coordinates( + return create_point_with_reduced_coordinates( x=linestring.xy[0][-1], y=linestring.xy[1][-1], _min_x=minx, _min_y=miny ) -def _get_coordinates_of_interpolated_points( +def get_coordinates_of_interpolated_points( linestring: LineString, _step_size: Number, min_x: float, min_y: float ) -> List[Tuple[float, float]]: interpolation_distance = _step_size intermediate_points = [] while interpolation_distance < linestring.length: point = linestring.interpolate(interpolation_distance) - reduced_point = _create_point_with_reduced_coordinates( + reduced_point = create_point_with_reduced_coordinates( x=point.x, y=point.y, _min_x=min_x, _min_y=min_y ) intermediate_points.append(reduced_point) diff --git a/jord/geopandas_utilities/serialisation/well_known_binary.py b/jord/geopandas_utilities/serialisation/well_known_binary.py index 92a9e5f..1a089f6 100755 --- a/jord/geopandas_utilities/serialisation/well_known_binary.py +++ b/jord/geopandas_utilities/serialisation/well_known_binary.py @@ -2,21 +2,26 @@ from pathlib import Path +import pandas +import shapely.geometry.base from shapely import wkb __all__ = ["load_wkbs_from_csv", "csv_wkt_generator"] -def load_wkbs_from_csv(csv_file_path: Path, geometry_column: str = "Shape") -> wkb: +def load_wkbs_from_csv( + csv_file_path: Path, geometry_column: str = "Shape" +) -> pandas.DataFrame: """ Well-Known Text """ - import pandas return pandas.read_csv(str(csv_file_path))[geometry_column].apply(wkb.loads) -def csv_wkt_generator(csv_file_path: Path, geometry_column: str = "Shape") -> wkb: +def csv_wkt_generator( + csv_file_path: Path, geometry_column: str = "Shape" +) -> shapely.geometry.base.BaseGeometry: """ :param csv_file_path: @@ -28,4 +33,4 @@ def csv_wkt_generator(csv_file_path: Path, geometry_column: str = "Shape") -> wk for idx, g in pandas.read_csv( str(csv_file_path), usecols=[geometry_column] ).iterrows(): - yield wkb.loads(g) + yield wkb.loads(g) # g is pandas Series? diff --git a/jord/geopandas_utilities/serialisation/well_known_text.py b/jord/geopandas_utilities/serialisation/well_known_text.py index cef1627..fabe368 100755 --- a/jord/geopandas_utilities/serialisation/well_known_text.py +++ b/jord/geopandas_utilities/serialisation/well_known_text.py @@ -5,6 +5,7 @@ from typing import Sequence import pandas +import shapely.geometry.base from pandas import DataFrame from shapely import wkt @@ -46,7 +47,9 @@ def load_wkts_from_csv( return df -def csv_wkt_generator(csv_file_path: Path, geometry_column: str = "Shape") -> wkt: +def csv_wkt_generator( + csv_file_path: Path, geometry_column: str = "Shape" +) -> shapely.geometry.base.BaseGeometry: """ :param csv_file_path: @@ -56,7 +59,7 @@ def csv_wkt_generator(csv_file_path: Path, geometry_column: str = "Shape") -> wk for idx, g in pandas.read_csv( str(csv_file_path), usecols=[geometry_column] ).iterrows(): - yield wkt.loads(g) + yield wkt.loads(g) # g is a pandas Series? if __name__ == "__main__": diff --git a/jord/shapely_utilities/__init__.py b/jord/shapely_utilities/__init__.py index 7728da5..018071e 100755 --- a/jord/shapely_utilities/__init__.py +++ b/jord/shapely_utilities/__init__.py @@ -26,3 +26,4 @@ from .mirroring import * from .selection import * from .uniformity import * +from .subdivision import * diff --git a/jord/shapely_utilities/lines.py b/jord/shapely_utilities/lines.py index f6f3b54..6635d1e 100755 --- a/jord/shapely_utilities/lines.py +++ b/jord/shapely_utilities/lines.py @@ -35,6 +35,7 @@ "remove_redundant_nodes", "split_line", "internal_points", + "snap_endings_to_points", ] import collections @@ -53,7 +54,6 @@ box, ) from shapely.geometry.base import BaseGeometry -from shapely.ops import linemerge as shapely_linemerge # from sorcery import assigned_names from warg import Number, pairs @@ -158,6 +158,9 @@ def to_lines( else: raise NotImplementedError(f"{geoms, type(geoms)}") + if not isinstance(lines, List): + lines = list(lines) + return lines @@ -462,18 +465,18 @@ def snappy_endings( # find all vertices within a radius of max_distance as possible target = nearest_neighbor_within(snapping_points, endpoint, max_distance) - # do nothing if no target point to snap to is found - if not target: + if not target: # do nothing if no target point to snap to is found continue - # find the LineString to modify within snapped_lines and update it - for i, snapped_line in enumerate(snapped_lines): + for i, snapped_line in enumerate( + snapped_lines + ): # find the LineString to modify within snapped_lines and update it if endpoint.touches(snapped_line): snapped_lines[i] = bend_towards(snapped_line, where=endpoint, to=target) break - # also update the corresponding snapping_points for i, snapping_point in enumerate(snapping_points): + # also update the corresponding snapping_points if endpoint.equals(snapping_point): snapping_points[i] = target break @@ -484,10 +487,65 @@ def snappy_endings( return snapped_lines -def bend_towards(line: LineString, where: Point, to: Point) -> LineString: +def snap_endings_to_points( + lines: Union[Iterable[LineString], MultiLineString], + snapping_points: Sequence[Point], + max_distance: float, +) -> Sequence[Union[LineString, MultiLineString]]: + """ + Snap endpoints of lines together if they are at most max_length apart. + + + :param lines: A list of LineStrings or a MultiLineString + :param max_distance: maximum distance two endpoints may be joined together + :return: + :rtype: Sequence[Union[LineString, MultiLineString]] + """ + + # initialize snapped lines with list of original lines + # snapping points is a MultiPoint object of all vertices + snapped_lines = to_lines(lines) + + if isinstance(snapping_points, MultiPoint): + snapping_points = list(snapping_points.geoms) + + # isolated endpoints are going to snap to the closest vertex + isolated_endpoints = find_isolated_endpoints(snapped_lines) + + # only move isolated endpoints, one by one + for endpoint in isolated_endpoints: + # find all vertices within a radius of max_distance as possible + target = nearest_neighbor_within(snapping_points, endpoint, max_distance) + + if not target: # do nothing if no target point to snap to is found + continue + + for i, snapped_line in enumerate( + snapped_lines + ): # find the LineString to modify within snapped_lines and update it + if endpoint.touches(snapped_line): + snapped_lines[i] = bend_towards(snapped_line, where=endpoint, to=target) + break + + for i, snapping_point in enumerate(snapping_points): + # also update the corresponding snapping_points + if endpoint.equals(snapping_point): + snapping_points[i] = target + break + + # post-processing: remove any resulting lines of length 0 + snapped_lines = [s for s in snapped_lines if s.length > 0] + + return snapped_lines + + +def bend_towards( + line: LineString, where: Point, to: Point, tolerance: float = 1e-6 +) -> LineString: """ Move the point where along a line to the point at location to. + :param tolerance: :param line: :param where: a point ON the line (not necessarily a vertex) :param to: a point NOT on the line where the nearest vertex will be moved to @@ -500,7 +558,7 @@ def bend_towards(line: LineString, where: Point, to: Point) -> LineString: coords = line.coords[:] # easy case: where is (within numeric precision) a vertex of line for k, vertex in enumerate(coords): - if where.almost_equals(Point(vertex)): + if where.equals_exact(Point(vertex), tolerance=tolerance): # move coordinates of the vertex to destination coords[k] = to.coords[0] return LineString(coords) @@ -546,7 +604,7 @@ def prune_short_lines( def linemerge( line_s: Union[ - LineString, MultiLineString, Sequence[LineString], Sequence[MultiLineString] + LineString, MultiLineString, Iterable[LineString], Iterable[MultiLineString] ] ) -> Union[LineString, MultiLineString]: """ @@ -558,15 +616,24 @@ def linemerge( :type line_s: LineString|MultiLineString :rtype:LineString|MultiLineString """ - lines = [] assert isinstance(line_s, (LineString, MultiLineString, Sequence)) if isinstance(line_s, LineString): return line_s - if isinstance(line_s, Sequence): - return shapely_linemerge([linemerge(l) for l in line_s]) + lines = [] + + if isinstance(line_s, Iterable): + for l in line_s: + a = linemerge(l) + if isinstance(a, LineString): + lines.append(a) + elif isinstance(a, MultiLineString): + lines.extend(a.geoms) + else: + raise NotImplementedError(f"{type(a)} is not supported") + return shapely.ops.linemerge(lines) for line in line_s.geoms: if isinstance(line, MultiLineString): @@ -576,7 +643,7 @@ def linemerge( # line is a line, so simply append it lines.append(line) - return shapely_linemerge(lines) + return shapely.ops.linemerge(lines) def are_incident(v1, v2) -> bool: @@ -703,10 +770,13 @@ def extend_line( def extend_lines( - lines: Union[LineString, MultiLineString, Iterable[LineString]], distance: Number + lines: Union[LineString, MultiLineString, Iterable[LineString]], + distance: Number, + simplify: bool = False, ) -> List[LineString]: """ + :param simplify: :param lines: :param distance: :return: @@ -720,7 +790,7 @@ def extend_lines( out_lines = [] for l in lines: - out_lines.append(extend_line(l, distance)) + out_lines.append(extend_line(l, distance, simplify=simplify)) return out_lines diff --git a/jord/shapely_utilities/morphology.py b/jord/shapely_utilities/morphology.py index 55f0cc6..65625e0 100755 --- a/jord/shapely_utilities/morphology.py +++ b/jord/shapely_utilities/morphology.py @@ -15,6 +15,8 @@ "close", "clean_shape", "zero_buffer", + "pro_opening", + "pro_closing", "BecameEmptyException", "collapse_duplicate_vertices", ] diff --git a/jord/shapely_utilities/points.py b/jord/shapely_utilities/points.py index d10ea6b..6830bf5 100755 --- a/jord/shapely_utilities/points.py +++ b/jord/shapely_utilities/points.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 -from typing import Sequence, List, Optional, Union, Tuple, Iterable, Generator +from typing import Generator, Iterable, List, Optional, Sequence, Tuple, Union import numpy import shapely -from shapely.geometry import LineString, Point, MultiPoint +from shapely.geometry import LineString, MultiPoint, Point from shapely.geometry.base import BaseGeometry from warg import Number @@ -54,6 +54,9 @@ def nearest_neighbor_within( elif isinstance(interesting_points, Point): closest_point = interesting_points else: + if isinstance(interesting_points, MultiPoint): + interesting_points = list(interesting_points.geoms) + distances = [ point.distance(ip) for ip in interesting_points if point.distance(ip) > 0 ] diff --git a/jord/shapely_utilities/polygons.py b/jord/shapely_utilities/polygons.py index 6a32a56..3d2686f 100755 --- a/jord/shapely_utilities/polygons.py +++ b/jord/shapely_utilities/polygons.py @@ -452,5 +452,5 @@ def ahsudh(): p.plot() pyplot.show() - ahsudh() + # ahsudh() # aishdjauisd() diff --git a/jord/shapely_utilities/projection.py b/jord/shapely_utilities/projection.py index 2b18a53..cf95938 100755 --- a/jord/shapely_utilities/projection.py +++ b/jord/shapely_utilities/projection.py @@ -194,10 +194,15 @@ def get_min_max_projected_line( if not isinstance(geom, shapely.LineString): geom = geom.boundary + if isinstance(other, shapely.MultiLineString): + geom = geom.boundary + min_v = max_v = 0.5 + other_coords = other.boundary.coords + # Find limits - for point_coords in other.boundary.coords: + for point_coords in other_coords: v = geom.project(Point(point_coords), normalized=True) if v < min_v: min_v = v diff --git a/jord/shapely_utilities/rings.py b/jord/shapely_utilities/rings.py index 7f7c30b..d860c52 100755 --- a/jord/shapely_utilities/rings.py +++ b/jord/shapely_utilities/rings.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 +import shapely from shapely import LinearRing +from shapely.geometry import LineString, MultiLineString, Point, Polygon +from shapely.ops import linemerge -__all__ = [ - "ensure_ccw_ring", - "ensure_cw_ring", -] +__all__ = ["ensure_ccw_ring", "ensure_cw_ring", "split_ring"] + +from typing import Union def ensure_ccw_ring(ring: LinearRing) -> LinearRing: @@ -18,3 +20,89 @@ def ensure_cw_ring(ring: LinearRing) -> LinearRing: if ring.is_ccw: return LinearRing(list(ring.coords)[::-1]) return ring + + +def split_ring( + ring: LinearRing, + split: Union[LinearRing, LineString, MultiLineString, shapely.GeometryCollection], +): + """Split a linear ring geometry, returns a [Multi]LineString + + See my PostGIS function on scigen named ST_SplitRing + """ + valid_types = ("MultiLineString", "LineString", "GeometryCollection") + if not hasattr(ring, "geom_type"): + raise ValueError("expected ring as a geometry") + elif not hasattr(split, "geom_type"): + raise ValueError("expected split as a geometry") + if ring.geom_type == "LinearRing": + ring = LineString(ring) + if ring.geom_type != "LineString": + raise ValueError( + "ring is not a LinearRing or LineString, found " + str(ring.geom_type) + ) + elif not ring.is_closed: + raise ValueError("ring is not closed") + elif split.is_empty: + return ring + elif not split.intersects(ring): + # split does not intersect ring + return ring + if split.geom_type == "LinearRing": + split = LineString(split) + if split.geom_type not in valid_types: + raise ValueError( + "split is not a LineString-like or GeometryCollection geometry, " + "found " + str(split.geom_type) + ) + + intersections = ring.intersection(split) + if intersections.is_empty: + # no intersections, returning same ring + return ring + elif intersections.geom_type == "Point": + # Simple case, where there is only one line intersecting the ring + result = Polygon(ring).difference(split).exterior + # If it is a coordinate of the ring, then the ring needs to be rotated + coords = result.coords[:-1] + found_i = 0 + for i, c in enumerate(coords): + if Point(c).almost_equals(intersections): + found_i = i + break + if found_i > 0: + result = Polygon(coords[i:] + coords[:i]).exterior + if result.interpolate(0).distance(intersections) > 0: + raise Exception( + "result start point %s to intersection %s is %s" + % (result.interpolate(0), intersections, result.distance(intersections)) + ) + elif result.geom_type != "LinearRing": + raise Exception("result is not a LinearRing, found " + result.geom_type) + elif not result.is_closed: + raise Exception("result is not closed") + return LineString(result) + + difference = ring.difference(split) + if difference.geom_type != "MultiLineString": + raise ValueError( + "expected MultiLineString difference, found " + difference.geom_type + ) + + start_point = ring.interpolate(0) + if start_point.distance(intersections) == 0: + # special case: start point is the same as an intersection + return difference + + # Otherwise the line where the close meets needs to be fused + fuse = [] + parts = list(difference.geoms) + for ipart, part in enumerate(parts): + if part.intersects(start_point): + fuse.append(ipart) + if len(fuse) != 2: + raise ValueError("expected 2 geometries, found " + str(len(fuse))) + # glue the last to the first + popped_part = parts.pop(fuse[1]) + parts[fuse[0]] = linemerge([parts[fuse[0]], popped_part]) + return MultiLineString(parts) diff --git a/jord/shapely_utilities/sliverin.py b/jord/shapely_utilities/sliverin.py index 2df37af..1730891 100644 --- a/jord/shapely_utilities/sliverin.py +++ b/jord/shapely_utilities/sliverin.py @@ -1,14 +1,31 @@ from collections import defaultdict -from typing import Collection, Dict, List, Sequence +from copy import copy +from typing import Collection, Dict, List, Sequence, Union import shapely +from shapely.constructive import simplify +from shapely.geometry.base import GeometrySequence +from tqdm import tqdm -from jord.shapely_utilities import dilate, iter_polygons +from jord.geometric_analysis import construct_centerline +from jord.shapely_utilities import dilate, extend_line, is_multi, iter_polygons from jord.shapely_utilities.desliver_wkt import a_wkt -from jord.shapely_utilities.morphology import pro_opening +from jord.shapely_utilities.lines import ( + ExtensionDirectionEnum, + find_isolated_endpoints, + linemerge, + snap_endings_to_points, +) +from jord.shapely_utilities.morphology import closing, erode, pro_opening +from jord.shapely_utilities.subdivision import subdivide __all__ = ["desliver"] +from jord.shapely_utilities.projection import ( + get_min_max_projected_line, + project_point_to_object, +) + def desliver( polygons: Collection[shapely.Polygon], buffer_size: float = 0.2 @@ -25,7 +42,158 @@ def desliver( return buffered_exterior -def desliver2( +def cut_polygon( + polygon: shapely.Polygon, line_split_collection: List[shapely.LineString] +) -> GeometrySequence: + line_split_collection.append( + polygon.boundary + ) # collection of individual linestrings for splitting in a list and add the polygon lines to it. + merged_lines = linemerge(line_split_collection) + border_lines = shapely.ops.unary_union(merged_lines) + return shapely.ops.polygonize(border_lines) + + +def multi_line_extend( + multi_line_string: Union[shapely.LineString, shapely.MultiLineString], + distance: float, +) -> shapely.MultiLineString: + isolated_endpoints = find_isolated_endpoints(multi_line_string) + + lines = [] + + if isinstance(multi_line_string, shapely.LineString): + ls = [multi_line_string] + else: + ls = multi_line_string.geoms + + for line in ls: + start_point, end_point = shapely.Point(line.coords[0]), shapely.Point( + line.coords[-1] + ) + + endpoint_in_isolated_points = end_point in isolated_endpoints + + direction = None + if start_point in isolated_endpoints: + if endpoint_in_isolated_points: + direction = ExtensionDirectionEnum.both + else: + direction = ExtensionDirectionEnum.start + elif endpoint_in_isolated_points: + direction = ExtensionDirectionEnum.end + + if direction is not None: + line = extend_line(line, offset=distance, simplify=False, side=direction) + + lines.append(line) + + return shapely.MultiLineString(lines) + + +def desliver_center_divide( + polygons: Collection[shapely.Polygon], buffer_size: float = 0.2 +) -> List[shapely.geometry.Polygon]: + buffered_exterior = [] + + if isinstance(polygons, Sequence): + polygons = list(polygons) + + for polygon in polygons: + polygon: shapely.Polygon + buffered_exterior.append(dilate(polygon, distance=buffer_size) - polygon) + + augmented_polygons = [] + + intersections = [] + for ith in range(len(buffered_exterior)): + a = buffered_exterior.copy() + b = a.pop(ith) + intersections.append(shapely.unary_union(a) & b) + + for ith, intersection in tqdm(enumerate(intersections)): + some_distance = intersection.minimum_clearance / 4.0 + + center_line = construct_centerline( + intersection, interpolation_distance=some_distance + ) + + center_line = simplify( + center_line, preserve_topology=False, tolerance=some_distance + ) + + center_line = simplify( + center_line, preserve_topology=True, tolerance=some_distance + ) + + center_line = multi_line_extend(center_line, distance=some_distance) + + if isinstance(intersection, shapely.Polygon): + snapping_points = [ + shapely.Point(c) for c in subdivide(intersection).exterior.coords + ] + else: + snapping_points = [ + shapely.Point(c) + for inter in intersection.geoms + for c in subdivide(inter).exterior.coords + ] + + snapped_center_line = snap_endings_to_points( + center_line, snapping_points=snapping_points, max_distance=some_distance + ) + + poly = polygons[ith] + + if True: + for line in snapped_center_line.copy(): + projected_line = get_min_max_projected_line(line, poly) + + start, end = shapely.Point(projected_line.coords[0]), shapely.Point( + projected_line.coords[-1] + ) + + start_line, end_line = ( + extend_line( + shapely.LineString( + (start, project_point_to_object(start, poly)) + ), + offset=some_distance, + ), + extend_line( + shapely.LineString((end, project_point_to_object(end, poly))), + offset=some_distance, + ), + ) + + snapped_center_line.extend((start_line, end_line)) + + res = cut_polygon(intersection, snapped_center_line) + + augmented = copy(poly) + for r in res: + un = r | poly + re = erode(dilate(un, distance=1e-10), distance=1e-9) + if is_multi(re): + continue + + f = closing(un) + + if True: + augmented |= f + else: + k = r & poly + if k: + if isinstance(k, shapely.LineString): + if k.length >= some_distance: + augmented |= r + + augmented = pro_opening(augmented, distance=some_distance) + augmented_polygons.append(augmented) + + return augmented_polygons + + +def desliver_least_intersectors_first( polygons: Collection[shapely.Polygon], buffer_size: float = 0.2 ) -> Dict[int, shapely.geometry.Polygon]: buffered_exterior = [] @@ -90,8 +258,8 @@ def desliver2( def sauihd2(): polygons = list(shapely.from_wkt(a_wkt).geoms) - once = desliver2(polygons) - out = desliver2(list(once.values())) + once = desliver_least_intersectors_first(polygons) + out = desliver_least_intersectors_first(list(once.values())) c = shapely.MultiPolygon(iter_polygons(out.values())) ... @@ -104,4 +272,13 @@ def sauihd(): c = shapely.MultiPolygon(iter_polygons(out.values())) ... - sauihd() + def sauihd3(): + polygons = list(shapely.from_wkt(a_wkt).geoms) + once = desliver_center_divide(polygons) + out = desliver_center_divide(list(once.values())) + + c = shapely.MultiPolygon(iter_polygons(out.values())) + print(c.wkt) + ... + + sauihd3() diff --git a/jord/shapely_utilities/subdivision.py b/jord/shapely_utilities/subdivision.py new file mode 100644 index 0000000..9d067f5 --- /dev/null +++ b/jord/shapely_utilities/subdivision.py @@ -0,0 +1,76 @@ +from typing import Union + +import shapely +from shapely.geometry.linestring import LineString +from shapely.geometry.polygon import LinearRing, Polygon +from warg import pairs + +from jord.shapely_utilities.lines import linemerge + +__all__ = ["subdivide", "subdivide_polygon", "subdivide_line", "subdivide_ring"] + + +def subdivide( + geom: Union[LineString, LinearRing, Polygon] +) -> Union[LineString, LinearRing, Polygon]: + if isinstance(geom, LineString): + return subdivide_line(geom) + elif isinstance(geom, LinearRing): + return subdivide_ring(geom) + elif isinstance(geom, Polygon): + return subdivide_polygon(geom) + + raise NotImplementedError(f"Subdivision for {type(geom)} not implemented") + + +def subdivide_line(line: LineString) -> LineString: + half_point = line.interpolate(0.5, normalized=True) + return shapely.LineString((line.coords[0], *(half_point.coords), line.coords[-1])) + + +def subdivide_ring(ring: LinearRing) -> LinearRing: + ring_segments = [] + for segment in pairs(list(ring.coords)): + ring_segments.append(subdivide_line(LineString(segment))) + + return shapely.LinearRing(linemerge(ring_segments)) + + +def subdivide_polygon(polygon: shapely.Polygon) -> shapely.Polygon: + exterior = subdivide_ring(polygon.exterior) + + interiors = [] + for interior in polygon.interiors: + interiors.append(subdivide_ring(interior)) + + return Polygon(exterior, holes=interiors) + + +if __name__ == "__main__": + + def uihasud(): + from jord.shapely_utilities import dilate + + a = dilate( + shapely.Point((0, 0)), distance=1, cap_style=shapely.BufferCapStyle.square + ) + + print(a.exterior) + print(subdivide_ring(a.exterior)) + + def uhasd(): + from jord.shapely_utilities import dilate + + a = dilate( + shapely.Point((0, 0)), distance=1, cap_style=shapely.BufferCapStyle.square + ) + b = dilate( + shapely.Point((0, 0)), distance=0.5, cap_style=shapely.BufferCapStyle.square + ) + c = a - b + + print(c) + print(subdivide_polygon(c)) + + uihasud() + uhasd() diff --git a/tests/shapely/test_polygon_utilities.py b/tests/shapely/test_polygon_utilities.py index 7fbc162..5e3662e 100644 --- a/tests/shapely/test_polygon_utilities.py +++ b/tests/shapely/test_polygon_utilities.py @@ -22,7 +22,7 @@ def third_point() -> shapely.Point: @pytest.fixture(scope="module") -def a_line(a_point: shapely.Point, another_point: shapely.Point) -> shapely.Point: +def a_line(a_point: shapely.Point, another_point: shapely.Point) -> shapely.LineString: return shapely.LineString([a_point, another_point]) diff --git a/tests/shapely/test_sliverin.py b/tests/shapely/test_sliverin.py index 0c635cc..1f6a547 100644 --- a/tests/shapely/test_sliverin.py +++ b/tests/shapely/test_sliverin.py @@ -1,7 +1,10 @@ import shapely from jord.shapely_utilities import dilate -from jord.shapely_utilities.sliverin import desliver2 +from jord.shapely_utilities.sliverin import ( + desliver_center_divide, + desliver_least_intersectors_first, +) def test_union_buf(): @@ -25,7 +28,7 @@ def test_union_buf(): print(res.wkt) -def test_desliver(): +def test_desliver_least_intersectors(): buffered_exterior = [ dilate( shapely.Point((0, 0)), cap_style=shapely.BufferCapStyle.square, distance=0.9 @@ -38,8 +41,42 @@ def test_desliver(): ), ] - print(shapely.MultiPolygon(buffered_exterior).wkt) - - res = desliver2(buffered_exterior, buffer_size=0.2) + res = desliver_least_intersectors_first(buffered_exterior, buffer_size=0.2) print(shapely.MultiPolygon(list(res.values())).wkt) + + +def test_desliver_intersection_center_distribute(): + buffered_exterior = [ + dilate( + shapely.Point((0, 0)), cap_style=shapely.BufferCapStyle.square, distance=0.9 + ), + dilate( + shapely.Point((2, 0)), cap_style=shapely.BufferCapStyle.square, distance=0.9 + ), + dilate( + shapely.Point((4, 0)), cap_style=shapely.BufferCapStyle.square, distance=0.9 + ), + ] + + res = desliver_center_divide(buffered_exterior, buffer_size=0.2) + + print(shapely.GeometryCollection(list(res)).wkt) + + +def test_desliver_intersection_center_distribute_circle(): + buffered_exterior = [ + dilate( + shapely.Point((0, 0)), cap_style=shapely.BufferCapStyle.round, distance=0.9 + ), + dilate( + shapely.Point((2, 0)), cap_style=shapely.BufferCapStyle.round, distance=0.9 + ), + dilate( + shapely.Point((4, 0)), cap_style=shapely.BufferCapStyle.round, distance=0.9 + ), + ] + + res = desliver_center_divide(buffered_exterior, buffer_size=0.3) + + print(shapely.GeometryCollection(list(res)).wkt)