EPSG
EPSG guidance note #7-2, http://www.epsg.org
2007-11-01
If the natural origin of the projection is at the topocentric origin, this is a special case of the Vertical Perspective (orthographic case) (method code 9839) in which the ellipsoid height of all mapped points is zero (h = 0).
Note: These formulas have been transcribed from EPSG Guidance Note #7-2. Users are encouraged to use that document rather than the text which follows as reference because limitations in the transcription will be avoided. The Orthographic Projection forward conversion from 2D geographic coordinates latitude and longitude (lat, lon) and the origin on the ellipsoid (latO, lonO) is given by: E = FE + nu cos(lat) sin (lon – lonO) N = FN + nu [sin(lat) cos(latO) – cos(lat) sin(latO) cos (lon – lonO)] + e^2 (nuO sin(latO) – nu sin(lat)) cos(latO) where nu is the prime vertical radius of curvature at latitude lat; nu = a /(1 – e^2 sin^2(lat))^0.5, nuO is the prime vertical radius of curvature at latitude of origin latO; nuO = a /(1 – e^2 sin^2(latO)^0.5, e is the eccentricity of the ellipsoid and e^2 = (a^2 – b^2)/a^2 = 2f – f^2 a and b are the ellipsoidal semi-major and semi-minor axes, 1/f is the inverse flattening, and the latitude and longitude of the projection origin are latO and lonO. These formulas are similar to those for the orthographic case of the vertical perspective (method code 9839) except that, for the Orthographic Projection given here, h = 0 and the term (nu + h) reduces to nu. The projection origin is at the topocentric system origin latO, lonO with false origin coordinates FE and FN. For the reverse formulas for latitude and longitude corresponding to a given Easting (E) and Northing (N), iteration is required as the prime vertical radius (nu) is a function of latitude. Begin by seeding the iteration with the center of projection (or some better guess): lat = latO lon = lonO Enter the iteration here with the (next) best estimates of lat and lon. Then solve for the radii of curvature in the prime vertical (nu) and meridian (rho): nu = a / (1 – e^2 sin^2(lat))^0.5 rho = a (1 – e^2) / (1 – e^2 sin^2(lat))^1.5 Compute test values of E and N (E' and N') using the forward equations: E' = FE + nu cos(lat) sin(lon – lonO) N' = FN + nu [sin(lat) cos(latO) – cos(lat) sin(latO) cos(lon – lonO)] + e^2 (nuO sin(latO) – nu sin(lat) ) cos(latO) Partially differentiate the forward equations to solve for the elements of the Jacobian matrix: J11 = dE/dlat = – rho sin(lat) sin(lon – lonO) J12 = dE/dlon = nu cos(lat) cos(lon – lonO) J21 = dN/dlat = rho [cos(lat) cos(latO) + sin(lat) sin(latO) cos(lon – lonO)] J22 = dN/dlon = nu sin(latO) cos(lat) sin (lon – lonO) Solve for the determinant of the Jacobian: D = J11 J22 – J12 J21 Solve the northerly and easterly differences this iteration: dE = E – E' dN= N – N' Adjust the latitude and longitude for the next iteration by inverting the Jacobian and multiplying by the differences: lat = lat + (J22 dE – J12 dN) / D lon = lon + (–J21 dE + J11 dN) / D Return to the entry point with new estimates of latitude and longitude and iterate until the change in lat and lon is not significant.