| Amazingly accurate UTM <-> Lat/Lon transformations [message #85905] |
Tue, 17 September 2013 01:50  |
tom.grydeland
Messages: 51
Registered: September 2012
|
Member |
|
|
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 >
|
|
|
|
| Re: Amazingly accurate UTM <-> Lat/Lon transformations [message #85909 is a reply to message #85905] |
Tue, 17 September 2013 11:16  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Tom Grydeland writes:
> 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.
Maybe it's just me, but it seems a little weird to me to claim that if
you use the software in a way that it is clearly not intended to be
used, you get bad results. Isn't this making a claim for clairvoyance on
the part of the developers?
I agree that IDL's map projection software is probably overdue for an
update. But, I'm not convinced this is any kind of a deal breaker.
Cheers,
David
P.S. For what it's worth, I, too, was confused by the initial reference
and had to spend 15 minutes or so figuring it out. That's about 14.5
minutes more than most readers of the newsgroup would have, I imagine.
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
|
|
|
| Re: Amazingly accurate UTM <-> Lat/Lon transformations [message #85915 is a reply to message #85909] |
Wed, 18 September 2013 02:15 |
tom.grydeland
Messages: 51
Registered: September 2012
|
Member |
|
|
On Tuesday, September 17, 2013 8:16:08 PM UTC+2, David Fanning wrote:
> Maybe it's just me, but it seems a little weird to me to claim that if
> you use the software in a way that it is clearly not intended to be
> used, you get bad results. Isn't this making a claim for clairvoyance on
> the part of the developers?
I'll admit to weird most Tuesdays and Saturdays, and the third Wednesday of the month, but clairvoyance I will leave to others (you know who you are!) :-)
> I agree that IDL's map projection software is probably overdue for an
> update. But, I'm not convinced this is any kind of a deal breaker.
My intention was not to throw dirt on the IDL mapping routines.
As I wrote (and you also point out), this is clearly abusing the UTM projection, and there are other projections more suited for this kind of application. I also wrote quite clearly that (depending on your requirements) the IDL routines produce acceptable results in the first and possibly second neighbouring zones.
That said, I am quite confident that I'm not the only one abusing UTM in this way, and indeed it turned out that my correspondent was also abusing UTM in exactly the same way, without being aware of the problems they were bringing upon themselves. Maybe it is the way of people to take what they know and see if it can be streched to fit over what they don't, what do you think?
So -- the information was intended for those who persist in doing this, even when told that it may not be the best of ideas, and to show them that the transformations _can_ be made in an invertible way, using the expressions from Wikipedia rather than IDL's mapping routines. Furthermore: IDL's routines will happily accept the task of transforming to UTM projections a quarter of a globe away, even if the morality of such action can be questioned. Given that the task will be performed, I think an accurate result is to be preferred over a wildly inaccurate one. Don't you agree?
> Cheers,
> David
>
> P.S. For what it's worth, I, too, was confused by the initial reference
> and had to spend 15 minutes or so figuring it out. That's about 14.5
> minutes more than most readers of the newsgroup would have, I imagine.
Maybe. I didn't find links to the particular section, and the question was meant to be about inverse hyperbolics anyway.
> David Fanning, Ph.D.
--T
|
|
|
|
|
|
| Re: Amazingly accurate UTM <-> Lat/Lon transformations [message #85920 is a reply to message #85914] |
Wed, 18 September 2013 06:19 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Fabien writes:
> I forwarded Tom's post to my colleagues, because many of them are not
> aware of these problems and Tom nicely showed the huge discrepancies
> that can occur when using UTM. Having said this, I am also concerned
> that even if the LonLat->EastingsNothings->LonLat transforms work better
> with Tom's code, it might still not be a good idea to use UTM for large
> maps.
>
The UTM projection is surely one of the most abused map projections in
history. Perhaps scientists should be made to take a quick course on it
in additional the the mandatory typing classes I've been advocating
forever.
> Regarding the map projections in IDL, it seems that IDL is far behind
> many other GIS softwares (correct me if I'm wrong):
> - the panel of "out of the box" projections is small
> - there is no datum shift transformation tool (!)
> - the accuracy, as shown by Tom, could be higher
>
> I believe this is due to the fact that Exelis wants to sell ENVI to IDL
> users in need for more complete/accurate projection engines...
Now, now, Fabien, I'm not dead yet. Please don't be usurping my role of
cynic on this newsgroup. ;-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: Amazingly accurate UTM <-> Lat/Lon transformations [message #85924 is a reply to message #85920] |
Wed, 18 September 2013 07:20 |
Andy Sayer
Messages: 127
Registered: February 2009
|
Senior Member |
|
|
On Wednesday, September 18, 2013 9:19:13 AM UTC-4, David Fanning wrote:
>
> The UTM projection is surely one of the most abused map projections in
>
> history. Perhaps scientists should be made to take a quick course on it
>
> in additional the the mandatory typing classes I've been advocating
>
> forever.
>
>
> David Fanning, Ph.D.
As a scientist and self-taught programmer, I completely agree with David that it would have been nice to formally have been taught programming, map projections (I know that I still don't know as much as I should about them), and many other things! As it is we learn as we go and sometimes I shudder to look back at old code.
I also agree with Tom, David etc that if IDL will do the coordinate conversion even when it is not 'sensible', it'd be good for it to be reasonably accurate, and/or have a big obvious warning in the documentation about that.
Andy
|
|
|
|
|
|
| Re: Amazingly accurate UTM <-> Lat/Lon transformations [message #85936 is a reply to message #85935] |
Thu, 19 September 2013 05:47 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Tom Grydeland writes:
>> My guess is that both of our efforts will
>> prove to be for naught for most users. :-)
>
> You're much too cynical.
Yeah, I think Coyote must have been typing. I was thinking "many" and
was shocked to see "most" when I re-read the article after it was
posted. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
