{"id":557,"date":"2009-11-24T15:45:59","date_gmt":"2009-11-24T14:45:59","guid":{"rendered":"http:\/\/sickel.net\/blogg\/?p=557"},"modified":"2012-06-27T12:10:55","modified_gmt":"2012-06-27T10:10:55","slug":"ecef-latlon-conversion","status":"publish","type":"post","link":"http:\/\/sickel.net\/blogg\/?p=557","title":{"rendered":"ECEF – lat\/lon conversion"},"content":{"rendered":"

I needed to do a ecef (earth centered earth fixed) to lat lon conversion in delphi. Following the formulaes found in Wikipedia<\/a> the following procedure was made:<\/p>\n

\r\nprocedure ecef2latlon(X,Y,Z : extended; var lat,lon,h : real);\r\n\/\/ returns lat with 6 digs accuracy, eg approx 0.1 m uncertainty.\r\n\/\/ Lon as accurate as the intrisic arctan2() function\r\n\r\nvar\r\n  r, e2, f, g,c,s,p,q,r0,u,v,z0,phi,lambda : extended;\r\nconst\r\n    \/\/ These constants are for wgs84\r\n   a=6378137.0;\r\n   b=6356752.3142;\r\n   ecc2=6.69437999014E-3; \/\/ \r\n   ecc22=6.73949674228E-3;\r\nbegin\r\n  r:=sqrt(x*x+y*y);\r\n  e2:=a*a-b*b;\r\n  f:=54.0*b*b*Z*Z;\r\n  g:=r*r+(1-ecc2)*Z*Z+ecc2*e2;\r\n  c:=ecc2*ecc2*f*r*r\/(g*g*g);\r\n  s:=power((1+c+sqrt(c*(c+2))),1\/3);\r\n  p:=f\/(3*power((s+1\/s+1),2)*g*g);\r\n  q:=sqrt(1+2*ecc2*ecc2*p);\r\n  r0:=-(p*ecc2*r)\/1+q+sqrt(0.5*a*a*(1+1\/q)-p*(1-ecc2)*z*z\/(q*(1+q))-0.5*p*r*r);\r\n  u:=sqrt(power((r-ecc2*r0),2)+z*z);\r\n  v:=sqrt(power((r-ecc2*r0),2)+(1-ecc2)*z*z);\r\n  z0:=b*b*z\/(a*v);\r\n  h:=u*(1-b*b\/(a*v));\r\n  phi:=arctan((z+ecc22*z0)\/r);\r\n  lambda:=arctan2(y,x);\r\n  lat:=phi\/pi*180;\r\n  lon:=lambda\/pi*180;\r\nend;\r\n<\/pre>\n
To check the result, I ported the code to bc:<\/p>\n
\r\ndefine cbrt(x) {\r\nreturn(e(l(x)\/3)) } \r\n\r\ndefine sgn(x) {\r\n        if (x == 0) {\r\n                return(0);\r\n        } else if (x < 0) {\r\n                return(-1);\r\n        } else if (x > 0) {\r\n                return(1);\r\n        }\r\n}\r\n\r\ndefine abs(x) {\r\n        if (x < 0) {\r\n                return(-1 * x);\r\n        } else {\r\n                return(x);\r\n        }\r\n}\r\n\r\ndefine atan2(y, x) {\r\n        auto pi, fi;\r\n        pi = 4 * a(1);\r\n        if (y == 0) {\r\n                if (x > 0) {\r\n                        return(0);\r\n                } else if (x == 0) {\r\n                        print \"undefined\\n\";\r\n                        halt;\r\n                } else if (x < 0) {\r\n                        return(pi);\r\n                }\r\n        }\r\n        fi = a(abs(y\/x));\r\n        if (x > 0) {\r\n                return(fi * sgn(y));\r\n        } else if (x == 0) {\r\n                return(pi * sgn(y) \/ 2);\r\n        } else if (x < 0) {\r\n                return((pi - fi) * sgn(y));\r\n        }\r\n}\r\n\r\n\r\npi=a(1)*4\r\n\r\nx=3147484.89\r\ny=589403.67\r\nz=5497533.82\r\n\r\na=6378137.00000000000000000;\r\nb=6356752.31420000000000000000;\r\necc2=6.694379990140000000000000E-3;\r\necc22=6.73949674228000000000000E-3;\r\n\r\n   r=sqrt(x*x+y*y);\r\n  e2=a*a-b*b;\r\n  f=54.0*b*b*z*z;\r\n  g=r^2+(1-ecc2)*z^2+ecc2*e2;\r\n  c=ecc2*ecc2*f*r*r\/(g*g*g);\r\n  s=cbrt(1+c+sqrt(c*(c+2)));\r\n  p=f\/(3*(s+1\/s+1)^2*g*g);\r\n  q=sqrt(1+2*ecc2*ecc2*p);\r\n  r0=-(p*ecc2*r)\/1+q+sqrt(0.5*a*a*(1+1\/q)-p*(1-ecc2)*z*z\/(q*(1+q))-0.5*p*r*r);\r\n  u=sqrt((r-ecc2*r0)^2+z*z);\r\n  v=sqrt((r-ecc2*r0)^2+(1-ecc2)*z*z);\r\n  z0=b*b*z\/(a*v);\r\n  h=u*(1-b*b\/(a*v));\r\n  phi=a((z+ecc22*z0)\/r);\r\n  lambda=atan2(y,x);\r\n  lat=phi\/pi*180;\r\n  lon=lambda\/pi*180;\r\n\r\nlat\r\nlon\r\n<\/pre>\n
Saving this as a file, eg ecef.bc and piping it through bc (i.e. bc -l < ecef.bc) should yield the result\n\n\n
\r\n59.94721465915728162420\r\n10.60647008394438896040\r\n<\/pre>\n
Using the above set x,y and z values in delphi should of cource yield the same result, give or take rounding errors.
\nDelphi (v7) and bc gives the same result to the 7th decimal place for latitude (59.9472147) that is approx mm presicion, which should be sufficient in most cases. A commersial, closed source, although mathlab based system I am using agrees with those two to the 5th decimal, eg 10cm presicion. As it is not possible to see what kind of data types are used in the commersial application, I would put the highest trust on bc. Although the code only has been tested in delphi, it is written in pretty general pascal and should work just as well for most pascal compilers.<\/p>\n
In php, I use the following function:<\/p>\n
\r\n\r\nfunction ecef2lla($x,$y,$z){\r\n\/\/ WGS84- parameters:\r\n$a = 6378137;\r\n$e = 8.1819190842622e-2;\r\n$b   = sqrt(pow($a,2)*(1-pow($e,2))); \/\/ ok 2517.0298622171\r\n$ep  = sqrt(($a*$a-$b*$b)\/pow($b,2)); \/\/ ok 2533.9931794567\r\n$p   = sqrt($x*$x+$y*$y); \/\/ ok\r\n$th  = atan2($a*$z,$b*$p);\r\n$lon = atan2($y,$x); \/\/ ok\r\n$lat = atan2(($z+($ep*$ep*$b*pow(sin($th),3))),($p-($e*$e*$a*pow(cos($th),3))));\r\n$N   = $a\/sqrt(1-$e*$e*pow(sin($lat),2));\r\n$alt = $p\/cos($lat)-$N;\r\n# return lon in range [0,2*pi)\r\n\/\/$lon = $lon%(2*pi());\r\n$lat=$lat\/pi()*180;\r\n$lon=$lon\/pi()*180;\r\n$ret['lat']=$lat;\r\n$ret['lon']=$lon;\r\n$ret['alt']=$alt;\r\nreturn $ret;\r\n}\r\n<\/pre>\n
This yields:
\nlat=59.947216007022
\nlon=10.606470083944<\/p>\n
I.e. about meter presicion  – maybe good enough, maybe not.<\/p>\n","protected":false},"excerpt":{"rendered":"
I needed to do a ecef (earth centered earth fixed) to lat lon conversion in delphi. Following the formulaes found in Wikipedia the following procedure was made: procedure ecef2latlon(X,Y,Z : extended; var lat,lon,h : real); \/\/ returns lat with 6 … Continue reading →<\/span><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_post_was_ever_published":false},"categories":[4],"tags":[],"class_list":["post-557","post","type-post","status-publish","format-standard","hentry","category-data"],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/pnVtD-8Z","jetpack_sharing_enabled":true,"jetpack_likes_enabled":true,"jetpack-related-posts":[],"_links":{"self":[{"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/posts\/557","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=557"}],"version-history":[{"count":7,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/posts\/557\/revisions"}],"predecessor-version":[{"id":559,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=\/wp\/v2\/posts\/557\/revisions\/559"}],"wp:attachment":[{"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=557"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=557"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/sickel.net\/blogg\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=557"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}