My inclination would be to simply start slicing off triangles. I don't see how anything else could avoid being awfully hairy.

Take three sequential points that comprise the polygon. Ensure the angle is less than 180. You now have a new triangle which should be no problem to calculate, delete the middle point from the polygon's list of points. Repeat until you have only three points left.

One way to do it would be to decompose the polygon into triangles, compute the area of the triangles, and take the sum as the area of the polygon.

1. Set a base point (the most convex point). This will be your pivot point of the triangles.
2. Calculate the most-left point (arbitrary), other than your base point.
3. Calculate the 2nd-most-left point to complete your triangle.
4. Save this triangulated area.
5. Shift over one point to the right each iteration.
6. Sum the triangulated areas

Or do a contour integral. Stokes' Theorem allows you to express an area integral as a contour integral. A little Gauss quadrature and Bob's your uncle.

A set of points without any other constraints don't necessarily uniquely define a polygon.

So, first you have to decide what polygon to build from these points - perhaps the convex hull? http://en.wikipedia.org/wiki/Convex_hull

Then triangulate and calculate area. http://www.mathopenref.com/polygonirregulararea.html

Here is the standard method, AFAIK. Basically sum the cross products around each vertex. Much simpler than triangulation.

Python code, given a polygon represented as a list of (x,y) vertex coordinates, implicitly wrapping around from the last vertex to the first:

def area(p):
return 0.5 * abs(sum(x0*y1 - x1*y0
for ((x0, y0), (x1, y1)) in segments(p)))


def segments(p):
return zip(p, p[1:] + [p[0]])

David Lehavi comments: It is worth mentioning why this algorithm works: It is an application of Green's theorem for the functions −y and x; exactly in the way a planimeter works. More specifically:

Formula above =
integral_over_perimeter(-y dx + x dy) =
integral_over_area((-(-dy)/dy+dx/dx) dy dx) =
2 Area

To expand on the triangulate and sum triangle areas, those work if you happen to have a convex polygon OR you happen to pick a point that doesn't generate lines to every other point that intersect the polygon.

For a general non-intersecting polygon, you need to sum the cross product of the vectors (reference point, point a), (reference point, point b) where a and b are "next" to each other.

Assuming you have a list of points that define the polygon in order (order being points i and i+1 form a line of the polygon):

