PDA

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