How do you calculate area of polygon where polygon is defined as array lon/lat points?

[Scott Howlett]

Hi everyone -

I'm certain to be embarrassed at the simplicity of the solution, but I can't for the life of me figure out how to construct a polygon from an array of lon/lat points and then get a meaningful area (i.e. in square meters, kilometers, miles or feet).  I'm guessing I need some type of transformation of coordinate systems?  I've spent a few hours on this but haven't got anywhere (I'm definitely not a GIS expert!) -> I thought this was going to be simple.

Anyway, here's the code:

        public static void TestPolygonArea()

        {

            string coordinates = "-79.525542519049552,43.691278124243432 -79.525382520578987,43.691281097414787 -79.525228855617627,43.69124858593392 -79.525096151437353,43.691183664769774 -79.52472799258571,43.690927163079735 -79.525379447437814,43.690771996666641 -79.525602330675355,43.691267524226838 -79.525542519049552,43.691278124243432";

            List<GeoAPI.Geometries.Coordinate> points = new List<GeoAPI.Geometries.Coordinate>();

            foreach (string coordinate in coordinates.Split(' '))

            {

                double lon = double.Parse(coordinate.Split(',')[0]);

                double lat = double.Parse(coordinate.Split(',')[1]);

                points.Add(new GeoAPI.Geometries.Coordinate(lon, lat));

            }

            NetTopologySuite.Geometries.LinearRing ring = new NetTopologySuite.Geometries.LinearRing(points.ToArray());

            NetTopologySuite.Geometries.Polygon polygon = new NetTopologySuite.Geometries.Polygon(ring);

            

            // this is always 0

            double ringArea = ring.Area;

            // this is 0.00000024442176540574213.  not sure what unit this is in.

            double polygonArea = polygon.Area;

        }

The polygon in question is in Toronto, Canada and I would expect the area to be somewhere in the range of 500 to 1500 square meters.  It plots super easy in Google maps, but I'd also like to be able to show the area.  Any help would be greatly appreciated and apologies in advance for such an simple question - as mentioned though, I couldn't find out how to do this.

[John Diss]

Hi Scott, your current polygon is described using a geographic coordinate system i.e Lon Lat. To get a useful area you need to transform it into a planar coordinate system and then the value of the polygon.Area property will depend on the unit of the chosen coordinate system. Depending on the coordinate system you transform to, you will get variations in the area returned.

 

HTH jd

--
You received this message because you are subscribed to the Google Groups "NetTopologySuite" group.
To view this discussion on the web visit https://groups.google.com/d/msg/nettopologysuite/-/6bQjSqBxm-MJ.
To post to this group, send email to nettopol...@googlegroups.com.
To unsubscribe from this group, send email to nettopologysui...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/nettopologysuite?hl=en.

[Scott Howlett]

Thanks John.  I had suspected that I needed some type of coordinate transformation.

I'm obviously quite new to this - any chance you might provide:

• a recommendation on what planar coordinate system to use?  As mentioned the shapes are primarily in Toronto, Canada but I expect to have them in other areas of Canada as well.  I'm not overly concerned with accuracy or performance as the calc will be done offline.
• some sample code as to how to do the actual transformation?

Given the extreme popularity of the lon/lat coordinate system, I'm surprised that the operation to do this isn't supported right in NTS:

something like:

string coordinates = "-79.525542519049552,43.691278124243432 -79.525382520578987,43.691281097414787 -79.525228855617627,43.69124858593392 -79.525096151437353,43.691183664769774 -79.52472799258571,43.690927163079735 -79.525379447437814,43.690771996666641 -79.525602330675355,43.691267524226838 -79.525542519049552,43.691278124243432";

var polygon = new LonLatPolygon(coordinates);

double areaInKm = polygon.AreaInKm;

I guess that's hoping for too much!

Thanks again.

/s.

To post to this group, send email to nettopo...@googlegroups.com.
To unsubscribe from this group, send email to nettopologysu...@googlegroups.com.

[prusinek74]

Hi everyone.

The truth is there it  is no need to transform these coordinates to planar coordinate system, to calculate polygon area. 

Here i found some theory on spherical geometry and there is also some words about calculating area. But it is not as simply as multiplying and i think some performance penalties are possible too.

Regards

Piotr

[prusinek74]

Forgot to link document

