Home , Post

Area calculation

Nikolai Shokhirev http://www.numericalexpert.com/blog/area_calculation/
In [1]:
import numpy as np

$$A=\sum_{i=1}^{N-1}I_{i+1,i}+I_{1,N}$$where e.g.

$$I_{2,1}=\intop_{x_{1}}^{x_{2}}ydx=\frac{1}{2}(y_{2}+y_{1})(x_{2}-x_{1})$$

General implementation

In [2]:
def poly_area(xy):
    """ 
        Calculates polygon area.
        x = xy[:,0], y = xy[:,1]
    """
    l = len(xy)
    s = 0.0
    # Python arrys are zer0-based
    for i in range(l):
        j = (i+1)%l  # keep index in [0,l)
        s += (xy[j,0] - xy[i,0])*(xy[j,1] + xy[i,1])
    return -0.5*s

Polygon plotting

In [3]:
def poly_plot(xy, titlestr = "", margin = 0.25):
    """ 
        Plots polygon. For arrow see:
        http://matplotlib.org/examples/pylab_examples/arrow_simple_demo.html
        x = xy[:,0], y = xy[:,1]
    """
    xmin = np.min(xy[:,0])
    xmax = np.max(xy[:,0])
    ymin = np.min(xy[:,1])
    ymax = np.max(xy[:,1])
    hl = 0.1
    l = len(xy)
    for i in range(l):
        j = (i+1)%l  # keep index in [0,l)
        dx = xy[j,0] - xy[i,0]
        dy = xy[j,1] - xy[i,1]
        dd = sqrt(dx*dx + dy*dy)
        dx = dx*(1 - hl/dd)
        dy = dy*(1 - hl/dd)
        arrow(xy[i,0], xy[i,1], dx, dy, head_width=0.05, head_length=0.1, fc='b', ec='b')
        xlim(xmin-margin, xmax+margin)
        ylim(ymin-margin, ymax+margin)
    title(titlestr)

3 x 2 Rectangle

In [4]:
xy = np.array([[1,1],[4,1],[4,3],[1,3]])
poly_plot(xy,'rectangle')
In [5]:
print(poly_area(xy))

6.0

Python-specific implementation

In [6]:
def poly_area1(xy):
    """ 
        Calculates polygon area.
        x = xy[:,0], y = xy[:,1]
    """
    xy1 = np.roll(xy,-1,axis = 0) # shift by -1
    return -0.5*np.inner(xy1[:,0] - xy[:,0],xy1[:,1] + xy[:,1])
In [7]:
print(poly_area1(xy))

6.0

Self-intersecting polygon

In [8]:
uv = np.array([[1,1],[4,3],[4,1],[1,3]])
poly_plot(uv,'Self-intersecting polygon')
In [9]:
print(poly_area(uv))

-0

More complex polygon

In [10]:
uv1 = np.array([[1,1],[2.5,2],[4,1],[4,3],[2.5,2],[1,3]])
poly_plot(uv1,'polygon')
In [11]:
print(poly_area(uv1))

3.0

$A(uv_1) = \frac{1}{2} A(uv)$

Polygon perimeter

In [12]:
def poly_perimeter(xy):
    """ 
        Calculates polygon perimeter.
        x = xy[:,0], y = xy[:,1]
    """
    xy1 = np.roll(xy,-1,axis = 0) # shift by -1
    return np.sum(sqrt((xy1[:,0] - xy[:,0])**2 + (xy1[:,1] - xy[:,1])**2))
In [13]:
print(poly_perimeter(xy))

10.0
In [14]:
print(poly_perimeter(uv))

11.2111025509

$P(uv) = 2*(2+\sqrt{13}) = 11.211102550927978586 $

In [15]:
print(poly_perimeter(uv1))

11.2111025509

Implementation using numpy.apply_along_axis

In [16]:
def poly_perimeter2(xy):
    """ 
        Calculates polygon perimeter.
        x = xy[:,0], y = xy[:,1]
    """
    xy1 = np.roll(xy,-1,axis = 0) # shift by -1
    return np.sum(apply_along_axis(linalg.norm, 1, xy1-xy))
In [17]:
print(poly_perimeter2(xy))

10.0
In [18]:
print(poly_perimeter2(uv1))

11.2111025509

 

Home , Post