CODEWITHSUNDEEP
pythonscipypolygonshapely
Mohamad
Mohamad

Reputation: 33

"ConvexHull", "alpha_shape" in python generate different area values why?

I am trying to compute the area of an irregular polygon.

I used "ConvexHull" and "alpha_shape" but the output of the surface(area) is not the same.

any suggestions please?

you will find the data, and the code used to compute the area below:

###################################################################

data = [(1.866200,4.379100),
(1.128700,4.814700),
(2.036000,3.074600),
(1.077500,6.100000),
(1.833300,7.126600),
(2.061200,8.399900),
(2.710500,9.418000),
(3.426000,9.141000),
(2.795400,2.644300),
(2.773300,1.400700),
(3.554100,1.084300),
(4.338900,2.025600),
(4.365600,9.981400),
(5.284800,9.657300),
(6.305400,10.194400),
(6.829000,9.124800),
(5.889000,4.120900),
(4.884500,1.854300),
(5.694300,2.821400),
(6.825100,6.487400),
(6.698000,5.189600),
(7.676500,8.884200),
(7.670500,7.610700)]

####################################################

import numpy as np

from shapely.geometry import Polygon

from scipy.spatial import ConvexHull, convex_hull_plot_2d

####################################################

hull = ConvexHull(data)

print('area of the surface = ', hull.volume)

####################################################

polygon = Polygon(data)

print('area of the surface = ', polygon.area)

Upvotes: 1

Views: 1244

Answers (2)

Mad Physicist
Mad Physicist

Reputation: 114330

If your points are well behaved and you want to get a convex hull that has approximately the same area as a polygon through the same points, you need to reorder the points. For example, you could order them by angle about the barycenter:

data = np.array(data)
center = data.mean(0)
angle = np.arctan2(*(data - center).T[::-1])
index = np.argsort(angle)
polygon2 = Polygon(data[index])

The plot looks much better now:

fig, ax = plt.subplots(1, 2)

ax[0].set_title("Original")
convex_hull_plot_2d(hull, ax=ax[0])
ax[0].plot(*polygon.exterior.xy, 'r--')

ax[1].set_title("Untangled")
convex_hull_plot_2d(hull, ax=ax[1])
ax[1].plot(*polygon2.exterior.xy, 'r--')
ax[1].plot(*center, 'kx')

plt.show()

enter image description here

And indeed you now get that the hull is only about 10% larger than the polygon:

>>> polygon2.area
37.25662339499999
>>> hull.volume
41.115790855
>>> polygon2.area / hull.volume * 100
90.61390434247064
>>> (hull.volume / polygon2.area - 1) * 100
10.35833929200876

Upvotes: 2

Mr. T
Mr. T

Reputation: 12410

Your shapely polygon is not a hull figure (it is not even convex), so unsurprisingly its area is lower:

import numpy as np    
from shapely.geometry import Polygon    
from scipy.spatial import ConvexHull, convex_hull_plot_2d
from matplotlib import pyplot as plt
   
hull = ConvexHull(data)    
print('area of the surface = ', hull.volume)
    
polygon = Polygon(data)    
print('area of the surface = ', polygon.area)

fig, ax = plt.subplots()
convex_hull_plot_2d(hull, ax=ax)
x, y = polygon.exterior.xy
ax.plot(x, y, c="red", ls="--", label="shapely")
ax.legend()
plt.show()

enter image description here

If we sort the data points according to the convex hull,

data_shapely = [
    (2.710500,9.418000),
    (4.365600,9.981400),
    (6.305400,10.194400),
    (7.676500,8.884200),
    (7.670500,7.610700),
    (6.698000,5.189600),
    (5.694300,2.821400),
    (4.884500,1.854300),
    (3.554100,1.084300),
    (2.773300,1.400700),
    (1.128700,4.814700),
    (1.077500,6.100000),
    (2.061200,8.399900)    
    ]

polygon = Polygon(data_shapely)

then the areas are the same:

area of the surface =  41.115790855
area of the surface =  41.11579085500001

You also misunderstand how shapely calculates areas of non-convex polygons.

data_shapely = [
    (0,0),
    (0,1),
    (1,1),
    (1,0)    
    ]

has unsurprisingly an area of 1. One might assume that

data_shapely = [
    (0,0),
    (1,1),
    (0,1),
    (1,0)    
    ]

has an area of 0.5, given its form. The result, however, is 0, as seemingly "flipped" areas in non-convex polygons are counted as negative values.

Upvotes: 2

Related Questions