Change of Internal polygon area after intersecting polygon

Question

I am checking the polygon change after the intersection(Using geopandas, gpd.overlay).
In the meantime, it was found that the area of the polygon inside the boundary was changed.
In my opinion, this is an inner shape that is not affected by the intersection, so I think there shouldn't be any change.

import geopandas as gpd
from shapely.geometry import Polygon

#the blue polygon epsg:5174

blue_polygon = Polygon(zip([391056.5249, 391054.4953, 391069.8387, 391071.7003, 391056.5249],
                           [188161.905, 188148.0892, 188145.6781, 188159.7887, 188161.905]))

#the red polygon epsg:5174

red_polygon =  Polygon(zip([391015.86137871415, 391018.6491016407, 391047.46615036216, 391052.876731906, 391053.1540432138, 391146.0636941614, 391134.1924880867, 391115.9424503843, 391062.8365774889, 391058.69462043804, 391015.86137871415],
                           [188027.04303140735, 188119.7504979779, 188158.7155616982, 188211.34962713026, 188214.0109712184, 188207.30340519012, 188085.12044081726, 188044.01717532205, 188053.51934601943, 188018.0247008443, 188027.04303140735]))

#intersect
blue_new = blue_polygon.intersection(red_polygon)

#area
print(blue_polygon.area)
217.47430880034142

print(blue_new.area)
217.4743088003414

#geometry
print(blue_polygon) 
POLYGON ((391056.5249 188161.905, 391054.4953 188148.0892, 391069.8387 188145.6781, 391071.7003 188159.7887, 391056.5249 188161.905))
print(blue_new) 
POLYGON ((391056.5249 188161.905, 391071.7003 188159.7887, 391069.8387 188145.6781, 391054.4953 188148.0892, 391056.5249 188161.905))

print(blue_polygon == blue_new)
False

I can think of a floating point problem, but I'm not sure why.
Why does this problem appear? Is it a bug?

And, what approach would be good for me to pick a polygon without change after intersection?
To solve this problem, I round to the third decimal place.
Then there seems to be no problem with the results.
But I don't know if this is the right approach.

Answer 1

This appears to be the underlying behaviour of GEOS. The resulting polygon is returned with coordinates in a clock-wise order, which is not correct. Exteriors should be CCW and holes should be CW, so it appears that result of an intersection is essentially the hole which blue cut in red.

If you do the same operation using PyGEOS, the result is exactly the same. I don't think there's much you can do about that in Python. I am not sure if this behaviour is expected in GEOS or if it is a bug.

Normally object.almost_equals(other) would resolve floating point issue, but since the orientation is opposite, not even that works for assessment of equality.

To conclude, rounding looks like your best option now.