Polygon area calculation using Latitude and Longitude generated from Cartesian space and a world file

Question

Given a series of GPS coordinate pairs, I need to calculate the area of a polygon (n-gon). This is relatively small (not larger than 50,000 sqft). The geocodes are created by applying an affine transform with data from a world file.

I have tried to use a two step approach by doing converting the geocodes to cartesian coordinates:

double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );
double yPos = (lat-latAnchor)*( Math.toRadians( 6378137 ) );

then I use a cross product calculation to determine the area.

The issue is that the results are a bit off in accuracy (around 1%). Is there anything I can look into to improve this?

Thanks.

Polo

Answer 1

1% error seems a bit high due to just your approximation. Are you comparing against actual measurements or some ideal calculation? Remember that there is error in the GPS as well that might be contributing.

If you want a more accurate method for doing this there's a good answer at this question. If you're going for a faster way you can use the WGS84 geoid instead of your reference sphere for converting to cartesian coordinates (ECEF). Here's the wiki link for that conversion.

Destroyer

I am comparing against real measurements of known areas. An interesting side note is that if I run the GPS coordinates through a Haversine method I get very accurate distance calculations which yield accurate perimeter values. – Polo

Sorry for the late response, ended up using the WGs84 geoid with the proj4 java library. Worked great, thanks for the help. – Polo

Answer 2

I checked on internet for various polygon area formulas(or code) but did not find any one good or easy to implement.

Now I have written the code snippet to calculate area of a polygon drawn on earth surface. The polygon can have n vertices with each vertex has having its own latitude longitude.

Few Important Points

1. The array input to this function will have "n + 1" elements. The last element will have same values as that of first one.
2. I have written very basic C# code, so that guys can also adapt it in other language.
3. 6378137 is the value of earth radius in metres.

4. The output area will have unit of square metres

   private static double CalculatePolygonArea(IList<MapPoint> coordinates)
   {
       double area = 0;
   
       if (coordinates.Count > 2)
       {
           for (var i = 0; i < coordinates.Count - 1; i++)
           {
               MapPoint p1 = coordinates[i];
               MapPoint p2 = coordinates[i + 1];
               area += ConvertToRadian(p2.Longitude - p1.Longitude) * (2 + Math.Sin(ConvertToRadian(p1.Latitude)) + Math.Sin(ConvertToRadian(p2.Latitude)));
           }
   
           area = area * 6378137 * 6378137 / 2;
       }
   
       return Math.Abs(area);
   }
   
   private static double ConvertToRadian(double input)
   {
       return input * Math.PI / 180;
   }
   

Klarrisa

I tried your code but something is wrong. any ideas? See: code – Sizing

you have put "area = area * R * R / 2;" inside the for loop while it should be outside loop. – Klarrisa

I think you should convert p1.Longitude and p2.Longitude to radians as well. After doing this modification, I got really similar result as I got from google.maps.geometry.spherical.computeArea function. – Resolute

After corrections this seems fine. And is very similar as getGeodesicArea in Open Layers (minus projection part). See: github.com/openlayers/openlayers/blob/v2.13.1/lib/OpenLayers/… – Devine

Answer 3

I am modifying a Google Map so that a user can calculate the area of a polygon by clicking the vertices. It wasn't giving correct areas until I made sure the Math.cos(latAnchor) was in radians first

So:

double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );

became:

double xPos = (lon-lonAnchor)*( 6378137*PI/180 ) )*Math.cos( latAnchor*PI/180 );

where lon, lonAnchor and latAnchor are in degrees. Works like a charm now.

Fowle

Answer 4

The reason for this "1%" discrepancy is The earth is very slightly ellipsoidal so by calculating using a spherical model gives errors typically up to 0.3%, give or take the location.

Rappee

Answer 5

Based on the solution by Risky Pathak here is the solution for SQL (Redshift) to calculate areas for GeoJSON multipolygons (with the assumption that linestring 0 is the outermost polygon)

create or replace view geo_area_area as 
with points as (
    select ga.id as key_geo_area
    , ga.name, gag.linestring
    , gag.position
    , radians(gag.longitude) as x
    , radians(gag.latitude) as y
    from geo_area ga
    join geo_area_geometry gag on (gag.key_geo_area = ga.id)
)
, polygons as (
    select key_geo_area, name, linestring, position 
    , x
    , lag(x) over (partition by key_geo_area, linestring order by position) as prev_x
    , y
    , lag(y) over (partition by key_geo_area, linestring order by position) as prev_y
    from points
)
, area_linestrings as (
    select key_geo_area, name, linestring
    , abs( sum( (x - prev_x) * (2 + sin(y) + sin(prev_y)) ) ) * 6378137 * 6378137 / 2 / 10^6 as area_km_squared
    from polygons
    where position != 0
    group by 1, 2, 3
)
select key_geo_area, name
, sum(case when linestring = 0 then area_km_squared else -area_km_squared end) as area_km_squared
from area_linestrings
group by 1, 2
;

Cupulate

Answer 6

Adapted RiskyPathak's snippet to PHP

