Computing a polygon region using Latitude and longitude generated from Cartesian space and the world file

Question

Given the sequence of GPS coordinate pairs, I need to calculate the area of ​​the polygon (n-gon). This is relatively small (not more than 50,000 square feet). Geocodes are created by applying affine transformations with data from a world file.

I tried using a two-step approach by converting geocodes to Cartesian coordinates:

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

then I use cross product calculation to determine the area.

The problem is that the results are slightly inaccurate (about 1%). Is there anything I can learn to improve this?

Thank you

Answer 1

A 1% error seems a little high due to your approximation. Do you compare actual measurements or some perfect calculation? Remember that there is also a bug in GPS that may contribute.

If you need a more accurate method for this, there is a good answer to this question. If you intend to use a faster method, you can use the WGS84 geoid instead of your reference sphere to convert to Cartesian coordinates (ECEF). Here is a link to the wiki for this conversion.

Answer 2

I am modifying the Google map so that the user can calculate the area of ​​the polygon by clicking the vertices. This did not give the correct ones until I was convinced that Math.cos (latAnchor) was the first in radians.

So:

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

become:

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

where lon, lonAnchor and latAnchor are in degrees. Now it works like a charm.

Answer 3

I checked various formulas of the polygon area (or code) on the Internet, but did not find anything good or easy to implement.

Now I have written a code snippet to calculate the area of ​​a polygon drawn on the surface of the earth. A polygon can have n vertices, each of which has its latitudinal longitude.

Some important points

• The input of the array to this function will contain the elements "n + 1". The last element will have the same values ​​as the first.
• I wrote very simple code in C #, so the guys can also adapt it in another language.
• 6378137 - value of the radius of the earth in meters.

• In the output area there will be a unit of square meters

   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; } 

Answer 4

Based on the Pathhaka Risk solution, here is a SQL solution (Redshift) for calculating areas for GeoJSON multipolygons (assuming that line row 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 ; 

Answer 5

Adapted RiskyPathak Snippet for 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; } 

Answer 6

Thanks Risks Pathak!

In the spirit of sharing, here is 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;