Sum(cross product ((point 0, point i), (point 0, point i + 1)) for i = 1 to n - 1.

Take the magnitude of that cross product and you have the surface area.

This will handle concave polygons without having to worry about picking a good reference point; any three points that generate a triangle that is not inside the polygon will have a cross product that points in the opposite direction of any triangle that is inside the polygon, so the areas get summed correctly.

Better than summing triangles is summing trapezoids in the Cartesian space:

area = 0;
for (i = 0; i < n; i++) {
i1 = (i + 1) % n;
area += (vertex[i].y + vertex[i1].y) * (vertex[i1].x - vertex[i].x) / 2.0;
}

The cross product is a classic.

If you have zillion of such computation to do, try the following optimized version that requires half less multiplications:

area = 0;
for( i = 0; i < N; i += 2 )
area += x[i+1]*(y[i+2]-y[i]) + y[i+1]*(x[i]-x[i+2]);
area /= 2;

I use array subscript for clarity. It is more efficient to use pointers. Though good compilers will do it for you.

The polygon is assumed to be "closed", which means you copy the first point as point with subscript N. It also assume the polygon has an even number of points. Append an additional copy of the first point if N is not even.

The algorithm is obtained by unrolling and combining two successive iterations of the classic cross product algorithm.

I'm not so sure how the two algorithms compare regarding numerical precision. My impression is that the above algorithm is better than the classic one because the multiplication tend to restore the loss of precision of the subtraction. When constrained to use floats, as with GPU, this can make a significant difference.

EDIT: "Area of Triangles and Polygons 2D & 3D" describes an even more efficient method

// "close" polygon
x[N] = x[0];
x[N+1] = x[1];
y[N] = y[0];
y[N+1] = y[1];


// compute area
area = 0;
for( size_t i = 1; i <= N; ++i )
area += x[i]*( y[i+1] - y[i-1] );
area /= 2;

language independent solution:

GIVEN: a polygon can ALWAYS be composed by n-2 triangles that do not overlap (n = number of points OR sides). 1 triangle = 3 sided polygon = 1 triangle; 1 square = 4 sided polygon = 2 triangles; etc ad nauseam QED

therefore, a polygon can be reduced by "chopping off" triangles and the total area will be the sum of the areas of these triangles. try it with a piece of paper and scissors, it is best if you ca visualize the process before following.

if you take any 3 consecutive points in a polygons path and create a triangle with these points, you will have one and only one of three possible scenarios:

1. resulting triangle is completely inside original polygon
2. resulting triangle is totally outside original polygon
3. resulting triangle is partially contained in original polygon

we are interested only in cases that fall in the first option (totally contained).

every time we find one of these, we chop it off, calculate its area (easy peasy, wont explain formula here) and make a new polygon with one less side (equivalent to polygon with this triangle chopped off). until we have only one triangle left.

how to implement this programatically:

create an array of (consecutive) points that represent the path AROUND the polygon. start at point 0. run the array making triangles (one at a time) from points x, x+1 and x+2. transform each triangle from a shape to an area and intersect it with area created from polygon. IF the resulting intersection is identical to the original triangle, then said triangle is totally contained in polygon and can be chopped off. remove x+1 from the array and start again from x=0. otherwise (if triangle is outside [partially or completely] polygon), move to next point x+1 in array.

additionally if you are looking to integrate with mapping and are starting from geopoints, you must convert from geopoints to screenpoints FIRST. this requires deciding a modelling and formula for earths shape (though we tend to think of the earth as a sphere, it is actually an irregular ovoid (eggshape), with dents). there are many models out there, for further info wiki. an important issue is whether or not you will consider the area to be a plane or to be curved. in general, "small" areas, where the points are up to some km apart, will not generate significant error if consider planar and not convex.

To calc the area of the polygon

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1#polygon_area

int cross(vct a,vct b,vct c)
{
vct ab,bc;
ab=b-a;
bc=c-b;
return ab.x*bc.y-ab.y*bc.x;
}
double area(vct p[],int n)
{
int ar=0;
for(i=1;i+1<n;i++)
{
vct a=p[i]-p[0];
vct b=p[i+1]-p[0];
area+=cross(a,b);
}
return abs(area/2.0);
}

This page shows that the formula

can be simplified to:

If you write out a few terms and group them according to common factors of xi, the equality is not hard to see.

The final summation is more efficient since it requires only n multiplications instead of 2n.

def area(x, y):
return abs(sum(x[i] * (y[i + 1] - y[i - 1]) for i in xrange(-1, len(x) - 1))) / 2.0

I learned this simplification from Joe Kington, here.

If you have NumPy, this version is faster (for all but very small arrays):

def area_np(x, y):
x = np.asanyarray(x)
y = np.asanyarray(y)
n = len(x)
shift_up = np.arange(-n+1, 1)
shift_down = np.arange(-1, n-1)
return (x * (y.take(shift_up) - y.take(shift_down))).sum() / 2.0

Implementation of Shoelace formula could be done in Numpy. Assuming these vertices:

import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

We can define the following function to find area:

def PolyArea(x,y):
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

And getting results:

print PolyArea(x,y)
# 0.26353377782163534

Avoiding loop makes this function ~50X faster than PolygonArea:

%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop

Note: I had written this answer for another question, I just mention this here to have a complete list of solutions.

C way of doing that:

float areaForPoly(const int numVerts, const Point *verts)
{
Point v2;
float area = 0.0f;


for (int i = 0; i<numVerts; i++){
v2 = verts[(i + 1) % numVerts];
area += verts[i].x*v2.y - verts[i].y*v2.x;
}


return area / 2.0f;
}

I'm going to give a few simple functions for calculating area of 2d polygon. This works for both convex and concave polygons. we simply divide the polygon into many sub-triangles.

//don't forget to include cmath for abs function
struct Point{
double x;
double y;
}
// cross_product
double cp(Point a, Point b){ //returns cross product
return a.x*b.y-a.y*b.x;
}


double area(Point * vertices, int n){  //n is number of sides
double sum=0.0;
for(i=0; i<n; i++){
sum+=cp(vertices[i], vertices[(i+1)%n]); //%n is for last triangle
}
return abs(sum)/2.0;
}

Python code

As described here: http://www.wikihow.com/Calculate-the-Area-of-a-Polygon

With pandas

import pandas as pd


df = pd.DataFrame({'x': [10, 20, 20, 30, 20, 10, 0], 'y': [-10, -10, -10, 0, 10, 30, 20]})
df = df.append(df.loc[0])


first_product = (df['x'].shift(1) * df['y']).fillna(0).sum()
second_product = (df['y'].shift(1) * df['x']).fillna(0).sum()


(first_product - second_product) / 2
600