Prev: RPN/RPL Calculator implementations, list of, regular post [long, FAQ]
Next: Graphics programming commands?
From: PremiumBlend on 6 May 2010 20:57
Hello,

I rewrote a program named AREA that I posted
previously on Sept. 9, 2009:

http://groups.google.com/group/comp.sys.hp48/browse_thread/thread/afc02915d1db1129/8d2e3bc581

Using suggestions by Jacob Wall and a code snippet
written by Virgil, I came up with a more concise
program that uses x, y values input via the INPUT
command instead of recalling them from a file. Here
it is:

\<< 28 MENU
"DEFINE AREA
BY COORDINATES:"
{ "" { -1 0 } v }
IFERR INPUT
THEN DROP2 0 MENU
ELSE OBJ\-> DEPTH
\->LIST DUP SIZE DUP
2 /
\<<
FOR x <-\list x
GET 2
STEP <-\k \->LIST
DUP DUP HEAD + TAIL
\>> \-> <-\list j <-\k
Filter
\<< 1 j 1 -
Filter EVAL 2 j
Filter EVAL 4 ROLL
* \GSLIST SWAP ROT *
\GSLIST - ABS 2 /
\>>
END
\>>

The program uses the cross coordinate method
which is employed by land surveyors and is
also known as the shoelace formula:

http://en.wikipedia.org/wiki/Shoelace_formula

A land surveyor might use this method to
calculate the area of a tract of land using
the following data:

State Plane Coordinates (m)

Northing (Y) Easting (X)

60,320.528 1,395,020.737
60,560.076 1,399,870.961
56,300.124 1,397,560.850
57,867.993 1,401,028.295
58,670.873 1,397,840.621

Area = 2589889 sq. meters

Elementary Surveying Ninth
Edition by Wolf and Brinker
page 274 states:

"A word of caution applies, however,
if coordinate values are extremely
large, as they would normally be-for
example, if state plane coordinates
are used (see Chapter 21). In those
cases, to ensure sufficient precision
and prevent serious round-off errors,
double precision should be used. Or,
as an alternative, the decimal place
in each coordinate can arbitrarily be
moved n places to the left, the area
calculated, and then multiplied by
10^2n."

Using the alternative suggestion, I ran
the program again but this time around
I moved all the decimal points three places
to the left. I then multiplied the result by
10^6. I got the same result as before.

Could someone explain to me why I
didn't get a different result?

Mark
From: Jacob Wall on 6 May 2010 22:07
> State Plane Coordinates (m)
>
> Northing (Y) Easting (X)
>
> 60,320.528 1,395,020.737
> 60,560.076 1,399,870.961
> 56,300.124 1,397,560.850
> 57,867.993 1,401,028.295
> 58,670.873 1,397,840.621
>
> Area = 2589889 sq. meters

Hello Mark,

Nothing to do with your question, but I looked at the coordinates and
the area you supplied and would like to point out that the sides of the
polygon you posted criss-cross, if joining the coordinates in the order
you posted. This is a problem that I have spent some time on, trying to
figure out a "intersection detection algorithm" but never did come up
with anything. For example, say you have 4 points,
P N E
1 0 0
2 10 0
3 10 10
4 0 10

Area = 100 units�

Now if you mix up the order, say:
P N E
1 0 0
3 10 10
2 10 0
4 0 10

Area = 0, with the methods used. Creating a 5th point at the center 5,5
and then calculating in the order 1 2 5 3 4 5 will give an answer of 50
units�, which is correct. The trick would be to come up with a way to
get the right answer when calculating in the order 1 3 2 4, although the
question quickly can become, "which are are you interested in" in these
cases.

The example I gave is pretty obvious of the limitations, but I've
encountered scenarios where it is not immediately obvious that an
incorrect area may be calculated if a polygon is not broken up into
smaller pieces, I don't have an example in front of me but they usually
involve multiple curves where segment areas are added or subtracted
after computing the polygon area. Even though a figure as drawn does
not criss-cross, when joining BC-EC and actually drawing the polygon it
becomes clear and generally the error becomes somewhat obvious when
looking at the answer and comparing to what would seem a logical estimate.

