Search code examples
numpypolygonareapoint-in-polygon

How to calculate area of lots of polygon where each polygon has different coordinates


Suppose I know 700 different values of x as well as another 700 different values of y coordinates. Every 7 points from my coordinates can construct a polygon. I can easily pick any formula (here I have used shoelace formula) to calculate the area of any polygon. But it's really tough to calculate area of each polygon in this way, is there any good way by which I can calculate the area of all polygon (for this case the total number of polygon is 100) at a time, e.g. by using any iteration.

A simple code is shown below:

import numpy as np

x = np.arange(0.5, 700)  #create random coordinates
y = np.arange(0.3, 700)  #create random coordinates

#first 7 coordinates that could form a polygon
x_1 = x[0:7]  
y_1 = y[0:7]

#Area of the first polygon
Area = 0.5 * np.array(np.dot(x_1, np.roll(x_1, 1)) - np.dot(y_1, np.roll(y_1, 1)))

Solution

  • Here's a vectorized way to get polygon areas for all -

    X = x.reshape(-1,7)
    Y = y.reshape(-1,7)
    
    Xr = np.roll(X,1,axis=1)
    Yr = np.roll(Y,1,axis=1)
    p1 = np.einsum('ij,ij->i',X,Xr)
    p2 = np.einsum('ij,ij->i',Y,Yr)
    Area_out = (p1-p2)/2
    

    Alternative, use np.matmul to get the sum-redcutions -

    p1 = np.matmul(X[:,None,:],Xr[:,:,None])[:,0,0]
    

    With equivalent @ operator (on Python3.x) -

    p1 = (X[:,None,:]@Xr[:,:,None])[:,0,0]
    

    Similarly, use it for p2.


    Avoid rolling with slicing

    We can make use of slicing to skip the rolling and hence get further perf. boost, like so -

    p1 = np.einsum('ij,ij->i',X[:,:-1],X[:,1:]) + X[:,0]*X[:,-1]
    p2 = np.einsum('ij,ij->i',Y[:,:-1],Y[:,1:]) + Y[:,0]*Y[:,-1]
    

    As explained earlier, we can bring np.matmul/@ operator on the same lines, i.e. -

    p1 = np.matmul(X[:,None,:-1],X[:,1:,None])[:,0,0] + X[:,0]*X[:,-1]
    

    and so on.