Finding an area is a common task in GIS. Using a GIS software like QGIS we can find area of a polygon using a measurement tool or calculate area of each polygon feature with field calculator tool. It's so convenient and that's the reason why we use the software. On the other hand, using a software to find a polygon area, in particular irregular polygon area is like a black box method. It means we can get a result without knowing how it works. This post is like unboxing the method. I will explain how to calculate area of polygon from a given unordered coordinate points.
Calculate Area of Polygon with an Ordered Coordinate Points
Let's start with finding area of an ordered polygon points. Suppose we have a polygon as in figure 1. From the figure we can see the ordered point from 1 to 4 with corresponding coordinates (2,2),(2,4),(5,4) and (9,2).
 |
| Figure 1. A polygon and ordered points |
Given an ordered points, we can calculate the area of polygon with the famous Shoelace method as in figure 2.
 |
| Figure 2. Shoelace Method |
Using the Shoelace method, let's calculate the area of polygon in figure 1 as following. \[A=\bigg| \frac{(2.4+2.4+5.2+9.2)-(2.2+4.5+4.9+2.2)}{2}\bigg |=\bigg |\frac{44-64}{2}\bigg |=10 \]
We get the area of polygon 10 in an area unit. Further more if you want to prove the result, simply divide the polygon into a rectangle and triangle along a dash line as in figure 1. Compute the area for each geometry and add both of them. The result will be the same.
Implementation in Python
We already calculated the area of polygon manually using the Shoelace method. Now let's do it in Python.
The following code is the code for computing a polygon area with a given ordered points. The code consist of two functions which is called explode_xy and shoelace_area. The first function is used to explode a given point coordinates into x and y list. The both list will be as input argument for shoelace_area function to calculate a polygon area using formula as in figure 2.
#ORDERED POINTS DATA xy=[(2,4),(2,2),(9,2),(5,4)] #EXPLODE X AND Y def explode_xy(xy): xl=[] yl=[] for i in range(len(xy)): xl.append(xy[i][0]) yl.append(xy[i][1]) return xl,yl def shoelace_area(x_list,y_list): a1,a2=0,0 x_list.append(x_list[0]) y_list.append(y_list[0]) for j in range(len(x_list)-1): a1 += x_list[j]*y_list[j+1] a2 += y_list[j]*x_list[j+1] l=abs(a1-a2)/2 return l xy_e=explode_xy(xy) A=shoelace_area(xy_e[0],xy_e[1]) print(A)
How to Sort a Random Polygon Points
Now, suppose you build a robot that can measure coordinates in a room with such an integrated sensor. Most probably the recorded coordinates not stored in an ordered structure. You want calculate the area of the room with the Shoelace method cause it's quite straight forward and easy to implement. But the problem is the method only works with an ordered points. So the question is: how to sort a random unordered polygon points?
Might be there are numerous approaches to tackle this problem, but here is mine. The main idea is to order the points based on the angle from a reference axis as illustrated in figure 3. From the figure can be seen that point p has the smallest angle and point s has the biggest angle. Point q's angle is bigger than p and r's angle is bigger than q. Overall the order of points will be p-q-r-s.
 |
