pyresample icon indicating copy to clipboard operation
pyresample copied to clipboard

intersection calculation never completes, infinite loop with 100% CPU

Open gerritholl opened this issue 4 years ago • 6 comments

Code Sample, a minimal, complete, and verifiable piece of code

This code sample is probably not minimal, but it's as minimal as I've been able to make it so far:

from pyresample.spherical import SphPolygon
from numpy import array
one = SphPolygon(array([[ 2.66335288, 1.46341369], [ 2.03996841, 1.49309089], [ 1.17181593, 1.46915573], [ 0.90909475, 1.43380709], [ 0.78241556, 1.40104485], [ 0.69729649, 1.36648793], [ 0.62513114, 1.3224743 ], [ 0.54830151, 1.24722376], [ 0.50817615, 1.186389 ], [ 0.23115846, 1.19331475], [-0.0367918 , 1.17324673], [-0.26683713, 1.12983734], [-0.45050788, 1.06906915], [-0.5933579 , 0.99637174], [-0.7047496 , 0.91564441], [-0.79300434, 0.82944525], [-0.86436317, 0.73941807], [-0.92325511, 0.64663193], [-0.97278461, 0.55180129], [-1.01513854, 0.45542036], [-1.05187503, 0.35784426], [-1.08411798, 0.25933845], [-1.11268612, 0.16010919], [-1.13817834, 0.06032253], [-1.16103001, -0.03988369], [-1.1815493 , -0.14039233], [-1.19993916, -0.24110214], [-1.21630767, -0.34192473], [-1.23066743, -0.44278231], [-1.24292205, -0.5436059 ], [-1.25283428, -0.64433317], [-1.25996324, -0.74490528], [-1.2635441 , -0.84526134], [-1.26224861, -0.94532789], [-1.25367377, -1.04499673], [-1.23313148, -1.14407292], [-1.1903402 , -1.2421318 ], [-1.09835852, -1.33802836], [-0.86484729, -1.42754387], [-0.12516873, -1.48708913], [ 0.85883216, -1.45406825], [ 1.20132521, -1.36982862], [ 1.32458409, -1.27528245], [ 1.3796474 , -1.17773082], [ 1.40617791, -1.07888876], [ 1.41813774, -0.97933736], [ 1.42156659, -0.87932642], [ 1.41945451, -0.77898624], [ 1.41341849, -0.67839964], [ 1.40438013, -0.57763123], [ 1.39287325, -0.47674082], [ 1.37919527, -0.37579029], [ 1.36348551, -0.27484729], [ 1.3457657 , -0.17398771], [ 1.32595882, -0.07329748], [ 1.30389374, 0.02712538], [ 1.27929839, 0.12716808], [ 1.25178184, 0.22669887], [ 1.220803 , 0.32556065], [ 1.18562164, 0.42356057], [ 1.14522367, 0.52045375], [ 1.09820814, 0.61591714], [ 1.04261663, 0.70950724], [ 0.97567814, 0.80059109], [ 0.89344091, 0.88823199], [ 0.79029754, 0.97100115], [ 0.65858077, 1.04667881], [ 0.48898223, 1.11184127], [ 0.27383753, 1.16152666], [ 0.0158546 , 1.1897051 ], [-0.26229105, 1.19159251], [-0.52485832, 1.16682671], [-0.74640942, 1.11975449], [-0.92191061, 1.05642095], [-1.0582518 , 0.98199522], [-1.16481114, 0.90009946], [-1.24954095, 0.81309195], [-1.31831344, 0.72248881], [-1.37527748, 0.6292802 ], [-1.42334132, 0.53413182], [-0.03925001, 1.19227661], [-0.01881986, 1.2542977 ], [ 0.02119208, 1.33180974], [ 0.06015454, 1.37794658], [ 0.10834974, 1.41501064], [ 0.18604488, 1.45150917], [ 0.3825057 , 1.49454965], [ 1.74820852, 1.53115399], [ 2.54470923, 1.48639946], [-1.9512296 , 0.61449627], [-1.94217602, 0.71548039], [-1.93622669, 0.81633576], [-1.93443368, 0.91699583], [-1.93868608, 1.01737187], [-1.95260477, 1.11731747], [-1.98390632, 1.21653505], [-2.05196883, 1.31427217], [-2.21983968, 1.4079696 ], [-2.77186451, 1.48364337], [ 2.35648575, 1.47571327], [ 1.89765851, 1.39514712], [ 1.75206211, 1.30051232], [ 1.69135059, 1.20244759], [ 1.66328334, 1.10305337], [ 1.65107612, 1.00297988], [ 1.64785574, 0.90249346], [ 1.65033852, 0.80172946], [ 1.65679154, 0.70077285], [ 1.66624369, 0.59968947], [ 1.67813872, 0.49854002], [ 1.69216801, 0.39738696], [ 1.70818546, 0.29629819], [ 1.72616348, 0.1953494 ], [ 1.74617212, 0.09462581], [ 1.76837282, -0.005776 ], [ 1.79302351, -0.10574562], [ 1.8204943 , -0.20515504], [ 1.85129542, -0.30385276], [ 1.88612129, -0.40165441], [ 1.92591725, -0.49832794], [ 1.97197994, -0.59357013], [ 2.02610762, -0.68696888], [ 2.09082321, -0.77794226], [ 2.16969495, -0.86563872], [ 2.26775487, -0.94877419], [ 2.39189047, -1.02537412], [ 2.55066221, -1.09240844], [ 2.75198818, -1.14544408], [ 2.99627758, -1.17884979], [-3.01632073, -1.18751634], [-2.75173224, -1.16985335], [-2.52151949, -1.12898792], [-2.33561891, -1.07052734], [-2.18990409, -0.99978609], [-2.07573196, -0.92068698], [-1.98501822, -0.83584712], [-1.9115487 , -0.74696423], [-1.85085136, -0.65514764], [-1.79976449, -0.56113971], [-1.756048 , -0.4654535 ], [-1.71809906, -0.36845673], [-1.68475714, -0.27042338], [-1.65517344, -0.17156563], [-1.62872363, -0.0720539 ], [-1.60495 , 0.02797055], [-1.58352355, 0.12838757], [-1.56422059, 0.22909311], [-1.54691068, 0.32999596], [-1.5315552 , 0.43101553], [-1.51821788, 0.53208015], [-1.507092 , 0.63312519], [-1.49855512, 0.7340904 ], [-1.49327489, 0.83491541], [-1.49241943, 0.93553111], [-1.4981059 , 1.03584173], [-1.51445846, 1.13568298], [-1.5504898 , 1.23470789], [-1.62973743, 1.33199796], [-1.83364104, 1.4241953 ], [-2.53106782, 1.49123105]]), radius=1)
other = SphPolygon( array([[-1.66503519, 1.3781822 ], [-1.65506795, 1.38197753], [-1.55416016, 1.41282574], [-1.40759998, 1.44151391], [-1.18656988, 1.46615086], [-0.85843546, 1.48316505], [-0.44132863, 1.4876831 ], [-0.04924721, 1.47785384], [ 0.23522308, 1.45736614], [ 0.42354922, 1.43086215], [ 0.54987815, 1.40118762], [ 0.63842834, 1.36979805], [ 0.64654151, 1.36627686], [ 0.58529545, 1.35999458], [ 0.44817004, 1.34115071], [ 0.33318 , 1.31862022], [ 0.23774442, 1.29326915], [ 0.15862172, 1.2657824 ], [ 0.09273496, 1.23668003], [ 0.03747186, 1.20635157], [-0.00927492, 1.17508948], [-0.04916714, 1.14311584], [-0.08350472, 1.11060209], [-0.11330436, 1.07768306], [-0.13936472, 1.04446697], [-0.16231727, 1.01104231], [-0.17053857, 0.99796649], [-0.17798056, 0.99929745], [-0.24129604, 1.00932545], [-0.30669568, 1.01736752], [-0.37377995, 1.02331562], [-0.44207279, 1.02708584], [-0.51103942, 1.02862315], [-0.5801109 , 1.02790469], [-0.64871283, 1.02494118], [-0.7162949 , 1.01977625], [-0.78235727, 1.01248377], [-0.84647102, 1.00316366], [-0.8533889 , 1.00201967], [-0.86135363, 1.01554005], [-0.88307817, 1.0492582 ], [-0.90780256, 1.08280163], [-0.93615495, 1.11608778], [-0.96893775, 1.1490156 ], [-1.00718497, 1.18145865], [-1.05223883, 1.21325519], [-1.10584836, 1.2441936 ], [-1.17028672, 1.27399133], [-1.2484639 , 1.30226485], [-1.34396043, 1.3284894 ], [-1.46079655, 1.35195201], [-1.60257299, 1.37171678]]), radius=1)
print(one.intersection(other))                                                                                           

