CODEWITHSUNDEEP
pythonshapely
Robbes
Robbes

Reputation: 145

Shapely - insert a single point into a Polygon from a projection

I am trying to insert a point from a projected line into a polygon (I need to do this a few times to 'cut' a section of the polygon exterior, but this the key issue).

There is close similarity with this question, but I find that answer does or doesn't work depending on the segment orientation. It works for simple vertical or horizontal projections, but in most other segment orientations results in a fail due to the projected point falling outside the polygon.exterior by ~ 1e-16 (numerical precision).
This question checks for a point being inside the polygon, and also fails when the point is on the polygon exterior.

I have tried extending the projected line (pline below) by a tiny bit then intersection, but that gives the exact same Point and numerical precision problem.
Avoiding the precision by adding a tiny .buffer(1e-15) does intercept, but that results in adding a (tiny) rounded line segment from the buffer which introduces errors elsewhere.

Example:

import matplotlib.pyplot as plt
import shapely as sh
from shapely import ops
from shapely.geometry import *

def insert_projected_point(point, poly):
    s = sh.normalize(Polygon(poly))
    pline = sh.shortest_line(point,s)
    p2 = sh.intersection(s, pline)
    # Per https://gis.stackexchange.com/questions/188594/how-can-i-add-points-to-a-linestring-in-shapely
    ring = LineString(list(s.exterior.coords))
    union = ring.union(pline)
    s2 = ops.polygonize(union)[0]
    return p2, s2

def do_results(p, s, p1, s1):
    print(f"s: {s}")
    print(f"p out: {p1}")
    print(f"s out: {s1}")


p = Point((3, 0.5))

# This works
poly1 = (0,0), (2,0), (2,2), (0,2), (0,0)
p1, s1 = insert_projected_point(p, poly1)
do_results(p, poly1, p1, s1)

# This doesn't work
poly2 = (0,0), (1,0), (2,2), (0,2), (0,0)
p2, s2 = insert_projected_point(p, poly2)
do_results(p, poly2, p2, s2)

Result:

s: ((0, 0), (2, 0), (2, 2), (0, 2), (0, 0))
p out: POINT (2 0.5)
sout: POLYGON ((0 0, 0 2, 2 2, 2 0.5, 2 0, 0 0))

s: ((0, 0), (1, 0),(2, 2), (0, 2), (0, 0))
p out: LINESTRING Z EMPTY
s out: POLYGON((0 0, 0 2, 2 2, 1 0, 0 0))

Upvotes: 2

Views: 452

Answers (1)

Robbes
Robbes

Reputation: 145

Solution to extend the intersect line by a tiny fraction works now (not sure where I went wrong earlier)

def insert_projected_point(point, poly, scaler=1):
    s = sh.normalize(Polygon(poly))
    pline = sh.shortest_line(point,s)
    
    # add this line to make it work
    pline = sh.affinity.scale(pline, xfact=scaler, yfact=scaler, origin='center')
    
    p2 = sh.intersection(s, pline)
    # Per https://gis.stackexchange.com/questions/188594/how-can-i-add-points-to-a-linestring-in-shapely
    ring = LineString(list(s.exterior.coords))
    union = ring.union(pline)
    s2 = ops.polygonize(union)[0]
    print(f"p out: {p2}")
    print(f"s out: {s2}")

p = Point((3, 0.5))
poly2 = (0,0), (1,0), (2,2), (0,2), (0,0)

# This doesn't work
insert_projected_point(p, poly2)
# p out: LINESTRING Z EMPTY
# s out: POLYGON ((0 0, 0 2, 2 2, 1 0, 0 0))

# This works
insert_projected_point(p, poly2, scaler=1 + 1e-12)
# p out: LINESTRING (1.6 1.2, 1.5999999999992998 1.2000000000003501)
# s out: POLYGON ((0 0, 0 2, 2 2, 1.6 1.2, 1 0, 0 0))

Upvotes: 0

Related Questions