Problem with point inside intersections polygon

[JanPolcher]

Hello,

When computing the overlap area of two polygons I some times find values which are far too large.
I could show that when increasing the number of points defining one of the polygons, then suddenly the overlap area becomes correct.
It turns out, that it is when the location of the center of the intersecting polygon is outside both original polygons (and thus also of the intersection) things go wrong.
Below is a simple piece of code which demonstrates the issues. It also displays the two polygons to show that as we change the number of points defining.

The output of this code :

> Version of SphericalGeometry :  1.2.19.dev0+gfbdc54a.d20200202
> Areas : large polygon [m^2]: 409982511 small polygon [m^2]: 669479
> ==================
> Small polygon of  5  points :
> Location of point inside intersecting polygon : 170.55 -38.74
> Area of small polygon [m^2]:  669479.69
> Do the two polygons intersect ? :  True
> Overlap with large polygon :  761642762.3568
> ==================
> Small polygon of  13  points :
> Location of point inside intersecting polygon : 170.55 -38.74
> Area of small polygon [m^2]:  669479.62
> Do the two polygons intersect ? :  True
> Overlap with large polygon :  761642844.3586
> ==================
> Small polygon of  21  points :
> Location of point inside intersecting polygon : 170.55 -38.74
> Area of small polygon [m^2]:  669479.77
> Do the two polygons intersect ? :  True
> Overlap with large polygon :  761642680.3551
> ==================
> Small polygon of  25  points :
> Location of point inside intersecting polygon : 350.55 38.74
> Area of small polygon [m^2]:  669480.05
> Do the two polygons intersect ? :  True
> Overlap with large polygon :  0.0061
> ==================
> Small polygon of  41  points :
> Location of point inside intersecting polygon : 350.55 38.74
> Area of small polygon [m^2]:  669480.05
> Do the two polygons intersect ? :  True
> Overlap with large polygon :  0.0061
> ==================

import sys
# 
import numpy as np
import spherical_geometry
from spherical_geometry import polygon
from spherical_geometry import vector
print("> Version of SphericalGeometry : ",spherical_geometry.__version__)
#
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from mpl_toolkits.basemap import Basemap
#
# A function to put a square polygon around a central point.
#
def boxit(cent, dlon, dlat, polyres) :
    boxll=[]
    loninc = dlon/polyres
    latinc = dlat/polyres
    # Side 1
    for pts in range(polyres) :
        boxll.append([cent[0]-dlon/2.0+loninc*pts, cent[1]+dlat/2.0])
    # Side 2
    for pts in range(polyres) :
        boxll.append([cent[0]+dlon/2.0, cent[1]+dlat/2.0-latinc*pts])
    # Side 3
    for pts in range(polyres) :
        boxll.append([cent[0]+dlon/2.0-loninc*pts, cent[1]-dlat/2.0])
    # Side 4
    for pts in range(polyres) :
        boxll.append([cent[0]-dlon/2.0, cent[1]-dlat/2.0+latinc*pts])
    # Close
    boxll.append(boxll[0])
    return boxll
#
#
#
EarthRadius=6370000.
polyll=[[-9.419769166656238, 38.653834491762225], \
        [-9.307308522254857, 38.677907248545154], \
        [-9.194765066594584, 38.701862566276745], \
        [-9.225426871471143, 38.78972563072393], \
        [-9.256171496793854, 38.87759011378188], \
        [-9.36899613899891, 38.85356451794185], \
        [-9.481737314455984, 38.82942115964239], \
        [-9.450711513172621, 38.74162715537378], \
        [-9.419769166656238, 38.653834491762225]]
centll=[[-9.33811076140455, 38.76573519339605]]
#
# Centre of small polygons
#
centsmall=[-9.45464, 38.73833]
#
sphpoly=polygon.SphericalPolygon.from_lonlat([p[0] for p in polyll],[p[1] for p in polyll], center=[centll[0][0],centll[0][1]])
#
DisplayLargeBox=False
DisplayIntersec=False
#
if DisplayLargeBox :
    # Display full large box
    dist=19000.0
    loncent=centll[0][0]
    latcent=centll[0][1]
else :
    # Only display mall Polygon
    dist=1000.0
    loncent=centsmall[0]
    latcent=centsmall[1]
