Heres an interesting one that I'm having a lot of difficulty with (i'm sure you maths boffins out there will scoff at my strugling with such a simple problem ) .

Given a set of lat\lon co-ordinates for a closed polygon shape how can I work out the surface area of that shape in sq miles or sq km (or feet or meters or something sensible!).

I've searched the internet but the solutions I find are all looking for the ultimate in accuracy and far to complex for me to understand! I'm really just looking at this from a 'Flat Earth' scenario it doesn't need to be extremly accurate (unless its just as easy to make it more accurate of course).

Hi,

No so simple but not really difficult also I try to explain. Take a pencil and draw a irregular polygon.

First step is, you take first point and draw a line to a second point skipping 1 point. What you now have is a triangle. Easy to calculated surface area of a triangle.

go on with a line to next point (skipping also just 1 point), do over and over again until you have the last resulting triangle.

I think I understand what your saying and I may be missing something but for the map shown here http://img384.imageshack.us/img384/1...acearea3kt.jpg

would the triangle shown with the red line not be outside the surface area of the shape?

I've wondered about this, too, but I've never found an answer that talks about polygons with (lat,long) vertices.

Here's how to do it with (x,y) vertices: http://astronomy.swin.edu.au/~pbourk...etry/polyarea/. This works for concave polygons like the one you showed. To use with (lat,long) vertices, I guess you could convert each vertex to an (x,y) in meters (or whatever), based on a local coordinate system, and assuming the world is locally flat. Clearly there are some approximation errors involved here.

Hi,

I think I understand what your saying and I may be missing something but for the map shown here http://img384.imageshack.us/img384/1...acearea3kt.jpg
would the triangle shown with the red line not be outside the surface area of the shape?

Yes that's what I meant in 'no so simple'. There are some complications here and there, but that's what we programmers do like eh dont we ?

Thanks for you help Wilfred, DrExcitement, it appears that this isn't as simple as I hoped (I didn't really think it would be to be honest!)

Looking at DrExcitements link it appears that this could be something I could use in my application if I had a way of converting Lat, Lon to some local Grid system, does anyone know of a way to do this?

Hi,

Just use lat/long pairs. eg multiply them with 10000 or so and make them integers. Using 1/10000 will give you accuratie around 11 meter. 1 minute is exact 1 sea mile.

You will have (small) calculation errors because 1 minute in longitude is 1 seamile on the equator and less of course if go more distance away for it. You can correct it by multiply it by Cos(Latitude).

A slight refinement of the reply above is to first compute the mean (lat,long) of all your vertices. This will be the origin of your local coordinate system, called (lat0,long0). Now for each (lat,long) vertext, the corresponding (x,y) is y = pi*earth_radius_meters*(lat - lat0)/180 and x = pi*earth_radius_meters*cos(lat0)*(long - long0)/180. This assumes (lat,long) is in degrees. earth_radius_meters is approximately 1000.0 * 6371.01. This will give you coordinates in meters, so your area will be in square meters. If you want another square area unit, use the corresponding linear unit for the radius of the earth.

Excellent, thank you both for your help, I should have the basis for a working algorithm now. I'll post a reply containing the class for any one else whos interested in this when i'm finished.

Not really a class for doing this but heres how I implemented this in my code (VB .NET), should give someone some ideas towards writing a class of there own. In this I have my coordinates held in the table 'Shapes' each shape has a unique ID.

Code:

    Imports System.Math

    Private Function CalculateArea(ByVal ID As String) As String

        Dim con As New System.Data.SqlClient.SqlConnection(ConnString)
        con.Open()
        Dim cmd As New System.Data.SqlClient.SqlCommand
        Dim da As New System.Data.SqlClient.SqlDataAdapter
        Dim ds As New DataSet
        Dim dr As DataRow
        Dim MyItem As ListViewItem
        Dim xCoord() As Double
        Dim yCoord() As Double
        Dim n As Integer
        Dim x, y, dist As Double

        cmd = con.CreateCommand
        cmd.CommandText = "SELECT * FROM Shapes WHERE ID = '" & ID & "' ORDER BY point"
        cmd.CommandType = CommandType.Text
        da.SelectCommand = cmd
        da.Fill(ds, "Coords")

        ReDim xCoord(ds.Tables("Coords").Rows.Count)
        ReDim yCoord(ds.Tables("Coords").Rows.Count)

        For n = 0 To UBound(xCoord)
            If n = UBound(xCoord) Then
                xCoord(n) = xCoord(0)
            Else
                xCoord(n) = ds.Tables("Coords").Rows(n).Item("Lat")
            End If

            If n = UBound(yCoord) Then
                yCoord(n) = yCoord(0)
            Else
                yCoord(n) = ds.Tables("Coords").Rows(n).Item("Lon")
            End If

        Next

        For n = 0 To UBound(xCoord)
            x += xCoord(n)
        Next

        x = x / (UBound(xCoord) + 1)

        dist = DistanceFrom(x, x, 0, 1)

        For n = 0 To UBound(xCoord)
            xCoord(n) = xCoord(n) * 69.11
        Next

        For n = 0 To UBound(yCoord)
            yCoord(n) = yCoord(n) * dist
        Next
        CalculateArea = Format(AreaByCoordinates(xCoord, yCoord), "#00.000 sqMiles")

    End Function
    Function AreaByCoordinates(ByVal Xcoord() As Double, ByVal Ycoord() As Double) As Double
        Dim I As Long
        Dim Xold As Double
        Dim Yold As Double
        Dim Yorig As Double
        Dim ArrayUpBound As Long
        Dim x, y As Double

        ArrayUpBound = UBound(Xcoord)
        Xold = Xcoord(ArrayUpBound)
        Yorig = Ycoord(ArrayUpBound)
        Yold = 0.0#

        For I = LBound(Xcoord) To ArrayUpBound
            x = Xcoord(I)
            y = Ycoord(I) - Yorig
            AreaByCoordinates = AreaByCoordinates + (Xold - x) * (Yold + y)
            Xold = x
            Yold = y
        Next

        AreaByCoordinates = Abs(AreaByCoordinates) / 2
    End Function
    Function DistanceFrom(ByVal lat1 As Double, ByVal lat2 As Double, ByVal lon1 As Double, ByVal lon2 As Double) As Double
        Dim theta, dist As Double
        theta = lon1 - lon2
        dist = Sin(deg2rad(lat1)) * Sin(deg2rad(lat2)) + Cos(deg2rad(lat1)) * Cos(deg2rad(lat2)) * Cos(deg2rad(theta))
        dist = Acos(dist)
        dist = rad2deg(dist)
        DistanceFrom = dist * 60 * 1.1515

    End Function
    Function deg2rad(ByVal deg As Double) As Double
        deg2rad = CDbl(deg * PI / 180)
    End Function
    Function rad2deg(ByVal rad As Double) As Double
        rad2deg = CDbl(rad * 180 / PI)
    End Function

I got the figure 69.11 from using this site

http://www.zodiacal.com/tools/lat_table.htm

I don't really have time to explain all of the details write now and this isn't entirely accurate due to the complications involved in making it so, but hopefully this should help anyone looking to do anything similar and perhaps they can expand on it and improve it and post there results back here. (Also if i've made any obvious mistakes could someone let me know !)