Problem description

The code apparently never completes and uses 100% CPU (I interrupted it after several minutes).

Expected Output

Expecting a result within a fraction of a second.

Actual Result, Traceback if applicable

As is there's not much useful information, but adding a logging message in spherical.SphPolygon._bool_oper and then running the MCVE reveals it appears to be in an infinite loop here.

https://github.com/pytroll/pyresample/blob/c440d1b1da71e94171d856af7990895ccd49e646/pyresample/spherical.py#L458-L460

Result with this addition:

$ python hanging-intersection.py
[DEBUG: 2021-02-10 10:32:27 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:27 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:27 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:27 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:28 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:28 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:29 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:29 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:29 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:32:30 : pyresample.spherical] still finding intersection
... (snip)
[DEBUG: 2021-02-10 10:34:50 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:34:51 : pyresample.spherical] still finding intersection
[DEBUG: 2021-02-10 10:34:51 : pyresample.spherical] still finding intersection
... (interrupt)
^CTraceback (most recent call last):
  File "hanging-intersection.py", line 218, in <module>
    print(one.intersection(other))
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/pyresample/spherical.py", line 495, in intersection
    return self._bool_oper(other, -1)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/pyresample/spherical.py", line 476, in _bool_oper
    inter, edge2 = edge1.get_next_intersection(narcs2, inter)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/pyresample/spherical.py", line 305, in get_next_intersection
    inter = self.intersection(arc)
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/pyresample/spherical.py", line 294, in intersection
    (abs(a__.hdistance(i) + b__.hdistance(i) - ab_) < EPSILON)) and
  File "/data/gholl/miniconda3/envs/py38b/lib/python3.8/site-packages/pyresample/spherical.py", line 87, in hdistance
    np.cos(point.lat) * np.cos(self.lat) *
