EPSG
Land Survey Office (ZU), Prague. www.cuzk.cz/zu
2018-08-29
Incorporates a polynomial transformation which is defined to be exact and for practical purposes is considered to be a map projection.
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. From the defining parameters the following constants for the projection may be calculated : A = a(1 - e^2)^0.5 / [1 - e^2 sin^2.(latC)] B = {1 + [e^2 * cos^4(latC) / (1 - e^2)]}^0.5 gammao = asin[sin(latC) / B] to = tan(pi/4 + gammao/2).[(1 + e sin(latC)) / (1 - e sin(latC))]^(e.B/2) / [tan(pi/4 + latC/2)]^B n = sin(latp) ro = kp.A / tan(latp) To derive the projected Southing and Westing coordinates of a point with geographical coordinates (lat, lon) the formulas for the Krovak are: U = 2(atan{to.tan^B(lat/2 + pi/4) / [(1 + e sin(lat)) / (1 - e sin(lat))]^[e.B/2]} - pi/4) V = B(lonO - lon) where lonO and lon must both be referenced to the same prime meridian. T = asin[cos(alphaC).sin(U) + sin(alphaC).cos(U). cos(V)] D = asin[cos(U).sin(V)/cos(T)] theta = n.D r = ro.tan^n(pi/4 + latp/2) / tan^n(T/2 + pi/4) Xp = r.cos(theta) Yp = r.sin(theta) Xr = Xp – Xo Yr = Yp – Yo dX = C1 + C3.Xr – C4.Yr – 2.C6.Xr.Yr + C5.(Xr^2 – Yr^2) + C7.Xr.(Xr^2 – 3.Yr^2) – C8.Yr.(3.Xr^2 – Yr^2) + 4.C9.Xr.Yr.(Xr^2 – Yr^2) + C10.(Xr^4 + Yr^4 – 6.Xr^2.Yr^2) dY = C2 + C3.Yr + C4.Xr + 2.C5.Xr.Yr + C6.(Xr^2 – Yr^2) + C8.Xr.(Xr^2 – 3.Yr^2)+ C7.Yr.(3.Xr^2 – Yr^2) – 4.C10.Xr.Yr.(Xr^2 – Yr^2) + C9.(Xr^4 + Yr^4 – 6.Xr^2.Yr^2) Southing X = FN + Xp – dX Westing Y = FE + Yp – dY The reverse formulas to derive the latitude and longitude of a point from its Southing and Westing values are: Xr' = (Southing – FN) – Xo Yr' = (Westing – FE) – Yo dX' = C1 + C3.Xr' – C4.Yr' – 2.C6.Xr'.Yr' + C5.(Xr'^2 – Yr'^2) + C7.Xr'.(Xr'^2 – 3.Yr'^2) – C8.Yr'.(3.Xr'^2 – Yr'^2) + 4.C9.Xr'.Yr'.(Xr'^2 – Yr'^2) + C10.(Xr'^4 + Yr'^4 – 6.Xr'^2.Yr'^2) dY' = C2 + C3.Yr' + C4.Xr' + 2.C5.Xr'.Yr' + C6.(Xr'^2 – Yr'^2) + C8.Xr'.(Xr'^2 – 3.Yr'^2) + C7.Yr'.(3.Xr'^2 – Yr'^2) - 4.C10.Xr'.Yr'.(Xr'^2 – Yr'^2) + C9.(Xr'^4 + Yr'^4 – 6.Xr'^2.Yr'^2) Xp' = (Southing – FN) + dX' Yp' = (Westing – FE) + dY' r' = [(Yp')^2 + (Xp')^2]^(1/2) theta' = atan2[Yp' , Xp'] (see GN7-2 implementation notes in preface for atan2 convention) D' = theta' / sin(latp) T' = 2{atan[((ro / r')^(1/n)).tan(pi/4 + latp/2)] - pi/4} U' = asin[cos(alphaC).sin(T') - sin(alphaC).cos(T').cos(D')] V' = asin(cos(T').sin(D') / cos(U')) Then latitude lat is found by iteration using U' as the value for lat(j-1) in the first iteration: lat(j) = 2*(atan{tO^(-1/B) tan^(1/B).(U'/2 + pi/4).[(1 + e sin(lat(j-1)) / (1 - e sin(lat(j-1))]^(e/2)} - pi/4) Then lon = lonO - V' / B where lon is referenced to the same prime meridian as lonO.
For Projected Coordinate Reference System: S-JTSK/05 (Ferro) / Modified Krovak Parameters: Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.15281 then e = 0.081696831 e^2 = 0.006674372 Latitude of projection centre = 49°30'00"N = 0.863937979 rad Longitude of Origin = 42°30'00"E of Ferro = 0.741764932 rad Co-latitude of cone axis = 30°17'17.30311" = 0.528627763 rad Latitude of pseudo standard parallel = 78°30'00"N 1.370083463 rad Scale factor on pseudo Standard Parallel (ko) = 0.9999 False Easting = 5000000.00 m False Northing = 5000000.00 m Ordinate 1 of evaluation point Xo = 1089000.00 m Ordinate 2 of evaluation point Yo = 654000.00 m C1 = 2.946529277E-02 C2 = 2.515965696E-02 C3 = 1.193845912E-07 C4 = -4.668270147E-07 C5 = 9.233980362E-12 C6 = 1.523735715E-12 C7 = 1.696780024E-18 C8 = 4.408314235E-18 C9 = -8.331083518E-24 C10 = -3.689471323E-24 Calculated projection constants: A = 6380703.611 B = 1.000597498 gammao = 0.863239103 to = 1.003419164 n = 0.979924705 ro = 1298039.005 Forward calculation for: Latitude = 50°12'32.442"N = 0.876312568 rad Longitude = 34°30'59.179"E of Ferro = 0.602425500 rad U = 0.875596951 V = 0.139422687 S = 1.386275051 D = 0.506554627 theta = 0.496385393 r = 1194731.002 Xp = 1050538.631 Yp = 568990.995 Xr = -38461.369 Yr = -85009.005 dX = -0.077 dY = 0.088 and Southing X = 6050538.71 m Westing Y = 5568990.91 m Reverse calculation for the same Southing and Westing: Xr' = -38461.292 Yr' = -85009.093 dX' = -0.077 dY' = 0.088 Xp' = 1050538.631 Yp' = 568990.995 r' = 1194731.002 theta' = 0.496385393 D' = 0.506554627 T' = 1.386275051 U' = 0.875596951 V' = 0.139422687 lat(iteration 1) = 0.876310603 lat(iteration 2) = 0.876312562 lat(iteration 3) = 0.876312568 Latitude = 0.876312568 rad = 50°12'32.442"N Longitude = 0.294083997 rad = 16°50'59.179"E