EPSG

EPSG guidance note #7-2, http://www.epsg.org

2022-01-06

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. For the calculation of easting and northing from latitude and longitude, first calculate constants for the projection: n = f / (2-f) B = [a/(1+n)] (1 + n^2/4 + n^4/64) h1 = n/2 – (2/3)n^2 + (5/16)n^3 + (41/180)n^4 h2 = (13/48)n^2 – (3/5)n^3 + (557/1440)n^4 h3 = (61/240)n^3 – (103/140)n^4 h4 = (49561/161280)n^4 Then the meridional arc distance from equator to the projection origin (Mo) is computed from: If LatO = 0 then Mo = 0 else if LatO ≡ 90°N ≡ π/2 radians Mo = B (π/2) else if LatO ≡ 90°S ≡ -π/2 radians Mo = B (-π/2) else Qo = asinh(tan LatO) – [e atanh(e sin LatO)] βo = atan(sinh Qo) ξO0 = asin (sin βo) Note: The previous two steps are taken from the generic calculation flow given below for latitude Lat, but here for LatO may be simplified to ξO0 = βo = atan(sinh Qo). ξO1 = h1 sin(2ξOo) ξO2 = h2 sin(4ξOo) ξO3 = h3 sin(6ξOo) ξO4 = h4 sin(8ξOo) ξO = ξO0+ ξO1+ ξO2+ ξO3+ ξO4 Mo = B ξO end Note: if the projection grid origin is very close to but not exactly at the pole (within 2" or 50m), in the equation above for Qo the tangent function is unstable and may fail. Rather than using Qo as above, Mo may instead be calculated from: Mo = a[(1 – e^2/4 – 3e^4/64 – 5e^6/256 –....)LatO – (3e^2/8 + 3e^4/32 + 45e^6/1024+....)sin(2.LatO) + (15e^4/256 + 45e^6/1024 +.....)sin(4.LatO) – (35e^6/3072 + ....)sin(6.LatO) + .....] with LatO in radians. Then Q = asinh(tan Lat) – [e atanh(e sin Lat)] β = atan(sinh Q) η0 = atanh [cos β sin(Lon – LonO)] ξ0 = asin (sin β cosh η0) ξ1 = h1 sin(2ξ0) cosh(2η0) η1 = h1 cos(2ξ0) sinh(2η0) ξ2 = h2 sin(4ξ0) cosh(4η0) η2 = h2 cos(4ξ0) sinh(4η0) ξ3 = h3 sin(6ξ0) cosh(6η0) η3 = h3 cos(6ξ0) sinh(6η0) ξ4 = h4 sin(8ξ0) cosh(8η0) η4 = h4 cos(8ξ0) sinh(8η0) ξ = ξ0 + ξ1 + ξ2 + ξ3 + ξ4 η = η0 + η1 + η2 + η3 + η4 and Easting, E = FE + ko B η Northing, N = FN + ko (B ξ – Mo) For the reverse formulas to convert Easting and Northing projected coordinates to latitude and longitude first calculate constants of the projection where n is as for the forward conversion, as are B and Mo: h1' = n/2 – (2/3)n^2 + (37/96)n^3 – (1/360)n^4 h2' = (1/48)n^2 + (1/15)n^3 – (437/1440)n^4 h3' = (17/480)n^3 – (37/840)n^4 h4' = (4397/161280)n^4 Then η' = (E – FE) / (B ko) ξ' = [(N – FN) + ko Mo] / (B ko) ξ1' = h1' sin(2ξ') cosh(2η') η1' = h1' cos(2ξ') sinh(2η') ξ2' = h2' sin(4ξ') cosh(4η') η2' = h2' cos(4ξ') sinh(4η') ξ3' = h3' sin(6ξ') cosh(6η') η3' = h3' cos(6ξ') sinh(6η') ξ4' = h4' sin(8ξ') cosh(8η') η4' = h4' cos(8ξ') sinh(8η') ξ0' = ξ' – (ξ1' + ξ2' + ξ3' + ξ4') η0' = η' – (η1' + η2' + η3' + η4') β' = asin(sin ξ0' / cosh η0') Q' = asinh(tan β') Q" = Q' + [e atanh(e tanh Q')] = Q' + [e atanh(e tanh Q")] which should be iterated until the change in Q" is insignificant. Then Lat = atan(sinh Q") Lon = LonO + asin(tanh(η0') / cos β')

For Projected Coordinate System OSGB36 / British National Grid Parameters: Ellipsoid Airy 1830 a = 6377563.396 m 1/f = 299.32496 then e'^2 = 0.00671534 and e^2 = 0.00667054 Latitude of natural origin (LatO) = 49°00'00"N = 0.85521133 rad Longitude of natural origin (LonO) = 2°00'00"W = -0.03490659 rad Scale factor (ko) = 0.9996013 False Eastings (FE) = 400000.00 m False Northings (FN) = -100000.00 m Forward calculation for: Latitude = 50°30'00.00"N = 0.88139127 rad Longitude = 00°30'00.00"E = 0.00872665 rad Constants of the projection: n = 0.00167322 B = 6366914.609 h1 = 0.0008347452 h2 = 0.0000007554 h3 = 1.18487E-09 h4 = 2.40864E-12 QO = 0.9787671618 ξO0 = 0.8518980373 ξO1 = 0.0008273732 ξO2 = -0.0000001986 ξO3 = -1.0918E-09 ξO4 = 1.2218E-12 Mo = 5429228.602 Q = 1.0191767215 β = 0.8781064142 η0 = 0.0278629616 ξ0 = 0.8785743280 η1 = -0.0000086229 ξ1 = 0.0008215669 η2 = -0.0000000786 ξ2 = -0.0000002768 η3 = 1.05551E-10 ξ3 = -1.01855E-09 η4 = 3.97791E-13 ξ4 = 1.67447E-12 η = 0.0278542603 ξ = 0.8793956171 Then Easting E = 577274.99 metres Northing N = 69740.50 metres Reverse calculation for same easting and northing first gives: h1' = 0.0008347455 h2' = 0.0000000586 h3' = 1.65563E-10 h4' = 2.13692E-13 Then ξ' = 0.87939562 η' = 0.0278542603 ξ1' = 0.0008213109 η1' = -0.0000086953 ξ2' = -0.0000000217 η2' = -0.0000000061 ξ3' = -1.41881E-10 η3' = 1.486E-11 ξ4' = 1.49609E-13 η4' = 3.50657E-14 ξ0' = 0.8785743280 η0' = 0.0278629616 β' = 0.8781064142 Q' = 1.0191767215 Q" 1st iteration = 1.0243166838 Q" 2nd iteration = 1.0243306667 Q" 3rd iteration = 1.0243307046 Q" 4th iteration = 1.0243307047 Then Latitude (Lat) = 50°30'00.000"N Longitude (Lon) = 00°30'00.000"E