KeyboardInterrupt

Versions of Python, package at hand and relevant dependencies

pyresample 1.17.0+13.g270b5f0 on Python 3.8.6.

Additional context

This is a cause (probably the cause) for the problem described at https://github.com/pytroll/trollflow2/issues/100

gerritholl avatar Feb 10 '21 09:02 gerritholl

The affected code seems to be originally added to pytroll-schedule by @mraspaud 7 years ago in 2014.

gerritholl avatar Feb 10 '21 09:02 gerritholl

could you plot the polygons to see what we are faced with?

mraspaud avatar Feb 10 '21 10:02 mraspaud

One polygon is normal (it's the area of interest), the other is all across the globe and then intersecting itself. It's a long swatch. Background: due to a long timeout in gatherer, granules from subsequent overpasses were grouped together and passed to trollflow2 (which in turn caused the hang in https://github.com/pytroll/trollflow2/issues/100).

polygons

Source code for plotting:

from pyresample.spherical import SphPolygon                                                                                                                                                                                                                   
import cartopy.crs as ccrs                                                                                                                                                                                                                                    
import matplotlib.pyplot as plt                                                                                                                                                                                                                               
import sattools.io                                                                                                                                                                                                                                            
import typhon.plots.common                                                                                                                                                                                                                                    
from satpy.utils import debug_on; debug_on()                                                                                                                                                                                                                  
from numpy import array, rad2deg                                                                                                                                                                                                                              
one = SphPolygon(array([[ 2.66335288, 1.46341369], [ 2.03996841, 1.49309089], [ 1.17181593, 1.46915573], [ 0.90909475, 1.43380709], [ 0.78241556, 1.40104485], [ 0.69729649, 1.36648793], [ 0.62513114, 1.3224743 ], [ 0.54830151, 1.24722376], [ 0.50817615, 1.186389 ], [ 0.23115846, 1.19331475], [-0.0367918 , 1.17324673], [-0.26683713, 1.12983734], [-0.45050788, 1.06906915], [-0.5933579 , 0.99637174], [-0.7047496 , 0.91564441], [-0.79300434, 0.82944525], [-0.86436317, 0.73941807], [-0.92325511, 0.64663193], [-0.97278461, 0.55180129], [-1.01513854, 0.45542036], [-1.05187503, 0.35784426], [-1.08411798, 0.25933845], [-1.11268612, 0.16010919], [-1.13817834, 0.06032253], [-1.16103001, -0.03988369], [-1.1815493 , -0.14039233], [-1.19993916, -0.24110214], [-1.21630767, -0.34192473], [-1.23066743, -0.44278231], [-1.24292205, -0.5436059 ], [-1.25283428, -0.64433317], [-1.25996324, -0.74490528], [-1.2635441 , -0.84526134], [-1.26224861, -0.94532789], [-1.25367377, -1.04499673], [-1.23313148, -1.14407292], [-1.1903402 , -1.2421318 ], [-1.09835852, -1.33802836], [-0.86484729, -1.42754387], [-0.12516873, -1.48708913], [ 0.85883216, -1.45406825], [ 1.20132521, -1.36982862], [ 1.32458409, -1.27528245], [ 1.3796474 , -1.17773082], [ 1.40617791, -1.07888876], [ 1.41813774, -0.97933736], [ 1.42156659, -0.87932642], [ 1.41945451, -0.77898624], [ 1.41341849, -0.67839964], [ 1.40438013, -0.57763123], [ 1.39287325, -0.47674082], [ 1.37919527, -0.37579029], [ 1.36348551, -0.27484729], [ 1.3457657 , -0.17398771], [ 1.32595882, -0.07329748], [ 1.30389374, 0.02712538], [ 1.27929839, 0.12716808], [ 1.25178184, 0.22669887], [ 1.220803 , 0.32556065], [ 1.18562164, 0.42356057], [ 1.14522367, 0.52045375], [ 1.09820814, 0.61591714], [ 1.04261663, 0.70950724], [ 0.97567814, 0.80059109], [ 0.89344091, 0.88823199], [ 0.79029754, 0.97100115], [ 0.65858077, 1.04667881], [ 0.48898223, 1.11184127], [ 0.27383753, 1.16152666], [ 0.0158546 , 1.1897051 ], [-0.26229105, 1.19159251], [-0.52485832, 1.16682671], [-0.74640942, 1.11975449], [-0.92191061, 1.05642095], [-1.0582518 , 0.98199522], [-1.16481114, 0.90009946], [-1.24954095, 0.81309195], [-1.31831344, 0.72248881], [-1.37527748, 0.6292802 ], [-1.42334132, 0.53413182], [-0.03925001, 1.19227661], [-0.01881986, 1.2542977 ], [ 0.02119208, 1.33180974], [ 0.06015454, 1.37794658], [ 0.10834974, 1.41501064], [ 0.18604488, 1.45150917], [ 0.3825057 , 1.49454965], [ 1.74820852, 1.53115399], [ 2.54470923, 1.48639946], [-1.9512296 , 0.61449627], [-1.94217602, 0.71548039], [-1.93622669, 0.81633576], [-1.93443368, 0.91699583], [-1.93868608, 1.01737187], [-1.95260477, 1.11731747], [-1.98390632, 1.21653505], [-2.05196883, 1.31427217], [-2.21983968, 1.4079696 ], [-2.77186451, 1.48364337], [ 2.35648575, 1.47571327], [ 1.89765851, 1.39514712], [ 1.75206211, 1.30051232], [ 1.69135059, 1.20244759], [ 1.66328334, 1.10305337], [ 1.65107612, 1.00297988], [ 1.64785574, 0.90249346], [ 1.65033852, 0.80172946], [ 1.65679154, 0.70077285], [ 1.66624369, 0.59968947], [ 1.67813872, 0.49854002], [ 1.69216801, 0.39738696], [ 1.70818546, 0.29629819], [ 1.72616348, 0.1953494 ], [ 1.74617212, 0.09462581], [ 1.76837282, -0.005776 ], [ 1.79302351, -0.10574562], [ 1.8204943 , -0.20515504], [ 1.85129542, -0.30385276], [ 1.88612129, -0.40165441], [ 1.92591725, -0.49832794], [ 1.97197994, -0.59357013], [ 2.02610762, -0.68696888], [ 2.09082321, -0.77794226], [ 2.16969495, -0.86563872], [ 2.26775487, -0.94877419], [ 2.39189047, -1.02537412], [ 2.55066221, -1.09240844], [ 2.75198818, -1.14544408], [ 2.99627758, -1.17884979], [-3.01632073, -1.18751634], [-2.75173224, -1.16985335], [-2.52151949, -1.12898792], [-2.33561891, -1.07052734], [-2.18990409, -0.99978609], [-2.07573196, -0.92068698], [-1.98501822, -0.83584712], [-1.9115487 , -0.74696423], [-1.85085136, -0.65514764], [-1.79976449, -0.56113971], [-1.756048 , -0.4654535 ], [-1.71809906, -0.36845673], [-1.68475714, -0.27042338], [-1.65517344, -0.17156563], [-1.62872363, -0.0720539 ], [-1.60495 , 0.02797055], [-1.58352355, 0.12838757], [-1.56422059, 0.22909311], [-1.54691068, 0.32999596], [-1.5315552 , 0.43101553], [-1.51821788, 0.53208015], [-1.507092 , 0.63312519], [-1.49855512, 0.7340904 ], [-1.49327489, 0.83491541], [-1.49241943, 0.93553111], [-1.4981059 , 1.03584173], [-1.51445846, 1.13568298], [-1.5504898 , 1.23470789], [-1.62973743, 1.33199796], [-1.83364104, 1.4241953 ], [-2.53106782, 1.49123105]]), radius=1)
other = SphPolygon( array([[-1.66503519, 1.3781822 ], [-1.65506795, 1.38197753], [-1.55416016, 1.41282574], [-1.40759998, 1.44151391], [-1.18656988, 1.46615086], [-0.85843546, 1.48316505], [-0.44132863, 1.4876831 ], [-0.04924721, 1.47785384], [ 0.23522308, 1.45736614], [ 0.42354922, 1.43086215], [ 0.54987815, 1.40118762], [ 0.63842834, 1.36979805], [ 0.64654151, 1.36627686], [ 0.58529545, 1.35999458], [ 0.44817004, 1.34115071], [ 0.33318 , 1.31862022], [ 0.23774442, 1.29326915], [ 0.15862172, 1.2657824 ], [ 0.09273496, 1.23668003], [ 0.03747186, 1.20635157], [-0.00927492, 1.17508948], [-0.04916714, 1.14311584], [-0.08350472, 1.11060209], [-0.11330436, 1.07768306], [-0.13936472, 1.04446697], [-0.16231727, 1.01104231], [-0.17053857, 0.99796649], [-0.17798056, 0.99929745], [-0.24129604, 1.00932545], [-0.30669568, 1.01736752], [-0.37377995, 1.02331562], [-0.44207279, 1.02708584], [-0.51103942, 1.02862315], [-0.5801109 , 1.02790469], [-0.64871283, 1.02494118], [-0.7162949 , 1.01977625], [-0.78235727, 1.01248377], [-0.84647102, 1.00316366], [-0.8533889 , 1.00201967], [-0.86135363, 1.01554005], [-0.88307817, 1.0492582 ], [-0.90780256, 1.08280163], [-0.93615495, 1.11608778], [-0.96893775, 1.1490156 ], [-1.00718497, 1.18145865], [-1.05223883, 1.21325519], [-1.10584836, 1.2441936 ], [-1.17028672, 1.27399133], [-1.2484639 , 1.30226485], [-1.34396043, 1.3284894 ], [-1.46079655, 1.35195201], [-1.60257299, 1.37171678]]), radius=1)
proj = ccrs.RotatedPole(pole_latitude=45, pole_longitude=180)
(f, ax) = plt.subplots(1, 1, subplot_kw={"projection": proj}, figsize=(18, 12))
ccg = ccrs.Geodetic()
ax.set_global()
ax.gridlines()
ax.coastlines()
ax.plot((deg:=rad2deg(one.vertices))[:, 0], deg[:, 1], 'x-', transform=ccg, color="blue")
ax.plot((deg:=rad2deg(other.vertices))[:, 0], deg[:, 1], 'x-', transform=ccg,
        color="red")
typhon.plots.common.write_multi(
        f,
        sattools.io.plotdir(create=True) / "polygons.png")

gerritholl avatar Feb 10 '21 13:02 gerritholl

Different example, tested with pyresample 1.22.3:

from pyresample.area_config import parse_area_file
from pyresample.boundary import AreaDefBoundary
from satpy.resample import get_area_file
from satpy.utils import debug_on; debug_on()
eurol, germ = parse_area_file(get_area_file(), "eurol", "germ")
AreaDefBoundary(eurol).contour_poly.intersection(AreaDefBoundary(germ).contour_poly)

gerritholl avatar Jan 04 '22 10:01 gerritholl

Even easier test: checking overlap with itself also runs into an infinite loop:

AreaDefBoundary(germ).contour_poly.intersection(AreaDefBoundary(germ).contour_poly)

gerritholl avatar Jan 04 '22 12:01 gerritholl

Recently, our operational geostationary production stalled on two servers simultaneously. The offending polygons in this case:

polygons

The line is coming from trollflow2s get_twilight_poly:

>>> list(poly.edges())
[((1.03176132689787, 0.0), (-0.5390349998970265, 1.2650196509486906)), ((-0.5390349998970265, 1.2650196509486906), (-2.109831326691923, 0.0)), ((-2.109831326691923, 0.0), (2.6025576536927666, -1.2650196509486906)), ((2.6025576536927666, -1.2650196509486906), (1.03176132689787, 0.0))]

gerritholl avatar May 22 '23 10:05 gerritholl