#
m = Basemap(projection='ortho',resolution='i', lon_0=loncent, lat_0=latcent,  llcrnrx=-dist, llcrnry=-dist, urcrnrx=dist, urcrnry=dist)
m.drawcoastlines()
sphpoly.draw(m, color='red')
#
#
color=['b', 'g', 'c', 'm', 'y', 'k']
dlon=0.004166499711573124*2.0
dlat=0.004166498314589262*2.0
lstnbpts=[1, 3, 5, 6, 10]
#
#
small = boxit(centsmall, dlon, dlat, 1)
smallpolyA=polygon.SphericalPolygon.from_lonlat([p[0] for p in small], [p[1] for p in small], center=centsmall)
print("> Areas : large polygon [m^2]:", int(sphpoly.area()*EarthRadius**2), "small polygon [m^2]:",int(smallpolyA.area()*EarthRadius**2))
print("> ==================")
#
#
for ic,nbpts in enumerate(lstnbpts) :
    small = boxit(centsmall, dlon, dlat, nbpts)
    smallpolyA=polygon.SphericalPolygon.from_lonlat([p[0] for p in small], [p[1] for p in small], center=centsmall)
    inter=smallpolyA.intersection(sphpoly)
    #
    print("> Small polygon of ",len(small)," points :")
    ininter=inter.polygons[0].inside
    llininter=vector.vector_to_lonlat(ininter[0],ininter[1],ininter[2])
    print("> Location of point inside intersecting polygon :", "%.2f"%llininter[0], "%.2f"%llininter[1])
    #
    if DisplayIntersec :
        print("> Area of intersecting polygon [m^2]: ", "%.2f"%(inter.area()*EarthRadius**2))
        inter.draw(m, color=color[ic])
    else :
        print("> Area of small polygon [m^2]: ", "%.2f"%(smallpolyA.area()*EarthRadius**2))
        smallpolyA.draw(m, color=color[ic])
    #
    print("> Do the two polygons intersect ? : ", smallpolyA.intersects_poly(sphpoly))
    overlap = smallpolyA.overlap(sphpoly)
    print("> Overlap with large polygon : ", "%.4f"%overlap)
    print("> ==================")
#
#
#
plt.show()

[JanPolcher]

This problem seems to come from the precision of the calculation of the inside point in function. polygon._find_new_inside returns a point which then in graph.intersection is detected as outside of the polygon and the function polygon.invert_polygon() is called to turn around the polygon.

As solution to this problem of precision in the calculation of the inside point can be to take the mean of all candidates. It can be achieved with this piece of code in polygon._find_new_inside :

            avgmid=np.mean(np.array(midpoint),axis=0)
            vector.normalize_vector(avgmid, output=avgmid)
            candidate = max(zip(orient, midpoint), key=lambda x: x[0])
            # candidate[1] is the original ouput but avgmid gives more stable results
            inside = avgmid

The example given above now works flawlessly.

[pllim]

3 years too late but are you interested in opening a PR? 😸

[mcara]

I do not believe there is an issue with "precision" of computations. I may be wrong, but it seems to me that the polygon is self-intersecting and, to the best of my understanding (@pllim please correct me), self-intersecting polygons are not supported.

Even if self-intersecting polygons are supported, proposed fix works simply because of a lucky set of points. It may work for this particular set of points and it will fail for another set. I would not recommend making a PR to implement proposed "fix".

[pllim]

I don't remember what is and is not supported here. Though generally speaking if something isn't supported, we should prevent it from being created in the first place, or throw a meaningful error later during computation if the former isn't possible. Or if neither of those are possible, then document it clearly in the docs. What do you think?

[JanPolcher]

Nice to see that this issue is being worked on !
I do not think that the polygon in my test is self-intersecting.
It is certainly small and that is why my assumption was a rounding error.

[mcara]

I don't remember what is and is not supported here.

I never knew (== I am not the original developer) what is and is not supported here. However, I do not believe self-intersecting polygons are supported. You will need to split them into independent polygons and form a union.

Here is how @JanPolcher's polygon looks when projected to a tangent plane perpendicular to a middle point (vector) of points 1-6:

It is also possible I am doing something wrong so please correct me if you get a different image of your polygon.

[mcara]

Though generally speaking if something isn't supported, we should prevent it from being created in the first place

How?

or throw a meaningful error later during computation if the former isn't possible.

Feel free to make a PR for this.

Or if neither of those are possible, then document it clearly in the docs.

That is doable (==something even I can do).

[pllim]

Theoretically, we could check if any of the lines intersect each other but my math is real rusty now.

ref: https://math.stackexchange.com/questions/149622/finding-out-whether-two-line-segments-intersect-each-other

[mcara]

On a plane - yes. On a sphere, math is more complicated and time consuming. Performance is an important consideration.

[pllim]

In that case, doc would be a reasonable solution for now.

[JanPolcher]

OK, and I need to check in my code which produces the intersecting polygone. It should be the result of the intersects_poly function. But that is already 3 years ago and I do not know if we did not solve it differently.

[JanPolcher]

On a plane - yes. On a sphere, math is more complicated and time consuming. Performance is an important consideration.

Yes, performance in Spherical Geometry is super important !

[mcara]

It looks like self-intersecting polygons should be supported. See #131