Thursday, 2 January 2020

scale - When does a local system introduce errors in a map due to the earth's curvature?


If the extent of a map is very small, like a building plan, there is no coordinate reference system used.


To which scale/size of building is that acceptable? Of course it depends on the precision requirements, but maybe there is a general rule of thumb or some fixed values?



Answer



This answer will derive a simple, useful rule of thumb. Its analysis suggests that within any region of diameter a, the relative error when using a suitably well-chosen map projection to measure distances within that region will typically be a/(2R) or less, where R is the earth's radius. Since R is large (over 6 million meters or 4 thousand miles), projections can yield accurate information for demanding studies over relatively large areas.





The amount of linear measurement error in a map depends on the projection used. However, there are some basic limits imposed by geometry.


To get a feel for these limits, consider a map that is concerned with representing point-to-point distances within a large triangular area, such as one covering the North Pole and two points on the Equator, say (0,0) and (90,0) (in longitude-latitude). Because all three distances are very nearly the same--they each subtend 90 degrees of arc--this is an equilateral triangle, with sides that are 10 megameters long. On a low-error map, it would therefore have to be drawn as an equilateral triangle ABC.


Figure


Consider a point A' near the middle of BC. That makes its (scaled) distances to B and C around 10/2 = 5 megameters, exactly as shown on the map. But A' happens to lie on the Equator (because BC is part of the Equator), so AA' has a distance of 10 megameters, too. The Euclidean distance on the map, though, is only around 8.7 megameters. That's a -13% error.


These discrepancies between the predictions of spherical and Euclidean geometries are why maps must err in representing some distances. (We needn't worry about the complexities of an ellipsoidal earth model because anything we conclude about a spherical model will be within about 0.3% of the truth.)


We can use a similar construction to get a feeling for the amount of error that must occur in smaller equilateral triangles. Since the right side of ABC is just the mirror image of its left side, we might as well focus on the right triangle ABA'. The inaccuracy in the projection can be related to the amount by which the Pythagorean Theorem fails: (AB)^2 is not exactly (A'B)^2 + (AA')^2 on the sphere. It really is


cos(AB) = cos(A'B)*cos(AA')

where the three spherical distances AB, A'B, and AA' are measured in terms of arclength (in radians).



(Don't worry--no trigonometry is going to be committed here. These cosines will be replaced by simple approximations in short order.)


The relationship between arclength on a spherical world of radius R and geodetic length on the surface is simple: multiply the arclength by R. To this end, then, consider a small spherical right triangle with legs of lengths a and b. Let the geodetic length of its hypotenuse be c. On the sphere the angles are a/R and b/R, whence


cos(c/R) = cos(a/R) * cos(b/R).

Since these lengths are presumed small, we may approximate the cosines of the very small angles they subtend as


cos(x) ~ 1 - x^2/2

with extremely high accuracy.


Plugging in this approximation gives


1 - (c/R)^2/2 ~ (1 - (a/R)^2/2)(1 - (b/R)^2/2)


with the solution


c^2 = a^2 + b^2 - 1/2 * (a*b/R)^2.

The first two terms on the right are what Euclidean geometry says c^2 ought to be: they would be the squared distance shown on a projected map at 1:1 scale. The last term, -(ab/r)^2/2, therefore can be considered a measure of the error in the squared distance. (Because it is negative, we learn that the tips of the legs are actually closer on the sphere than they appear on the map. The sphere's positive curvature is drawing the legs inward toward each other.)


Because measurement errors are relative--a one meter error is huge in a building plan and inconsequential in a map of a continent--let's express the size of this error relative to the length of the hypotenuse on the map:


(relative error)^2 =  1/2 * (a*b/R)^2 / (a^2 + b^2);

relative error = (ab/sqrt(2(a^2 + b^2))) / R


Because these values will tend to be small, let's multiply them by 10^6 to express them as parts per million (ppm). In the case of two equal legs, where the error shows up the most,


relative error = (a^2/sqrt(4a^2)) / R = a/(2R) = 10^6 a/(2R) ppm.

This ought to be a good guide concerning how the errors in the best possible projections will tend to scale with the size of the projected region. Since the WGS authalic radius is R = 6.3710072 megameters, it suggests we should expect no more than 1 ppm relative error for every 12.74 meters (about 40 feet) covered by a good projection.


For reference, a coordinate system commonly used for relatively high-accuracy calculations is UTM, which has a built-in error of 400 ppm along its central meridians. Compared to that inaccuracy, we needn't start worrying about the errors in any good projected map until the distances exceed 400 * 12.74 meters = 5 kilometers or so. That certainly covers maps of most buildings. At almost ten times this distance, around 40 kilometers, the error estimate first reaches 3000 ppm. That is the discrepancy between a spherical earth model and an ellipsoidal one. This, then, should be the scale at which the two earth models just barely begin to diverge visibly.


For larger distances this error approximation is too conservative: the quadratic approximations to the cosines no longer hold. For instance, at 10 megameters--as in the original example--the formula gives 10/12.74 = 78% (too low), whereas we saw that the error was likely around -13% or so. Nevertheless, the error estimate is of the right order of magnitude. It's a useful rule of thumb.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...