| Figure 3. Angle of points from a reference axis |
Following is step by step approach to order a random polygon points.
- Find a centroid coordinate by taking the mean of x and y.
- Subtract each points coordinate with centroid coordinate.
- Put of subtracted coordinates into a corresponding quadrant.
- Calculate the angle of each point with arc tangent.
- Adjust the calculated angles with respected to quadrant and reference axis.
- Order the points based on the adjusted angle from the smallest to the largest.
Implementing The Sorting Algorithm in Python
Based on the steps above, let's implement it in python.
Calculate Centroid Coordinate
To calculate the centroid coordinate can be done by taking the mean of x and y. From the previous code we split each coordinate with explode_xy function. The function returns x and y list. Then to calculate the mean, simply we sum up the list and divide it with the length of the list.
The following code is to calculate the centroid coordinate. At the First line the numpy library is imported. This library will be used at later steps.
import numpy as np xl=xy_e[0] yl=xy_e[1] xc=sum(xl)/len(xl) yc=sum(yl)/len(yl)
Subtract Each Coordinate with Centroid
To subtract each point coordinate with centroid can be done easly using numpy array. For that we need to transform the x list and y list into numpy array and then subtract it with centroid coordinate as the following code.
xa=np.array(xl)-xc ya=np.array(yl)-yc
Grouping Subtracted Coordinates Into a Quadrant
At this step we will group each coordinate into a corresponding quadrant. Why do we need this step? Because at the later step when calculating the angle of each point, the reference axis is not the same. And at the end we have to reference all the angles from the same axis. When doing the transformation the condition will be difference for each quadrant.
In this tutorial I determined the first quadrant form the bottom left to the bottom right with clockwise direction as shown in figure 4.
 |
| Figure 4. Quadrant |
Each point will be grouped into a corresponding quadrant based on its x and y value. From figure 4 can be seen that, if a point has x and y negative, it will be in the first quadrant. If x negative but y positive then it will be in the second quadrant and so on.
The following code is used to group a point coordinate based on the condition above. The result is stored into a list which consist of four dictionaries data structure for each corresponding quadrant. The dictionary structure is used to map each point based on it's index in the unordered coordinate list. Therefore the index will be the key for each values in a dict items.
q=[{},{},{},{}] for i in range(len(xa)): if (xa[i]<0 and ya[i]<0): q[0][i]=(xa[i],ya[i]) elif (xa[i]<0 and ya[i]>0): q[1][i]=(xa[i],ya[i]) elif (xa[i]>0 and ya[i]>0): q[2][i]=(xa[i],ya[i]) else: q[3][i]=(xa[i],ya[i])
Calculating and Adjusting Angle
To calculate a point angle from the centroid can be used arc tangent function by dividing y with x. The resulting angle will have different reference axis for different quadrant ($\beta_i$), but mainly along x axis as illustrated in figure 5. What we need is the angle from a common reference axis, which is -y axis (the highlighted blue line). To get the angle from the common reference axis ($\alpha_i$), it need to be adjusted with adding or subtract with $90^o$ or $270^o$.
 |
| Figure 5. Point angle and adjustment |
The following code is to calculate the angle of points in each quadrant and perform the required adjustment.
alpha={} for i in range(len(q)): for j in q[i].keys(): beta=np.degrees(np.arctan(abs(q[i][j][1])/abs(q[i][j][0]))) if i==0: s=90-beta elif i==1: s=90+beta elif i==2: s=270-beta else: s=270+beta alpha[j]=s
Sorting The Angles
From the previous step, we get a variable alpha which is a dictionary that contain adjusted angle from a common reference axis. The last step is to sort the angle from smallest to largest one. This can be done with the next code.
The first line of the code is declaring an empty list to store the sorted angle. Then the angle is sorted with the sorted built in function. Next in a loop, the point coordinate will be stored into the list which has equal key.
re_xy=[] sorted_alpha=sorted(alpha.values()) for i in sorted_alpha: for k,l in alpha.items(): if i==l: check_p=xy[k] in re_xy if not check_p: re_xy.append(xy[k])
After running the code, you should get an ordered polygon coordinate in a list variable which is called re_xy. To calculate the area using the Shoelace method, explode the list into x and y again and use shoelace_area function as in the first listing code above.
I've tested the algorithm for several polygon shapes and it works well. If you find it not works for a certain shape, please provide me the coordinates of the shape in sequence order, so I can check it. I also very welcome if you have any suggestion to improve the algorithm.
That's all this tutorial how to calculate a polygon area from unordered coordinate points. In this tutorial we discuss about the famous shoelace method to calculate area from points coordinate. Then construct an algorithm to sort coordinates based on the angle from centroid coordinate. The algorithm is explaining step by step with implementation in Python code. Hope it useful for you and thanks for reading!
