Sunday 19 July 2015

coordinate system - Calculating areal distortion outside UTM zone?


One of my colleagues is working with data that is spread across two UTM zone. The majority of the data is in one zone, with a few outliers in another zone. He would like to know what the area distortion of those outliers would be if they were in the main UTM zone.


Is there a formula to calculate the areal distortion knowing how far into the other UTM zone the features were?



Answer



UTM uses a transverse Mercator projection with a scale factor of 0.9996 at the central meridian. In the Mercator, the distance scale factor is the secant of the latitude (one source: http://en.wikipedia.org/wiki/Mercator_projection), whence the area scale factor is the square of this scale factor (because it applies in all directions, the Mercator being conformal). Understanding the latitude as the spherical distance to the equator, and approximating the ellipsoid with a sphere, we can apply this formula to any aspect of the Mercator projection. Thus:



The scale factor is 0.9996 times the secant of the (angular) distance to the central meridian. The area scale factor is the square of this quantity.



To find this distance, consider the spherical triangle formed by traveling along a geodesic from an arbitrary point at (lon, lat) = (lambda, phi) straight towards the central meridian at longitude mu, along that meridian to the nearest pole, and then back along the lambda meridian to the original point. The first turn is a right angle and the second one is an angle of lambda-mu. The amount traveled along the last portion is 90-phi degrees. The Spherical Law of Sines applied to this triangle states




sin(lambda-mu) / sin(distance) = sin(90 degrees) / sin(90-phi)



with solution



distance = ArcSin(sin(lambda-mu) * cos(phi)).



This distance is given as an angle, which is convenient for computing the secant.


Example


Consider UTM zone 17, with central meridian at -183 + 17*6 = -81 degrees. Let the outlying location be at longitude -90 degrees, latitude 50 degrees. Then



Step 1: Spherical distance from (-90, 50) to the -81 degree meridian equals ArcSin(sin(9 degrees) * cos(50 degrees)) = 0.1007244 radians.


Step 2: The area distortion equals (0.9996 * sec(0.1007244 radians))^2 = 1.009406.


(Numerical calculations with the GRS 80 ellipsoid give the value as 1.009435, showing that the answer we computed is 0.3% too low: that's the same order of magnitude as the flattening of the ellipsoid, indicating the error is due to the spherical approximation.)


Approximations


To get a feel for how the area changes, we can use some trig identities to simplify the overall expression and expand that as a Taylor series in lambda-mu (the displacement between the point's longitude and the longitude of the UTM central meridian). It works out to



Area scale factor ~ 0.9992 * (1 + cos(phi)^2 * (lambda-mu)^2).



As with all such expansions, the angle lambda-mu must be measured in radians. The error is less than 0.9992 * cos(phi)^4 * (lambda-mu)^4, which is close to the square of the difference between the approximation and 1--that is, the square of the value after the decimal point.


In the example with phi = 50 degrees (with a cosine of 0.642788) and lambda-mu = -9 degrees = -0.15708 radians, the approximation gives 0.9992 * (1 + 0.642788^2 * (-0.15708)^2) = 1.009387. Looking past the decimal point and squaring, we deduce (even without knowing the correct value) that its error cannot be greater than (0.009387)^2 = less than 0.0001 (and in fact the error is only one-fifth that size).



From this analysis it is evident that at high latitudes (where cos(phi) is small), scale errors will always be small; and at lower latitudes, area scale errors will behave like the square of the difference in longitudes.


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...