diff --git a/src/c3nav/mapdata/models/geometry/level.py b/src/c3nav/mapdata/models/geometry/level.py index c9c62052..0d0784cb 100644 --- a/src/c3nav/mapdata/models/geometry/level.py +++ b/src/c3nav/mapdata/models/geometry/level.py @@ -19,7 +19,8 @@ from c3nav.mapdata.models import Level from c3nav.mapdata.models.access import AccessRestrictionMixin from c3nav.mapdata.models.geometry.base import GeometryMixin from c3nav.mapdata.models.locations import SpecificLocation -from c3nav.mapdata.utils.geometry import assert_multilinestring, assert_multipolygon, clean_geometry +from c3nav.mapdata.utils.geometry import (assert_multilinestring, assert_multipolygon, clean_geometry, + cut_polygon_with_line) class LevelGeometryMixin(GeometryMixin): @@ -193,7 +194,7 @@ class AltitudeArea(LevelGeometryMixin, models.Model): if not area.geometry_prep.intersects(stair): continue - divided = assert_multipolygon(area.geometry.difference(stair.buffer(0.0001))) + divided = assert_multipolygon(cut_polygon_with_line(area.geometry, stair)) if len(divided) > 2: raise ValueError area.geometry = divided[0] diff --git a/src/c3nav/mapdata/utils/geometry.py b/src/c3nav/mapdata/utils/geometry.py index c8ed0902..a494cd81 100644 --- a/src/c3nav/mapdata/utils/geometry.py +++ b/src/c3nav/mapdata/utils/geometry.py @@ -1,11 +1,14 @@ +from collections import deque, namedtuple from itertools import chain +from typing import List import matplotlib.pyplot as plt from django.core import checks from matplotlib.patches import PathPatch from matplotlib.path import Path from shapely import speedups -from shapely.geometry import LinearRing, LineString, Polygon +from shapely.geometry import LinearRing, LineString, Point, Polygon +from shapely.ops import unary_union if speedups.available: speedups.enable() @@ -107,3 +110,86 @@ def get_rings(geometry): return (geometry, ) return () + + +cutpoint = namedtuple('cutpoint', ('point', 'polygon', 'ring')) + + +def cut_line_with_point(line, point): + distance = line.project(point) + pointlist = [(point.x, point.y)] + for i, p in enumerate(line.coords): + subdistance = line.project(Point(p)) + if subdistance >= distance: + return (LineString(line.coords[:i] + pointlist), + LineString(pointlist + line.coords[i+(1 if subdistance == distance else 0):])) + + +def cut_polygon_with_line(polygon, line): + polygons: List[List[LinearRing]] = [[polygon.exterior, *polygon.interiors] + for polygon in assert_multipolygon(polygon)] + + points = deque() + for i, polygon in polygons: + for j, ring in polygon: + intersection = ring.intersection(line) + for item in getattr(intersection, 'geoms', (intersection, )): + if isinstance(item, Point): + points.append(cutpoint(item, i, j)) + elif isinstance(item, LineString): + points.append(cutpoint(item.coords[0], i, j)) + points.append(cutpoint(item.coords[-1], i, j)) + else: + raise ValueError + + points = deque(sorted(points, key=lambda p: line.project(p.point))) + + last = points.popleft() + while points: + current = points.popleft() + if current.polygon == last.polygon: + polygon = polygons[current.polygon] + if current.ring != last.ring: + ring1 = cut_line_with_point(polygon[last.ring], last.point) + ring2 = cut_line_with_point(polygon[current.ring], current.point) + new_ring = LinearRing(ring1[:-1] + ring1[0:] + ring2[:-1] + ring2[0:] + ring1[:1]) + if current.ring == 0 or last.ring == 0: + polygon[0] = new_ring + interior = current.ring if last.ring == 0 else last.ring + polygon[interior] = None + mapping = {interior: 0} + else: + new_i = len(polygon) + mapping = {last.ring: new_i, current.ring: new_i} + polygon.append(new_ring) + polygon[last.ring] = None + polygon[current.ring] = None + + points = deque((cutpoint(item.point, item.polygon, mapping[item.ring]) + if (item.polygon == current.polygon and item.ring in mapping) else item) + for item in points) + elif current.ring == 0: + new_i = len(polygons) + exterior = cut_line_with_point(polygon[0], current.point) + exterior = cut_line_with_point(LinearRing(exterior[:-1] + exterior[0:]), last.point) + exterior1 = LinearRing(exterior[0][0:] + exterior[0][:1]) + exterior2 = LinearRing(exterior[1][0:] + exterior[1][:1]) + geom = Polygon(exterior1) + polygon[0] = exterior1 + new_polygon = [exterior2] + polygons.append(new_polygon) + mapping = {} + for i, interior in enumerate(polygon[1:]): + if not geom.contains(interior): + mapping[i] = len(new_polygon) + new_polygon.append(interior) + + points = deque((cutpoint(item.point, new_i, mapping[item.ring]) + if (item.polygon == current.polygon and item.ring in mapping) else item) + for item in points) + points = deque((cutpoint(item.point, new_i, 0) + if (item.polygon == current.polygon and item.ring == 0 and + not exterior1.contains(item.point)) else item) + for item in points) + + return unary_union(tuple(Polygon(polygon[0], polygon[1:]) for polygon in polygons))