| |
 |
|
|
Science Forum Index » Geology - Satellite Navigation Forum » polygonal area on a sphere
Page 2 of 2 Goto page Previous 1, 2
|
| Author |
Message |
| Guest |
Posted: Tue Feb 27, 2007 4:39 am |
|
|
|
|
Quote: If you wish constructive answers, you need to provide more information:
Purpose of calculation
Accuracy required
"diameter" of polygon. i.e. what is the maximum of distance between
any two points on the polygon.
How is your math? Calculus? Trig?
If you intend to program the solution, what language?
I need to calculate the area in a rescue alpine region.
I suppose that the area will be few km2 and the maximum distance
between adjacent points will be some hundreds meters, therefore I can
suppose the region planar.
I know trig and I've read some articles on internet about this problem
(spherycal trigonometry, spherycal excess, etc.) but I don't find a
complete solution for my problem
For programming the language used is vb.net (or C#).
Thanks. |
|
|
| Back to top |
|
| J. J. Lodder |
Posted: Tue Feb 27, 2007 6:13 am |
|
|
|
Guest
|
Quote: Hi,
please, someone can tell me how can I to calculate the area of a
polygon on the earth's surface where I know the coordinates (lat/lon)
of its points?
The area is equal to the spherical excess,
which can be calculated directly from the coordinates.
For small polygons use plane trig,
Please, can you teach me or indicate an url where I can learn to
calculate the spherical excees from (lon, lat) coordinates?
Calculating an angle of a spherical triangle with given corners
is standard spherical trigonometry.
(Look up the formula for yourself,
there's no point in me retyping it here)
Just calculate them all and add them up.
Subtracht the plane sum of angles,
and you have the spherical excess.
By a standard theorem it equals the area
of your polygon.
Example: From the North pole walk to the equator,
then 90 degrees east, then back north to the pole.
This gives you a triangle with three straight angles,
hence (subtracting pi) with a sperical excess of pi/2.
As it should since by inspection the area is 4pi/8.
Note however that the subtraction is ill-conditioned
for small polygons.
A pocket calculator won't do.
Best,
Jan |
|
|
| Back to top |
|
| J. J. Lodder |
Posted: Tue Feb 27, 2007 6:13 am |
|
|
|
Guest
|
<luca.beltramo@gmail.com> wrote:
Quote: Hi,
please, someone can tell me how can I to calculate the area of a
polygon on the earth's surface where I know the coordinates (lat/lon)
of its points?
The area is equal to the spherical excess,
which can be calculated directly from the coordinates.
For small polygons use plane trig,
Please, can you teach me or indicate an url where I can learn to
calculate the spherical excees from (lon, lat) coordinates?
Thanks
Hi,
please, someone can tell me how can I to calculate the area of a
polygon on the earth's surface where I know the coordinates (lat/lon)
of its points?
The area is equal to the spherical excess,
which can be calculated directly from the coordinates.
For small polygons use plane trig,
Please, can you teach me or indicate an url where I can learn to
calculate the spherical excees from (lon, lat) coordinates?
Calculating an angle of a spherical triangle with given corners
is standard spherical trigonometry.
(Look up the formula for yourself,
there's no point in me retyping it here)
Just calculate them all and add them up.
Subtracht the plane sum of angles,
and you have the spherical excess.
By a standard theorem it equals the area
of your polygon.
Example: From the North pole walk to the equator,
then 90 degrees east, then back north to the pole.
This gives you a triangle with three straight angles,
hence (subtracting pi) with a sperical excess of pi/2.
As it should, since by inspection the area is 4pi/8.
Note however that the subtraction is ill-conditioned
for small polygons. (but there you can use plane geo)
A pocket calculator won't do.
Best,
Jan |
|
|
| Back to top |
|
| Ted Edwards |
Posted: Tue Feb 27, 2007 3:21 pm |
|
|
|
Guest
|
luca.beltramo@gmail.com wrote:
Quote: I need to calculate the area in a rescue alpine region.
I suppose that the area will be few km2 and the maximum distance
between adjacent points will be some hundreds meters, therefore I can
suppose the region planar.
I know trig and I've read some articles on internet about this problem
(spherycal trigonometry, spherycal excess, etc.) but I don't find a
complete solution for my problem
For programming the language used is vb.net (or C#).
Thanks.
Choose one of your points as origin and use the formulae from
http://williams.best.vwh.net/avform.html#flat to give you distance north
and east for each point. You can then use the formula I posted earlier
in this thread for the area of a polygon in the X-Y plane.
L¢1000 õ Let L = 1000 meters
LL¢45 100 õ Let LL be lat = 45, lon = 100 as reference position
«LL1¢45 100 SodanoD 1000 0 õ LL1 is a point 1000m north of LL
45.00899832 100
«LL2¢45 100 SodanoD 1000 90 õ LL1 is a point 1000m east of LL
44.9999993 100.0126828
LL Sodano LL1 õ distance on the WGS84 ellipsoid betwen LL and LL1
1000.000036 360
LL Sodano LL2 õ distance on the WGS84 ellipsoid betwen LL and LL2
1000.000036 90
If you can tolerate 0.036mm error in a length of 1000m, the flat earth
approximation should do the job for you. :-)
As to the bulge, look up the formula for the area in a figure bounded by
an arc and its chord. For a 1000m chord and a 6000Km radius (roughly
that of Earth), this is (in sq.m)
6E6 ArcArea1 1000
13.88888884
for a 500m chord, this is
6E6 ArcArea1 500
1.736111164
As you can see, this decreases as the cube of the chord length
÷/6E6 ArcArea1 1000 500
7.999999725
A "few km2" is ~3E6m2 so neglecting the above bulge would give an error
of only a few parts per million. Also, note that the area calculated
above is the area bounded by the arc of the Earth's surface and the
chord through the Earth. This is always greater that the area of the
projection of that figure onto the plane of interest so the error
estimate is an upper bound.
I would happily give you my code but it is in APL2. Do you have access
to APL2?
Ted |
|
|
| Back to top |
|
| Guest |
Posted: Wed Feb 28, 2007 11:43 am |
|
|
|
|
Quote: Choose one of your points as origin and use the formulae fromhttp://williams.best.vwh.net/avform.html#flatto give you distance north
and east for each point. You can then use the formula I posted earlier
in this thread for the area of a polygon in the X-Y plane.
L¢1000 õ Let L = 1000 meters
LL¢45 100 õ Let LL be lat = 45, lon = 100 as reference position
«LL1¢45 100 SodanoD 1000 0 õ LL1 is a point 1000m north of LL
45.00899832 100
«LL2¢45 100 SodanoD 1000 90 õ LL1 is a point 1000m east of LL
44.9999993 100.0126828
LL Sodano LL1 õ distance on the WGS84 ellipsoid betwen LL and LL1
1000.000036 360
LL Sodano LL2 õ distance on the WGS84 ellipsoid betwen LL and LL2
1000.000036 90
If you can tolerate 0.036mm error in a length of 1000m, the flat earth
approximation should do the job for you. :-)
As to the bulge, look up the formula for the area in a figure bounded by
an arc and its chord. For a 1000m chord and a 6000Km radius (roughly
that of Earth), this is (in sq.m)
6E6 ArcArea1 1000
13.88888884
for a 500m chord, this is
6E6 ArcArea1 500
1.736111164
As you can see, this decreases as the cube of the chord length
÷/6E6 ArcArea1 1000 500
7.999999725
A "few km2" is ~3E6m2 so neglecting the above bulge would give an error
of only a few parts per million. Also, note that the area calculated
above is the area bounded by the arc of the Earth's surface and the
chord through the Earth. This is always greater that the area of the
projection of that figure onto the plane of interest so the error
estimate is an upper bound.
I would happily give you my code but it is in APL2. Do you have access
to APL2?
Unfortunately I don't have access to APL2 and I don't know it.
However thank you very much for the explanation. It's clear and
precise!
.... but I have still a questio about an affermation that I don't
understand.
Why you say " A "few km2" is ~3E6m2"? From which consideration derives
"3E6m2"?
Thanks,
Luca |
|
|
| Back to top |
|
| Ted Edwards |
Posted: Wed Feb 28, 2007 1:47 pm |
|
|
|
Guest
|
luca.beltramo@gmail.com wrote:
Quote: ... but I have still a questio about an affermation that I don't
understand.
Why you say " A "few km2" is ~3E6m2"? From which consideration derives
"3E6m2"?
You stated, "I suppose that the area will be few km2 and the maximum
distance between adjacent points will be some hundreds meters,..."
"a few" is ~3 and I assume that by km2 you mean square kilometers. I
square kilometer is a figure 1km by 1km or any shape of the same area.
1km by 1km is 1000m by 1000m and has an area of 1000x1000=1E6 meters
squared. Thus you are interested in an area of the order of 3E6 square
meters. Since the potential errors are in meters or square meters, I
prefer to keep to the same units throughout.
Ted |
|
|
| Back to top |
|
| Dale DePriest |
Posted: Wed Feb 28, 2007 2:49 pm |
|
|
|
Guest
|
luca.beltramo@gmail.com wrote:
Quote: If your polygon is small enough to allow you to ignore the Earth's
curvature, you can use the following formula. You need to convert
lat/lon to x-y coordinates. If the polygon lies within one UTM region,
you can just convert lat/lon to UTM and use the UTM coordinates.
area = sum_{i=0}^{n-1} ((x_i + x_{i+1}) (y_{i+1} - y_i)) / 2
where point 0 = point N.
Is there a general formula in case of the area is between two UTM
regions?
No. If you need to work between two regions then use lat/lon degree
values instead of UTM values. If the point is near to another region you
can often use UTM values for one region by extending the values into the
second region. Even maps show these extensions with secondary markers
when the region switches in a populated area.
--
_ _ Dale DePriest
/`) _ // http://users.cwnet.com/dalede
o/_/ (_(_X_(` For GPS and GPS/PDAs |
|
|
| Back to top |
|
| ed williams |
Posted: Fri Mar 02, 2007 12:16 pm |
|
|
|
Guest
|
On Feb 28, 7:43 am, luca.beltr...@gmail.com wrote:
Quote: Choose one of your points as origin and use the formulae fromhttp://williams.best.vwh.net/avform.html#flattogive you distance north
and east for each point. You can then use the formula I posted earlier
in this thread for the area of a polygon in the X-Y plane.
L¢1000 õ Let L = 1000 meters
LL¢45 100 õ Let LL be lat = 45, lon = 100 as reference position
«LL1¢45 100 SodanoD 1000 0 õ LL1 is a point 1000m north of LL
45.00899832 100
«LL2¢45 100 SodanoD 1000 90 õ LL1 is a point 1000m east of LL
44.9999993 100.0126828
LL Sodano LL1 õ distance on the WGS84 ellipsoid betwen LL and LL1
1000.000036 360
LL Sodano LL2 õ distance on the WGS84 ellipsoid betwen LL and LL2
1000.000036 90
If you can tolerate 0.036mm error in a length of 1000m, the flat earth
approximation should do the job for you. :-)
As to the bulge, look up the formula for the area in a figure bounded by
an arc and its chord. For a 1000m chord and a 6000Km radius (roughly
that of Earth), this is (in sq.m)
6E6 ArcArea1 1000
13.88888884
for a 500m chord, this is
6E6 ArcArea1 500
1.736111164
As you can see, this decreases as the cube of the chord length
÷/6E6 ArcArea1 1000 500
7.999999725
A "few km2" is ~3E6m2 so neglecting the above bulge would give an error
of only a few parts per million. Also, note that the area calculated
above is the area bounded by the arc of the Earth's surface and the
chord through the Earth. This is always greater that the area of the
projection of that figure onto the plane of interest so the error
estimate is an upper bound.
I would happily give you my code but it is in APL2. Do you have access
to APL2?
Unfortunately I don't have access to APL2 and I don't know it.
However thank you very much for the explanation. It's clear and
precise!
... but I have still a questio about an affermation that I don't
understand.
Why you say " A "few km2" is ~3E6m2"? From which consideration derives
"3E6m2"?
Thanks,
Luca
(1) If the area is small enough that you can consider the earth as
flat:
Choose a point as origin.
Calculate the (x,y) (ie distance_East, distance_North)
coordinates of the vertices with respect to your origin.
http://williams.best.vwh.net/avform.htm#flat
Calculate the area of the plane polygon (by Method 1 in
http://www.efg2.com/Lab/Graphics/PolygonArea.htm )
(2) If the sphericity of the earth is important, and you have a
polygonal figure bounded by great circular arcs:
Triangulate the polygon.
Find the lengths of the sides of the triangles using
http://williams.best.vwh.net/avform.htm#Dist
Sum the areas of the individual triangles using l"Huiller's
formula http://williams.best.vwh.net/avform.htm#Triangle
Ed |
|
|
| Back to top |
|
| Ted Edwards |
Posted: Fri Mar 02, 2007 5:49 pm |
|
|
|
Guest
|
ed williams wrote:
The formula given there is
area = ( sum_{i=0...n-1} (x[i] y[i+1])-(x[i+1] y[i]) )/2
an alternate which, IIRC, is less subject to round off errors is
area = ( sum_{i=0...n-1} (x[i]+x[i+1]) (y[i+1] - y[i]) )/2
which is readily expressed in APL2 as
A{<-}0.5{times}+/(2+/x){times}{neg}2-/y
where x is the x-coordinates and y, the corresponding y-coordinates.
Ted |
|
|
| Back to top |
|
| ed williams |
Posted: Sat Mar 03, 2007 11:56 am |
|
|
|
Guest
|
On Mar 2, 1:49 pm, Ted Edwards <Ted_Espaml...@telus.net> wrote:
Quote: ed williams wrote:
Calculate the area of the plane polygon (by Method 1 in
http://www.efg2.com/Lab/Graphics/PolygonArea.htm)
The formula given there is
area = ( sum_{i=0...n-1} (x[i] y[i+1])-(x[i+1] y[i]) )/2
an alternate which, IIRC, is less subject to round off errors is
area = ( sum_{i=0...n-1} (x[i]+x[i+1]) (y[i+1] - y[i]) )/2
which is readily expressed in APL2 as
A{<-}0.5{times}+/(2+/x){times}{neg}2-/y
where x is the x-coordinates and y, the corresponding y-coordinates.
Ted
You can reduce the number of arithmetic operations and have rounding
benefits by using:
sum_{i=1..n} x(i) (y(i+1)-y(i-1))/2
where we have further extended the array so that (x(n+1),y(n
+1))=(x(1),y(1)) as well as (x(n),y(n))=(x(0),y(0))
See, for instance:
http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#2D%20Polygons |
|
|
| Back to top |
|
| Ted Edwards |
Posted: Sat Mar 03, 2007 2:38 pm |
|
|
|
Guest
|
ed williams wrote:
Quote: On Mar 2, 1:49 pm, Ted Edwards <Ted_Espaml...@telus.net> wrote:
area = ( sum_{i=0...n-1} (x[i] y[i+1])-(x[i+1] y[i]) )/2
area = ( sum_{i=0...n-1} (x[i]+x[i+1]) (y[i+1] - y[i]) )/2
You can reduce the number of arithmetic operations and have rounding
benefits by using:
sum_{i=1..n} x(i) (y(i+1)-y(i-1))/2
where we have further extended the array so that (x(n+1),y(n
+1))=(x(1),y(1)) as well as (x(n),y(n))=(x(0),y(0))
In what follows z is a 182 point GPS track created by combining the
tracks from several walks in the nearby hills. In order to have one
continuous track, I combined them such that there are some legs that go
back over themselves. The result is a very "messy" polygon.
A1 and A2 are the two formulae from my previous post (repeated above)
and A3 is the formula you just gave. Areas are in square meters. Eps
is machine epsilon. A couple of examples simple enough to work by hand
showed zero difference in the results.
z{<-}meadow
A1{<-}0.5{times}+/({times}/0 1{rot_}z)-{times}/1 0{rot_}z
A2{<-}0.5{times}+/(2+/z[;1]){times}{neg}2-/z[;2]
A3{<-}.5{times}+/z[;1]{times}{first}-/1 {neg}1{rot}{each}{enclose}z[;2]
(2-/A1 A2 A3 A1){div}A3 {rem} relative differences
1.914861534E-15 -5.106297424E-16 -1.404231792E-15
{rem} A2-A1 A3-A2 A1-A3
A3
455967.65
Eps
2.220446049E-16
It seems the closest agreement is between A3 and A2 and, as you point
out, A3 involves the least computation.
Ted |
|
|
| Back to top |
|
| ed williams |
Posted: Sun Mar 04, 2007 12:04 pm |
|
|
|
Guest
|
On Mar 3, 10:38 am, Ted Edwards <Ted_Espaml...@telus.net> wrote:
Quote: ed williams wrote:
On Mar 2, 1:49 pm, Ted Edwards <Ted_Espaml...@telus.net> wrote:
area = ( sum_{i=0...n-1} (x[i] y[i+1])-(x[i+1] y[i]) )/2
area = ( sum_{i=0...n-1} (x[i]+x[i+1]) (y[i+1] - y[i]) )/2
You can reduce the number of arithmetic operations and have rounding
benefits by using:
sum_{i=1..n} x(i) (y(i+1)-y(i-1))/2
where we have further extended the array so that (x(n+1),y(n
+1))=(x(1),y(1)) as well as (x(n),y(n))=(x(0),y(0))
In what follows z is a 182 point GPS track created by combining the
tracks from several walks in the nearby hills. In order to have one
continuous track, I combined them such that there are some legs that go
back over themselves. The result is a very "messy" polygon.
A1 and A2 are the two formulae from my previous post (repeated above)
and A3 is the formula you just gave. Areas are in square meters. Eps
is machine epsilon. A couple of examples simple enough to work by hand
showed zero difference in the results.
z{<-}meadow
A1{<-}0.5{times}+/({times}/0 1{rot_}z)-{times}/1 0{rot_}z
A2{<-}0.5{times}+/(2+/z[;1]){times}{neg}2-/z[;2]
A3{<-}.5{times}+/z[;1]{times}{first}-/1 {neg}1{rot}{each}{enclose}z[;2]
(2-/A1 A2 A3 A1){div}A3 {rem} relative differences
1.914861534E-15 -5.106297424E-16 -1.404231792E-15
{rem} A2-A1 A3-A2 A1-A3
A3
455967.65
Eps
2.220446049E-16
It seems the closest agreement is between A3 and A2 and, as you point
out, A3 involves the least computation.
Ted
For lurkers, the area computed by this formula is signed. The area
is positive if its interior lies to the left of the point sequence and
vice-versa. With a figure-8 loop, you'd get the difference of the area
of the two sub-loops. |
|
|
| Back to top |
|
| Ted Edwards |
Posted: Sun Mar 04, 2007 1:48 pm |
|
|
|
Guest
|
ed williams wrote:
Quote: For lurkers, the area computed by this formula is signed. The area
is positive if its interior lies to the left of the point sequence
and vice-versa. With a figure-8 loop, you'd get the difference of the
area of the two sub-loops.
Also very useful if you have something with holes. e.g. a sheet metal
layout for which you wish to calculate weight. Get area by going left
around the outline and right around the holes. Also useful for area of
a lake excluding islands.
Ted |
|
|
| Back to top |
|
| ed williams |
Posted: Sun Mar 04, 2007 11:06 pm |
|
|
|
Guest
|
On Mar 4, 9:48 am, Ted Edwards <Ted_Espaml...@telus.net> wrote:
Quote: ed williams wrote:
For lurkers, the area computed by this formula is signed. The area
is positive if its interior lies to the left of the point sequence
and vice-versa. With a figure-8 loop, you'd get the difference of the
area of the two sub-loops.
Also very useful if you have something with holes. e.g. a sheet metal
layout for which you wish to calculate weight. Get area by going left
around the outline and right around the holes. Also useful for area of
a lake excluding islands.
Ted
There is an exact, but very sensitive to rounding error, formula for
the area of an n-sided spherical polygon, generalizing the more well-
known relation to the spherical excess:
Area = (Sum (Interior_angles - pi) + 2*pi) * R^2
This will be horrible for small polygons, good for large ones
(compared to the surface area of the earth) (because the formula gives
zero in the limit of a plane polygon). ll'Huiller's formula, that I
recommended earlier, does not have this problem.
Ed
http://historical.library.cornell.edu/cgi-bin/cul.math/docviewer?did=00640001&seq=79&frames=0&view=50 |
|
|
| Back to top |
|
| |
Page 2 of 2 Goto page Previous 1, 2
All times are GMT - 5 Hours
The time now is Tue Oct 07, 2008 8:35 pm
|
|