Here it is: http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/40409/1/07-03.pdf

[FObermaier]

Thanks for the input

[FObermaier]

if the data you are using is more or less static, I'd go about and transform those geometries to a planar coordinate system and not on the fly.

You can easily do that using e.g. proj4 directly or other utilities that use it like OGR, PostGis, SpatiaLite.

The coordinate system to use depends on the area of interst. For the toronto area, maybe UTM ZONE 17N is a good choice? You can find the definition for e.g. Proj4, Well known text on www.spatialreference.org. You MUST know what Geographic coordinate system your coordinates are in (WGS72, WGS84?) to get the best result.

If you want/need to do that on the fly, checkout the GeometryTransformTest method of the GoogleGroupsTest class
(http://code.google.com/p/nettopologysuite/source/browse/trunk/NetTopologySuite.Samples.Console/Tests/Various/GoogleGroupTests.cs#386)

Hth FObermaier

[FObermaier]

Forward from Ted Macy:

If all you want to do is determine true area, transform to a local Transverse Mercator projection, using the centroid of the geometry as the origin of the projection.    We run into this frequently... discrepancies in area for a field located near the border of a UTM zone.   In agricultural standards efforts, the trend is to always use local transforms for measuring distance, area and parallel lines.  When working with a piece of equipment in a field, the equipment is a fixed width, and will produce the same number of paths for a given field, but if you don’t you a local projection, you will mathematically get a different number of paths, depending upon the projection you use.

______________________

Ted Macy

MapShots, Inc.

Web:     www.mapshots.com

On Monday, January 7, 2013 9:16:53 PM UTC+1, Scott Howlett wrote:

[Scott Howlett]

For what it's worth, this is what I ended up doing.  It does a transform and gives a result that looks about right.  Would appreciate confirmation from you knowledgeable folks if this looks about right.  I ended up using DotSpatial - there seems to be a lot of overlap in SharpMap, NetTopologySuite, DotSpatial, Proj4 and I'm not certain which is best for my use.  Again, I was thinking this would be easy so I'm just getting something that works (or that seems to work).

---------------

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using DotSpatial.Projections;

using DotSpatial.Data;

using DotSpatial.Topology;

namespace Test

{

    class Program

    {

        static void Main(string[] args)

        {

            TestPolygonArea();

        }

        public static void TestPolygonArea()

        {

            // this feature can be see visually here http://www.allhx.ca/on/toronto/westmount-park-road/25/

            string feature = "-79.525542519049552,43.691278124243432 -79.525382520578987,43.691281097414787 -79.525228855617627,43.69124858593392 -79.525096151437353,43.691183664769774 -79.52472799258571,43.690927163079735 -79.525379447437814,43.690771996666641 -79.525602330675355,43.691267524226838 -79.525542519049552,43.691278124243432";

            string[] coordinates = feature.Split(' ');

            // dotspatial takes the x,y in a single array, and z in a separate array.  I'm sure there's a 

            // reason for this, but I don't know what it is.

            double[] xy = new double[coordinates.Length*2];

            double[] z = new double[coordinates.Length];

            for (int i = 0; i < coordinates.Length; i++)

            {

                double lon = double.Parse(coordinates[i].Split(',')[0]);

                double lat = double.Parse(coordinates[i].Split(',')[1]);

                xy[i*2] = lon;

                xy[i * 2 + 1] = lat;

                z[i] = 0;

            }

            ProjectionInfo pStart = KnownCoordinateSystems.Geographic.World.WGS1984;

            //which UTM zone to use can be found here http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system

            //the feature is in Toronto, Canada which is Zone17.  Technically, should use 17T I think but this isn't defined

            //so am just using 17N - not sure what impact this has.

            ProjectionInfo pEnd = KnownCoordinateSystems.Projected.UtmNad1927.NAD1927UTMZone17N;

            // do the actual reprojection

            Reproject.ReprojectPoints(xy, z, pStart, pEnd, 0, coordinates.Length);

            // build up a list of Coordinate, to create the polygon

            List<Coordinate> co = new List<Coordinate>();

            for (int i = 0; i < coordinates.Length; i++)

            {

                co.Add(new Coordinate(xy[i*2],xy[i*2+1]));

            }

            Polygon polygon = new Polygon(co);

            double area = polygon.Area;

            // area is 2188 sq m