Calculation of intersection polygon is stuck, does not even bring up an error.

[JanPolcher]

When two polygons have nearly overlapping segments it is impossible to compute their intersecting polygon.

The example below shows that we can test if the two polygons intersect but when we ask to compute the intersecting polygons the code is stuck. It does not even throw up an error.

import numpy as np
import matplotlib.pyplot as plt
import os,sys

sys.path += ["./spherical_geometry"]
from spherical_geometry import polygon

# Large Cell
lc_lons=[2.1145180961198093, 2.2167924672941695, 2.2166355474157626,  2.1142279214148516,  2.1145180961198093]
lc_lats=[45.59579076517927, 45.595946981135235,  45.66751087193568, 45.66735444850059, 45.59579076517927]
lc_cent=[2.1655435080610346, 45.6316621889731]

#Small Cell
sc_lons=[2.200007549639062, 2.2166742168709477, 2.2166742168709477, 2.200007549639062, 2.200007549639062]
sc_lats=[45.66666615789904, 45.66666615789904, 45.649999491515025, 45.649999491515025, 45.66666615789904]
sc_cent=[2.208341, 45.658333]

lc=polygon.SphericalPolygon.from_lonlat(lc_lons, lc_lats, center=lc_cent)
sc=polygon.SphericalPolygon.from_lonlat(sc_lons, sc_lats, center=sc_cent)

plt.plot(sc_lons, sc_lats, c="r")
plt.plot(sc_cent[0], sc_cent[1], marker="x", c="r")
plt.plot(lc_lons, lc_lats, c="g")
plt.plot(lc_cent[0], lc_cent[1], marker="x", c="g")
plt.show()

print("Is there an intersection ? ", lc.intersects_poly(sc))
inter=lc.intersection(sc)
print("Area of intersection :", inter.area())

It is sufficient to shift one of polygons slightly to get it working again ... but then obviously the overlapping area is wrong !
I hope this description of the bug helps solve it.

[JanPolcher]

I just noted that if I move one of the polygons by 3e-07 (in longitude/latitude system) the intersection polygon is correctly computed. Below this shift the above code is stuck.
So it appears to me that if two segments are too close to each other (and parallel !) then SphericalGeometry has an issue with precision.
Is this related to the ticket #297 ???

[JanPolcher]

I have identified the point in the code which is related to the issue I face

        # Any nodes whose Cartesian coordinates are closer together
        # than 2 ** -32 will cause numerical problems in the
        # intersection calculations, so we merge any nodes that
        # are closer together than that.

2 ** -32 (about 2.3 10^-10) seems to be too generous !
When I go to 2 ** -26 (1.5 10^-8) my example works.
Relevant here is that I am running in numpy.float64. Obviously this should be modulated by the type of reel we are dealing with.

[mcara]

I have no problem running your example:

Area of intersection : 5.90753739260208e-08

and:

[mcara]

I actually was able to reproduce your issue with the released version 1.3.3. Could you please update to 1.3.5 and test again?

[JanPolcher]

Thanks for looking into it. I just updated my copy to the HEAD (commit 0469ca0 (HEAD -> main, origin/main, origin/HEAD)) and I still have the same problem.
It is stuck in computing the area of the overlap polygon.
But the reduction of criteria for the proximity of two nodes still works.
Probably this criteria needs to be a function of the precision of the machine and python libraries. In FORTRAN there is an EPSILON which is the smallest number which can be represented for the type considered. Something like that would work.