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 befor example, if state plane coordinates are used (see Chapter 21). In those cases, to ensure sufficient precision and prevent serious roundoff 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 crisscross, 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 crisscross, when joining BCEC 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 crisscross, 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 crisscross, when joining BCEC 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 crisscross 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 roundoff 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? 