EPSG
"Map Projections - A Working Manual" by John P. Snyder, USGS Professional Paper 1395.
2023-06-29
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 formulas to convert geodetic latitude and longitude (lat, lon) to Easting (E) and Northing (N) are: Easting (E) = EF + (rho . sin(theta)) Northing (N) = NF + rhoO – (rho . cos(theta)) where theta = n . (lon - lonO) rho = [a . (C – n.alpha)^0.5] / n rhoO = [a . (C – n.alphaO)^0.5] / n and C = m1^2 + (n . alpha1) n = (m1^2 – m2^2) / (alpha2 - alpha1) m1 = cos lat1 / (1 – e^2 sin^2(lat1))^0.5 m2 = cos lat2 / (1 – e^2 sin^2(lat2))^0.5 alpha = (1 – e^2) . {[sin(lat) / (1 – e^2 sin^2(lat))] – [1/(2e)] . ln [(1 – e sin(lat)) / (1 + e sin(lat))]} alphaO = (1 – e^2) . {[sin(latO) / (1 – e^2 sin^2(latO))] – [1/(2e)] . ln [(1 – e sin(latO)) / (1 + e sin(latO))]} alpha1 = (1 – e^2) . {[sin(lat1) / (1 – e^2 sin^2(lat1))] – [1/(2e)] . ln [(1 – e sin(lat1)) / (1 + e sin(lat1))]} alpha2 = (1 – e^2) . {[sin(lat2) / (1 – e^2 sin^2(lat2))] – [1/(2e)] . ln [(1 – e sin(lat2)) / (1 + e sin(lat2))]} The reverse formulas to derive the geodetic latitude and longitude of a point from its Easting and Northing values are: lat = ß' + [(e^2/3 + 31e^4/180 + 517e^6/5040) . sin 2ß'] + [(23e^4/360 + 251e^6/3780) . sin 4ß'] + [(761e^6/45360) . sin 6ß'] lon = lonO + (theta / n) where C, n and rhoO are as in the forward equations, ß' is the authalic latitude, and: ß' = asin(alpha' / {1 – [(1 – e^2) / 2e] . ln [(1 – e) / (1 + e)]}) alpha' = [C – (rho'^2 . N^2 / a^2)] / n rho' = {(E – EF)^2 + [rhoO – (N – NF)]^2 }^0.5 If n is positive, theta = atan2 {(E – EF) , [rhoO – (N – NF)]} but if n is negative, the signs of both arguments of the atan2 function must be reversed: theta = atan2 {– (E – EF), – [rhoO – (N – NF)]} (see implementation notes in GN7-2 preface for atan2 convention).
Example 1 (Northern Hemisphere): For Projected Coordinate Reference System: NAD83 / Great Lakes Albers (EPSG CRS code 3174) Parameters: Ellipsoid: GRS 1980 a = 6378137.00 metres 1/f = 298.2572221 then e = 0.081819191 and e^2 = 0.00669438 Latitude of false origin = 45°34'08.3172"N = 0.48578331 rad Longitude of false origin = 84°27'21.4380"W = -1.72787596 rad Latitude of 1st standard parallel = 42°07'21.9864"N = 0.49538262 rad Latitude of 2nd standard parallel = 49°00'54.6480"N = 0.52854388 rad Easting at false origin = 1000000.000 metre Northing at false origin = 1000000.000 metre Forward calculation for: Latitude = 42°45'00.000"N = 0.746128255 rad Longitude = 78°45'00.000"W = -1.374446786 rad first gives : alpha = 1.351293970 alpha0 = 1.4218650778 alpha1 = 1.3351453325 alpha2 = 1.5034868510 m1 = 0.7428286888 m2 = 0.6571136237 n = 0.7128137342 C = 1.5035043911 rho = 6577011.9827 rho0 = 6263350.4332 theta = 0.0709874815 Then Easting = 1466493.492 metre Northing = 702903.006 metre Reverse calculation for same easting and northing first gives: theta' = 0.0709874815 alpha' = 1.3512939695 rho' = 6577011.983 beta' = 0.7438962839 Then Latitude = 42°45'00.000"N Longitude = 78°45'00.000"W Example 2 (Southern Hemisphere): Parameters: Ellipsoid: GRS 1967 Modified a = 6378160.0 metres 1/f = 298.25 then e = 0.081820180 and e^2 = 0.006694542 Latitude of false origin = 32°00'00.000"S = -0.5585053606 rad Longitude of false origin = 60°00'00.000"W = -1.0471975512 rad Latitude of 1st standard parallel = 5°00'00.000"S = -0.0872664626 rad Latitude of 2nd standard parallel = 42°00'00.000"S = -0.7330382858 rad Easting at false origin = 0.000 metre Northing at false origin = 0.000 metre Forward calculation for: Latitude = 18°30'02.016"S = -0.322895686 rad Longitude = 46°00'01.538"W = -0.802858912 rad first gives : alpha = -0.630662757 alpha0 = -1.054065016 alpha1 = -0.173150420 alpha2 = -1.331965641 m1 = 0.9962200286 m2 = 0.7442610814 n = -0.378429434 C = 1.0579795608 rho = -15255863.89 rho0 = -13683051.84 theta = -0.092464933 Then Easting = 1408623.196 metre Northing = 1507641.482 metre Reverse calculation for same easting and northing first gives: (Since n is negative, we reverse the signs of both atan2 arguments to calculate theta') theta' = -0.092464933 alpha' = -0.630662757 rho' = 15255863.888 beta' = -0.321550075 Then Latitude = 18°30'02.016"S Longitude = 46°00'01.538"W