2

First, I needed to convert geographic coordinates (like those received from GPS) to global cartesian coordinates. I found an earlier question (How do I calculate a xyz-position of a gps-position relative to an other gps-position?) which worked perfectly. Now I need to go the other way around.

In Wikipedia, there are equations to convert between those, but the equations from geographic to cartesian are different, and not indicating which is latitude and which is longitude. I know it's because it's not geographic per se, and I also know that those angles have different constraints than those of latitude/longitude. This makes me not trust the reverse equations there. I searched elsewhere but I couldn't find the answer.

Precision is not important. Consider a perfect sphere.

Community
  • 1
GuiRitter
  • 141
  • 1
  • 1
  • 5
  • What software will you be using or are you trying to write the maths yourself? P.S. I haven't noticed any censoring on this site. – Michael Stimson Nov 03 '14 at 02:49

2 Answers2

2

The wiki reference you cite is based on a mathematical, not a geographic/cartographic, convention, as you say. It does however try to relate the geographic/cartographic equivalents. Here's my interpretation of the equivalencies:

Mathematics             Geography/Cartography
r, radial distance      R + h, radius + altitude
φ, polar angle          90 - φ, colatitude (latitude = 90 - colatitude)
θ, azimuthal angle      λ, longitude
radians                 degrees = radians * 360 / 2π

The conversion equations in the wiki article appear to be using the physics convention, not the mathematical one and confused me, too. (Maths and physics conventions have θ and φ swapped!) Another reference, mathworld.wolfram.com/SphericalCoordinates, helps us out.

The equations you seek are

r = sqrt(x*x + y*y + z*z)
h = r - R
λ = [arctan(y / x) * 360 / 2π]
  = [atan2(y, x) * 360 / 2π] ;; so as to account for quadrant or div-by-zero issues
φ = 90 - arccos(z / r) * 360 / 2π
  = arcsin(z / r) * 360 / 2π
Martin F
  • 8,948
  • 36
  • 58
2

I declared one point in each octant of the globe, transformed it to XYZ using the equations I had, and then tested martin f's answer. It didn't returned the same points. Then I delved deeper into Wikipedia's equations, and I finally understood how they worked. Then I adapted them.

latitude(r, x, y, z) = arcsin(z/r)(180/π)

longitude(r, x, y, z) = 
    if (x > 0) {
        arctan(y/x)(180/π)
    } else if (y > 0) {
        arctan(y/x)(180/π) + 180
    } else {
        arctan(y/x)(180/π) - 180
    }

(I wanted to write a piecewise function, but LaTeX isn't supported here) This is considering North and East as positive and right hand rule for XYZ, where Y is transverse to the plane that contains the Greenwhich meridian, positive to the East as well, and Z is aligned to North-South.

Thanks for helping me get there.

Edit: I just checked on Wikipedia, and it seems I learned the wrong right hand rule. Here, I mean X in the thumb, Y in the index finger and Z in the middle finger.

GuiRitter
  • 141
  • 1
  • 1
  • 5