Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
Amazingly accurate UTM <-> Lat/Lon transformations
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
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.")
On 09/17/2013 08:16 PM, David Fanning wrote: > 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. 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. 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...
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
Tom Grydeland writes: > 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? I agree that wildly inaccurate results should be avoided. ;-) I'm not totally convinced the way to avoid them in this case is to roll your own software. I'd rather see folks fix the problem by using map projection software correctly. My guess is that both of our efforts will prove to be for naught for most users. :-) 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.")
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.")
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
On Wednesday, September 18, 2013 3:04:13 PM UTC+2, David Fanning wrote: > I'm not totally convinced the way to avoid them in this case is to roll > your own software. I'd rather see folks fix the problem by using map > projection software correctly. In other words, don't use UTM for this kind of task. I don't disagree. On the other hand, I cannot demonstrate how much more accurate these formulas are without implementing them (i.e. rolling my own), and the improvements are not available to anyone else without rolling their own until Exelis takes notice. > My guess is that both of our efforts will > prove to be for naught for most users. :-) You're much too cynical. The Coyote library is a great benefit to all those using it, and my little note has enlightened at least one reader. Trying to be of benefit to _most_ users is a tall order. > David --T
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.")