make polygon valid

[zjhthu]

Please note

I got a polygon but it is invalid and the area is zero. The polygon looks like this,

How can I get a valid polygon? More specifically, a valid rectangle in this case. I have tried several methods but none of them works. The code is attached below:

import numpy as np
import shapely
from shapely.geometry import Polygon, MultiPolygon, MultiPoint
from shapely.validation import make_valid


polygon_points = np.array( [
    [-3.63173005, -3.25323995],
       [-3.63173005, -3.25323994],
       [-3.63173   , -3.25323943],
       [-3.63173   , -3.25323939],
       [-3.63173   ,  0.91996939],
       [-3.63173   ,  0.91996939],
       [-3.63173   ,  0.91996939],
       [-3.63173   , -3.25319091],
       [-3.63173024, -3.25319572],
       [-3.63173094, -3.25320049],
       [-3.63173211, -3.25320516],
       [-3.63173374, -3.2532097 ],
       [-3.6317358 , -3.25321405],
       [-3.63173827, -3.25321818],
       [-3.63174114, -3.25322205],
       [-3.63174438, -3.25322562],
       [-3.63174795, -3.25322886],
       [-3.63175182, -3.25323173],
       [-3.63175595, -3.2532342 ],
       [-3.6317603 , -3.25323626],
       [-3.63176484, -3.25323789],
       [-3.63176951, -3.25323906],
       [-3.63177428, -3.25323976],
       [-3.63177909, -3.25324   ],
       [-5.42948091, -3.25324   ],
       [-5.42948572, -3.25323976],
       [-5.42949049, -3.25323906],
       [-5.42949516, -3.25323789],
       [-5.4294997 , -3.25323626],
       [-5.42950405, -3.2532342 ],
       [-5.42950818, -3.25323173],
       [-5.42951205, -3.25322886],
       [-5.42951562, -3.25322562],
       [-5.42951886, -3.25322205],
       [-5.42952173, -3.25321818],
       [-5.4295242 , -3.25321405],
       [-5.42952626, -3.2532097 ],
       [-5.42952789, -3.25320516],
       [-5.42952906, -3.25320049],
       [-5.42952976, -3.25319572],
       [-5.42953   , -3.25319091],
       [-5.42953   ,  1.17583091],
       [-5.42952976,  1.17583572],
       [-5.42952906,  1.17584049],
       [-5.42952789,  1.17584516],
       [-5.42952626,  1.1758497 ],
       [-5.4295242 ,  1.17585405],
       [-5.42952173,  1.17585818],
       [-5.42951886,  1.17586205],
       [-5.42951562,  1.17586562],
       [-5.42951205,  1.17586886],
       [-5.42950818,  1.17587173],
       [-5.42950405,  1.1758742 ],
       [-5.4294997 ,  1.17587626],
       [-5.42949516,  1.17587789],
       [-5.42949049,  1.17587906],
       [-5.42948572,  1.17587976],
       [-5.42948091,  1.17588   ],
       [-3.63177909,  1.17588   ],
       [-3.63177428,  1.17587976],
       [-3.63176951,  1.17587906],
       [-3.63176484,  1.17587789],
       [-3.6317603 ,  1.17587626],
       [-3.63175595,  1.1758742 ],
       [-3.63175182,  1.17587173],
       [-3.63174795,  1.17586886],
       [-3.63174438,  1.17586562],
       [-3.63174114,  1.17586205],
       [-3.63173827,  1.17585818],
       [-3.6317358 ,  1.17585405],
       [-3.63173374,  1.1758497 ],
       [-3.63173211,  1.17584516],
       [-3.63173094,  1.17584049],
       [-3.63173024,  1.17583572],
       [-3.63173   ,  1.17583091],
       [-3.63173   ,  0.91998061],
       [-3.63173   ,  0.91998061],
       [-3.63173   ,  0.91998061],
       [-3.63173   ,  1.17587939],
       [-3.63173   ,  1.17587944],
       [-3.63173005,  1.17587991],
       [-3.63173005,  1.17587995],
       [-3.6317306 ,  1.17588   ],
       [-3.63173061,  1.17588   ],
       [-5.42952939,  1.17588   ],
       [-5.42952943,  1.17588   ],
       [-5.42952995,  1.17587995],
       [-5.42952996,  1.17587985],
       [-5.42952999,  1.17587955],
       [-5.42953   ,  1.17587933],
       [-5.42953   , -3.25323933],
       [-5.42952999, -3.25323954],
       [-5.42952996, -3.25323988],
       [-5.42952995, -3.25323994],
       [-5.42952989, -3.25323996],
       [-5.42952945, -3.25324   ],
       [-5.42952939, -3.25324   ],
       [-3.63173061, -3.25324   ],
       [-3.63173057, -3.25324   ],
       [-3.63173006, -3.25323995],
       [-3.63173005, -3.25323995]
       ])

polygon = Polygon(polygon_points)
print('polygon area ' + str(polygon.area))


def simplify_polygon(polygon_points, eps=1e-3, mode=1):
    polygon = Polygon(polygon_points)
    
    if mode == 1:
        polygon_new = polygon.buffer(0)
    elif mode == 2:
        polygon_new = shapely.simplify(polygon, eps)
    elif mode == 3:
        polygon_new = polygon.buffer(eps).buffer(-eps)
    elif mode == 4:
        polygon_new = shapely.validation.make_valid(polygon)
        polygon_new = list(polygon_new.geoms)[0]


    if isinstance(polygon_new, MultiPolygon):
        # polygon_new = polygon_new[0]
        print('MultiPolygon areas ' + str([p.area for p in polygon_new.geoms]))
        # get a list of polygons
        polygons = [p for p in polygon_new.geoms]
        return [np.array(p.exterior.coords) for p in polygons]
    else:
        print('Single Polygon area ' + str(polygon_new.area))
        return np.array(polygon_new.exterior.coords)

simplify_polygon(polygon_points, mode=1)

simplify_polygon(polygon_points, mode=2)

simplify_polygon(polygon_points, mode=3)

simplify_polygon(polygon_points, mode=4)

[jorisvandenbossche]

The specific problem seems to be with the corners of the rectangle. Zooming in on the top left one:

So you can see that the actual "area" are only those tiny parts in the corner. I suppose this is also the reason why a typical method of taking the buffer(0) doesn't work. It actually does "work", but it gives you a valid MultiPolygon consisting of those 4 tiny corner polygons.

When we ask for the reason the polygon is invalid, we also get:

>>> shapely.is_valid_reason(polygon)
'Ring Self-intersection[-5.42948091 -3.25324]'

Given the structure of the current polygon you have, and that this doesn't lend itself to the typical methods to turn it into a valid polygon, the easiest way to get the actual rectangular polygon might be to get the bounding box (envelope):

>>> print(shapely.envelope(polygon))
POLYGON ((-5.42953 -3.25324, -3.63173 -3.25324, -3.63173 1.17588, -5.42953 1.17588, -5.42953 -3.25324))

[zjhthu]

What if my data is not rectangular？ e.g., it is a real 'polygon'.

[jorisvandenbossche]

That's very difficult to say in general, as it will probably depend case by case (the algorithms like make_valid can fix certain invalidities, but they also can't fix everything. For example in this case, the algorithm doesn't necessarily know you want the area in the middle).

I might be more useful to check why you got such data in the first place (where do those coordinates come from? Can that be fixed to return correct coordinates instead?)

[zjhthu]

The data comes from a mesh projection. I project a 3d mesh to a plane to get the outlines of the mesh projection. I will check the projection function.