`SphericalPolygon.intersection` returns empty region despite `intersects_poly` returns `True`

[hgao-astro]

For example, the following two images overlap with each other and intersects_poly returns True, but intersection returns an empty polygon. Interestingly, when I increase the steps to >=2, the results are correct. Let me know if you need more info to reproduce.

path1 = "xxx.fits"
path2 = "yyy.fits"
hdr1 = fits.getheader(path1)
hdr2 = fits.getheader(path2)
poly1 = SphericalPolygon.from_wcs(hdr1, steps=1)
poly2 = SphericalPolygon.from_wcs(hdr2, steps=1)
poly1.intersects_poly(poly2)    # True
poly_overlap = poly1.intersection(poly2)
print(poly_overlap)                   # []
print(list(poly1.to_radec()))     # [(array([269.2131081 , 268.71714565, 268.71739619, 269.21561182, 269.2131081 ]), array([67.84928366, 67.85019998, 67.95553246, 67.95461129, 67.84928366]))]
print(list(poly2.to_radec()))    # [(array([269.2131081 , 268.80975469, 268.81043559, 269.21564764, 269.2131081 ]), array([67.84928366, 67.85014272, 67.95697488, 67.95611121, 67.84928366]))]

[mcara]

Could you, please, provide files xxx.fits and yyy.fits?

[hgao-astro]

Sorry for the late reply. These were temporary resampled images by SWarp, and I have already lost them. I wonder if the issue can be reproduced from the (ra, dec) list in the example.