Hi all,
In this newsgroup, I recently inquired as to whether inverse hyperbolics
would be available in IDL anytime soon. I added a link to the Wikipedia
entry on UTM (The Universal Transverse Mercator map projection) by way of
explanation. I got an exasperated reply:
"If you want this post to be usefull to at least one more person (I cant
talk only for me) you might want to give more information than a link to
a wikipedia page. I am using UTM (in IDL) and I am asking myself: what's
the point?"
To my mind, the question (on inverse hyperbolics) was self-explanatory. My
pocket calculator is older than IDL, and it has inverse hyperbolics. It's
just weird that IDL doesn't. (I'm glad to see that Mark Piper agrees.)
If you, like me, abuse the UTM coordinate system to map regions that cover much more than the intended 6 degrees of latitude (or 12 degrees in zones 31-37X around Svalbard), and you want your conversions to be reversible (in case you should ever need to reproject), then you should be aware that the projections performed by IDL's map system are _not_ reversible without significant errors outside their target zone.
In our application, we have used projection in UTM zone 33 (centered at 15°E)
to cover most of Europe. Iceland is the westmost part, extending to
westward of 21°W, and belonging in UTM zone 27.
My benchmark point in the tables below is Surtsey (63.303°N, 20.6047°W), which
is close to the middle of zone 27. I use the IDL routines MAP_PROJ_INIT,
MAP_PROJ_FORWARD and MAP_PROJ_INVERSE to convert these coordinates to UTM and
back to Lat/Lon. Whatever discrepancy arises is multiplied with a Earth
circumference of 4e7 metres and the longitude scaled by COS(63.303°) to give the
round-trip error in metres. (This computation is just for the scale of the
error, no need to drag in geodesy for that).
Here is my IDL code for this (and if I'm doing something wrong, PLEASE let
me know!):
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
zones = 27 + indgen(13) - 6
lonlat = [-20.6047d0, 63.303d0] ; Surtsey
scale = 1e7/90 * [cos(!dtor*lonlat[0]), 1d0]
print, 'L/L -> UTM -> L/L roundtrip error, IDL'
print, '******'
if 1 then $
foreach zone, zones do begin
map = map_proj_init('UTM', zone=zone, ellipsoid=12)
uco = map_proj_forward(lonlat[0], lonlat[1], map=map)
lola = map_proj_inverse(uco[0], uco[1], map=map)
print, format='("zone ", i3, " lon error: ", f0.4, " m, lat error: ", f0.4, " m")', $
zone, (lonlat - lola)*scale
endforeach
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
And the results:
L/L -> UTM -> L/L roundtrip error, IDL
******
zone 21 lon error: -55646.5530 m, lat error: 4697.2484 m
zone 22 lon error: -12510.4027 m, lat error: 889.8227 m
zone 23 lon error: -2241.8025 m, lat error: 129.2282 m
zone 24 lon error: -271.5019 m, lat error: 11.9819 m
zone 25 lon error: -15.5833 m, lat error: 0.4807 m
zone 26 lon error: -0.1433 m, lat error: 0.0027 m
zone 27 lon error: -0.0000 m, lat error: 0.0000 m
zone 28 lon error: 0.0566 m, lat error: 0.0010 m
zone 29 lon error: 9.7307 m, lat error: 0.2835 m
zone 30 lon error: 196.7202 m, lat error: 8.3330 m
zone 31 lon error: 1744.4105 m, lat error: 97.4464 m
zone 32 lon error: 10130.8167 m, lat error: 702.6532 m
zone 33 lon error: 46134.3743 m, lat error: 3814.3222 m
So in "home zone" and its two neighbours, the error is within 15 cm, but
it increases rapidly to unusability. In zone 33 (which we use) the
round-trip error is close to 50 km!
I have implemented the expressions from Wikipedia, in IDL, with the
following results on exactly the same test (notice that the scale of the
error is now micrometers):
L/L -> UTM -> L/L roundtrip error, Karney
******
zone 21 lon error: -342.8242 um, lat error: 448.4638 um
zone 22 lon error: -258.8382 um, lat error: 455.0497 um
zone 23 lon error: -189.8761 um, lat error: 456.3019 um
zone 24 lon error: -133.0246 um, lat error: 455.2139 um
zone 25 lon error: -84.9120 um, lat error: 453.5347 um
zone 26 lon error: -42.3714 um, lat error: 452.2083 um
zone 27 lon error: -2.5883 um, lat error: 451.6786 um
zone 28 lon error: 37.0292 um, lat error: 452.0868 um
zone 29 lon error: 79.0518 um, lat error: 453.3294 um
zone 30 lon error: 126.2577 um, lat error: 455.0000 um
zone 31 lon error: 181.7635 um, lat error: 456.2474 um
zone 32 lon error: 248.9413 um, lat error: 455.4381 um
zone 33 lon error: 330.8319 um, lat error: 449.7830 um
I'm doing nothing except what is in the Wikipedia entry here, and I'm using
Mike Galloy's MG_ATANH implementation of tanh^{-1}.
Adding terms O(n^4) (The coefficients for n^4 terms can be found in papers
Karney [2011] and Kawase [2011, 2013]) improves the accuracy further:
L/L -> UTM -> L/L roundtrip error, Karney
******
zone 21 lon error: 1.0656 um, lat error: -2.9938 um
zone 22 lon error: 0.6832 um, lat error: -2.8951 um
zone 23 lon error: 0.4382 um, lat error: -2.8106 um
zone 24 lon error: 0.2768 um, lat error: -2.7395 um
zone 25 lon error: 0.1652 um, lat error: -2.6969 um
zone 26 lon error: 0.0794 um, lat error: -2.6685 um
zone 27 lon error: 0.0041 um, lat error: -2.6606 um
zone 28 lon error: -0.0691 um, lat error: -2.6661 um
zone 29 lon error: -0.1530 um, lat error: -2.6914 um
zone 30 lon error: -0.2598 um, lat error: -2.7324 um
zone 31 lon error: -0.4131 um, lat error: -2.7987 um
zone 32 lon error: -0.6451 um, lat error: -2.8848 um
zone 33 lon error: -1.0050 um, lat error: -2.9859 um
In my mind, this suggests that the people who work on map projections for
IDL would spend their time usefully if they were to take a look at the
formulas given in Wikipedia for the conversion between Lat/Lon and UTM.
Also anyone who uses UTM with IDL, especially if they abuse the UTM
projection well outside its intended area (instead of switching to a
projection more suited for large regions) should have a look at the
Wikipedia article on UTM and its references, and try implementing these
formulas.
--T (tom snail norut period no)
References:
Karney, C. F. F. (2011). Transverse Mercator with an accuracy of a few nanometers. Journal of Geodesy, 85(8), 475–485. doi:10.1007/s00190-011-0445-3
<URL: http://arxiv.org/abs/1002.1417 >
KAWASE, K. (2011). A General Formula for Calculating Meridian Arc Length and its Application to Coordinate Conversion in the Gauss-Krüger Projection. Bulletin of the Geospatial Information Authority of Japan,, 59, 1–13.
<URL: www.gsi.go.jp/common/000062452.pdf >
KAWASE, K. (2013). Concise Derivation of Extensive Coordinate Conversion Formulae in the Gauss-Krüger Projection. Bulletin of the Geospatial Information Authority of Japan, 60, 1–6.
<URL: www.gsi.go.jp/common/000065826.pdf >
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