Anyways, just an aside.

Jacob
From: PremiumBlend on 6 May 2010 23:26
On May 6, 10:07 pm, Jacob Wall <jac...(a)surv50.ca> wrote:
> > State Plane Coordinates (m)
>
> > Northing (Y)  Easting (X)
>
> > 60,320.528   1,395,020.737
> > 60,560.076   1,399,870.961
> > 56,300.124   1,397,560.850
> > 57,867.993   1,401,028.295
> > 58,670.873   1,397,840.621
>
> > Area = 2589889 sq. meters
>
> Hello Mark,
>
> Nothing to do with your question, but I looked at the coordinates and
> the area you supplied and would like to point out that the sides of the
> polygon you posted criss-cross, if joining the coordinates in the order
> you posted.  This is a problem that I have spent some time on, trying to
> figure out a "intersection detection algorithm" but never did come up
> with anything.  For example, say you have 4 points,
> P  N   E
> 1  0   0
> 2  10  0
> 3  10  10
> 4  0   10
>
> Area = 100 units²
>
> Now if you mix up the order, say:
> P  N   E
> 1  0   0
> 3  10  10
> 2  10  0
> 4  0   10
>
> Area = 0, with the methods used.  Creating a 5th point at the center 5,5
> and then calculating in the order 1 2 5 3 4 5 will give an answer of 50
> units², which is correct.  The trick would be to come up with a way to
> get the right answer when calculating in the order 1 3 2 4, although the
> question quickly can become, "which are are you interested in" in these
> cases.
>
> The example I gave is pretty obvious of the limitations, but I've
> encountered scenarios where it is not immediately obvious that an
> incorrect area may be calculated if a polygon is not broken up into
> smaller pieces, I don't have an example in front of me but they usually
> involve multiple curves where segment areas are added or subtracted
> after computing the polygon area.  Even though a figure as drawn does
> not criss-cross, when joining BC-EC and actually drawing the polygon it
> becomes clear and generally the error becomes somewhat obvious when
> looking at the answer and comparing to what would seem a logical estimate..
>
> Anyways, just an aside.
>
> Jacob

Oops! I grabbed those coordinates from a book called Geodesy
for Geomatics and GIS Professionals by Elithorp and
Findorff, page 154. I didn't realize they didn't make a polyon by
the sequence in which they were listed. I was just trying to find
a table of state plane coordinates in order to show what I was
talking about. Do you have any tables of very large state plane
coordinates that I could post? Will get back to you on the
criss-cross issue.
From: Jacob Wall on 7 May 2010 00:41
> Do you have any tables of very large state plane
> coordinates that I could post? .

Could just swap 3 and 4:

State Plane Coordinates (m)

Northing (Y) Easting (X)

1 60,320.528 1,395,020.737
2 60,560.076 1,399,870.961
4 57,867.993 1,401,028.295
3 56,300.124 1,397,560.850
5 58,670.873 1,397,840.621

I get an area of 12,055,386.998 m�

With the approach I posted back in Sept. 2009:

\<< Nlist DUP HEAD - DUP TAIL SWAP HEAD + Elist DUP HEAD - * \GSLIST
Elist DUP HEAD - DUP TAIL SWAP HEAD + Nlist DUP HEAD - * \GSLIST - ABS
2. / \>>

the entire group of coordinates is "shifted" down to near the coordinate
system origin for the purpose of the calculation. This essentially does
the same thing as moving the decimal place over.

By calculating the same area without the "shift", I get an answer of:
12,055,387.000 m�, which really is the same thing. The program below is
without the "shift":

\<< Nlist DUP TAIL SWAP HEAD + Elist * \GSLIST Elist DUP TAIL SWAP HEAD
+ Nlist * \GSLIST - ABS 2. / \>>

So the reason for this is that there are only 5 points within the
polygon, if we add a bunch more, say add another point half way between
each of the points already listed, maintaining the shape and area:

