|
|
| GeoCommunity Mailing List |
| |
| Mailing List Archives |
| Subject: | Re: [gislist] Calculate Area |
| Date: |
04/04/2005 08:15:01 PM |
| From: |
Bill Thoen |
|
|
On Mon, 4 Apr 2005, Douglas F. Wunneburger wrote:
> Click on the Help menu button, then ArcGIS desktop Help ...
>
> (Regarding recent complaints about using this list for ESRI support, I
> consider this to be a general question as the same answer should apply to
> your software package of choice.)
I just can't resist. The following response should point out out how far
we've come from actually talking about GIS algorithms to just looking for
help off a menu. When you find the menu option to calculate area, will
you know (or care) how it's actually done? The algorithm used really does
matter.
Here's how this question used to be answered when the GIS forum was called
GIS-L (and it's a really great answer):
---
Early in July 1997, John Chioles wrote:
> Can anyone offer the calculation/source code for calculating
> the area of a [spherical] polygon?
Bill Huber replies:
There are some relatively simple ways, depending on your needs for
accuracy, the kinds of polygons you are considering, and how they are
represented.
If:
* Areas to about four sig. figs are good enough
* Polygons do not sprawl all over the globe
* They are represented in "spaghetti code" that provides the coordinates
of vertices in decimal degrees then there are some nice algorithms.
This posting provides an algorithm expressed in working, tested, Avenue
code (for ArcView 3.0 and later). The code is relatively simple and short
(less than 100 executable lines, plus many comments), and so is readily
amenable to transcription into most programming languages (including those
independent of any GIS platform). It has to use double precision floating
point arithmetic or better.
Description
--------------
The algorithm assumes that a polygon is represented as a collection of
connected components. Each component is represented as a sequence of
vertices, traversed so that the interior is kept to the left
(mathematically "positive" orientation) for positive areas and to the
right for negative areas (representing "holes," for example). The
components must "wrap around" by including the same point both first and
last in the list. All component areas are assumed to be less than one
half the area of the earth. The line segments between vertices are
assumed to follow geodesics along the GRS 80 ellipsoid and the vertices
are assumed to represent the termini of these geodesics in decimal degrees
as (longitude, latitude): longitudes increase to the east and can be any
real value: the latitudes must be between -90 and 90.
The output is a double-precision floating point number providing the
signed area of the polygon, in square miles. The asymptotic storage and
computational demands of the algorithm are directly proportional to the
number of vertices.
The algorithm is presented in Avenue code with extensive comments, copied
directly from scripts that have been compiled and run successfully. It
has been coded for clarity and "transcribability," rather than speed and
efficiency of resources.
Results
---------
Extensive testing with world countries, U.S. states, and U.S. zip code-3
boundaries provided interesting and illuminating results. Some are:
1. The alternative of computing areas in projected coordinates, even
using an equal-area projection, is not as accurate. The algorithm given
here is at least 20 times more precise for the polygons that were tested.
2. The accuracy of the presented algorithm is probably about one part in
10^6 or 10^7. Beyond that, there may be floating point roundoff errors.
3. The accuracy of vector data sets that meet usual standards will
provide between one part in 10^5 and one part in 10^3 accuracy, depending
on how good the data are and how many vertices the polygon has.
Therefore the algorithm presented should be more than adequate for
practical applications.
4. The presented algorithm is theoretically optimal in its use of
computing resources, for detailed polygons. With a little work its
running time could probably be halved, but further improvements are
unlikely.
Limitations
--------------
There are limitations to this algorithm. Basically, if
* Your polygon has millions of vertices, or
* It has thin arms that snake into outlying latitudes, or
* It encircles a pole (without explicitly visiting the pole), or
* Comes very close to both the north and the south poles, or
* Is extremely small (less than a 100 square miles-but I don't know how
much less),
then it could give results that are off by as much as 0.3% or it could
fail entirely. Fortunately, these are conditions that (a) are rare (b)
can be checked in code and (c) can be corrected by slicing the polygon
into smaller pieces and summing the areas of the pieces or by using a
plane-geometry algorithm for really small pieces, or both.
Th
|
|
Sponsored by:
For information
regarding
advertising rates
Click Here!
|