function CalculatePolygonArea($coordinates) {
    $area = 0;
    $coordinatesCount = sizeof($coordinates);
    if ($coordinatesCount > 2) {
      for ($i = 0; $i < $coordinatesCount - 1; $i++) {
        $p1 = $coordinates[$i];
        $p2 = $coordinates[$i + 1];
        $p1Longitude = $p1[0];
        $p2Longitude = $p2[0];
        $p1Latitude = $p1[1];
        $p2Latitude = $p2[1];
        $area += ConvertToRadian($p2Longitude - $p1Longitude) * (2 + sin(ConvertToRadian($p1Latitude)) + sin(ConvertToRadian($p2Latitude)));
      }
    $area = $area * 6378137 * 6378137 / 2;
    }
    return abs(round(($area)));
}

function ConvertToRadian($input) {
    $output = $input * pi() / 180;
    return $output;
}

// mssing clossing )

Cairistiona

Answer 7

Thank you Risky Pathak!

In the spirit of sharing, here's my adaptation in Delphi:

interface

uses 
  System.Math; 

TMapGeoPoint = record
  Latitude: Double;
  Longitude: Double;
end;


function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;

implementation

function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;
var
  Area: Double;
  i: Integer;
  P1, P2: TMapGeoPoint;
begin
 Area := 0;

 // We need at least 2 points
 if (AGeoPoints.Count > 2) then
 begin
   for I := 0 to AGeoPoints.Count - 1 do
   begin
     P1 := AGeoPoints[i];
     if i < AGeoPoints.Count - 1  then
       P2 := AGeoPoints[i + 1]
     else
       P2 := AGeoPoints[0];
     Area := Area + DegToRad(P2.Longitude - P1.Longitude) * (2 + 
        Sin(DegToRad(P1.Latitude)) + Sin(DegToRad(P2.Latitude)));
    end;

    Area := Area * 6378137 * 6378137 / 2;

  end;

  Area := Abs(Area); //Area (in sq meters)

  // 1 Square Meter = 0.000247105 Acres
  result := Area * 0.000247105;
end;

Entirely

Answer 8

Adapted RiskyPathak's snippet to Ruby

def deg2rad(input)
  input * Math::PI / 180.0
end

def polygone_area(coordinates)
  return 0.0 unless coordinates.size > 2

  area = 0.0
  coor_p = coordinates.first
  coordinates[1..-1].each{ |coor|
    area += deg2rad(coor[1] - coor_p[1]) * (2 + Math.sin(deg2rad(coor_p[0])) + Math.sin(deg2rad(coor[0])))
    coor_p = coor
  }

  (area * 6378137 * 6378137 / 2.0).abs # 6378137 Earth's radius in meters
end

Embezzle

Answer 9

Tried to do this in swift playground and got results that are way off Example coord: (39.58571008386715,-104.94522892318253) that I am plugging into the function

func deg2rad(_ number: Double) -> Double {
        return number * .pi / 180
    }
    
func areaCalc(lat: [Double]?, lon: [Double]?){
    guard let lat = lat,
       let lon = lon
    else { return }
    var area: Double = 0.0
    if(lat.count > 2){
        for i in stride(from: 0, to: lat.count - 1, by: 1) {
            
            let p1lon = lon[i]
            let p1lat = lat[i]
            let p2lon = lon[i+1]
            let p2lat = lat[i+1]
            
            area = area + (deg2rad(p2lon - p1lon)) * (2 + sin(deg2rad(p1lat))) + (sin(deg2rad(p2lat)))
        }
        area = area * 6378137.0 * 6378137.0
        area = abs(area / 2)
    }
}

Wreckage

As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. – Hippy

Answer 10

I'm not sure why, but the results I'm getting using the formula above , are nowhere near the results google maps is returning when measuring the same area.

So digging around, I've came up with this javascript method to calculate a polygon area defined by GPS coordinates:

const PI = Math.PI;
const EARTH_RADIUS_IN_METERS = 6378137;
const EARTH_CIRCUMFERENCE_IN_METERS = 2*EARTH_RADIUS_IN_METERS*PI;

function areaClaculator(points) {
    let area = null;
    if (!isValueEmpty(points) && points.length > 2) {
        let p0 = points[0]
        let newPoints = [];
        for (let i=1; i<points.length; i++) {
            let p = points[i];
            
            let y = (p.lat - p0.lat) / 360 * EARTH_CIRCUMFERENCE_IN_METERS;
            let x = (p.lng - p0.lng) / 360 * EARTH_CIRCUMFERENCE_IN_METERS * Math.cos(rad(p.lat));
            let entry = {};
            entry.x = x;
            entry.y = y;
            newPoints.push(entry);
        }
        
        if (!isValueEmpty(newPoints) && newPoints.length > 1) {
            area = 0;
            for (let i=0;i< newPoints.length - 1; i++) {
                let p1 = newPoints[i];
                let p2 = newPoints[i+1];
                
                area += ((p1.y * p2.x) - (p1.x*p2.y))/2;
            }
            area = Math.abs(area);
        }
    }
    return area;
}

function rad(degrees) {
  return degrees * PI / 180;
}

Where points stores values like: {lng: -73.462556156587, lat: 45.48566183708046}

Leafy

Answer 11

Google Maps Utils library provide a method to calculate the Area.

You have to add below dependency

implementation("com.google.maps.android:android-maps-utils:3.8.2")

And then call below method and give it yours Latitude and Longitude list. It will return area in square meters.

fun surfaceArea(list: List<LatLng>): Double {
    if (list.size < 3) {
        return 0.0
    }
    return SphericalUtil.computeArea(list)
}

Selfinterest