1 60,320.5280 1,395,020.7370
6 60,440.3020 1,397,445.8490
2 60,560.0760 1,399,870.9610
7 59,214.0345 1,400,449.6280
4 57,867.9930 1,401,028.2950
8 57,084.0585 1,399,294.5725
3 56,300.1240 1,397,560.8500
9 57,485.4985 1,397,700.7355
5 58,670.8730 1,397,840.6210
10 59,495.7005 1,396,430.6790

Using the approach posted in 2009, I still get the same answer:
12,055,386.998 m�

Now without the "shift", the answer is:
12,055,388.000 m�

And it only gets worse from there on, the more points you add to the mix.

You can see how the longer the list of high numbered coordinate values
is, the larger a number \GSLIST will return, and with only 12 available
digits you run out of digits pretty quickly.

Jacob
From: PremiumBlend on 7 May 2010 23:34
On May 7, 12:41 am, Jacob Wall <jac...(a)surv50.ca> wrote:
> > Do you have any tables of very large state plane
> > coordinates that I could post? .
>
> Could just swap 3 and 4:
>
> State Plane Coordinates (m)
>
> Northing (Y)  Easting (X)
>
> 1 60,320.528   1,395,020.737
> 2 60,560.076   1,399,870.961
> 4 57,867.993   1,401,028.295
> 3 56,300.124   1,397,560.850
> 5 58,670.873   1,397,840.621
>
> I get an area of 12,055,386.998 m²
>
> With the approach I posted back in Sept. 2009:
>
> \<< Nlist DUP HEAD -  DUP TAIL SWAP HEAD + Elist DUP HEAD - * \GSLIST
> Elist DUP HEAD - DUP TAIL SWAP HEAD + Nlist DUP HEAD - * \GSLIST - ABS
> 2. / \>>
>
> the entire group of coordinates is "shifted" down to near the coordinate
> system origin for the purpose of the calculation.  This essentially does
> the same thing as moving the decimal place over.
>
> By calculating the same area without the "shift", I get an answer of:
> 12,055,387.000 m², which really is the same thing.  The program below is
> without the "shift":
>
> \<< Nlist DUP TAIL SWAP HEAD + Elist * \GSLIST Elist DUP TAIL SWAP HEAD
> + Nlist * \GSLIST - ABS 2. / \>>
>
> So the reason for this is that there are only 5 points within the
> polygon, if we add a bunch more, say add another point half way between
> each of the points already listed, maintaining the shape and area:
>
> 1 60,320.5280 1,395,020.7370
> 6 60,440.3020 1,397,445.8490
> 2 60,560.0760 1,399,870.9610
> 7 59,214.0345 1,400,449.6280
> 4 57,867.9930 1,401,028.2950
> 8 57,084.0585 1,399,294.5725
> 3 56,300.1240 1,397,560.8500
> 9 57,485.4985 1,397,700.7355
> 5 58,670.8730 1,397,840.6210
> 10 59,495.7005 1,396,430.6790
>
> Using the approach posted in 2009, I still get the same answer:
> 12,055,386.998 m²
>
> Now without the "shift", the answer is:
> 12,055,388.000 m²
>
> And it only gets worse from there on, the more points you add to the mix.
>
> You can see how the longer the list of high numbered coordinate values
> is, the larger a number \GSLIST will return, and with only 12 available
> digits you run out of digits pretty quickly.
>
> Jacob

Your "shift" code works very well. In an example from
the Wolf and Brinker text, the authors shift the axes
so that the most southerly and westerly stations are
equal to zero. The way the author presented it, I was
lead to believe that insufficient precision and serious
round-off errors occur unless double precision is used.
I think the problem is that due to how the coordinate
method was developed, unless your values are very
near to the coordinate system origin you will introduce
error into the calculations. Since this method uses
areas of individual trapezoids and triangles, maybe
the more oblique some of the vertices are, the more
error is introduced into the calculations.
 |  Next  |  Last
Pages: 1 2
Prev: RPN/RPL Calculator implementations, list of, regular post [long, FAQ]
Next: Graphics programming commands?