View Full Version : Surface Area

Dazzer

12-06-2005, 08:09 AM

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 :wink: ) .

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).

Wilfried

12-06-2005, 09:21 AM

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.

Dazzer

12-06-2005, 09:56 AM

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

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

DrExcitement

12-06-2005, 01:21 PM

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/~pbourke/geometry/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.

Wilfried

12-06-2005, 02:03 PM

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/143/surfacearea3kt.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 ?

Dazzer

12-08-2005, 03:54 AM

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?

Wilfried

12-08-2005, 07:37 AM

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).

DrExcitement

12-08-2005, 09:50 AM

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.

Dazzer

12-08-2005, 05:27 PM

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.

Dazzer

01-05-2006, 04:13 AM

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.

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 :wink: !)

Powered by vBulletin® Version 4.2.2 Copyright © 2014 vBulletin Solutions, Inc. All rights reserved.

Search Engine Friendly URLs by

vBSEO 3.6.1