jon.recoil.org

Source file utm.ml

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
(* UTM projection utilities — WGS84 ellipsoid, Northern hemisphere *)

let a = 6378137.0
let f = 1.0 /. 298.257223563
let b_ = a *. (1.0 -. f)
let e2 = 2.0 *. f -. f *. f
let e'2 = e2 /. (1.0 -. e2)
let n = (a -. b_) /. (a +. b_)
let k0 = 0.9996

let zone_of_lon lon =
  Float.to_int (Float.floor ((lon +. 180.0) /. 6.0)) + 1

let central_meridian zone =
  Float.of_int (zone * 6 - 183)

let wgs84_to_utm ~zone lon lat =
  let lon0 = central_meridian zone in
  let lat_r = lat *. Float.pi /. 180.0 in
  let lon_r = lon *. Float.pi /. 180.0 in
  let lon0_r = lon0 *. Float.pi /. 180.0 in
  let cos_lat = Float.cos lat_r in
  let sin_lat = Float.sin lat_r in
  let tan_lat = Float.tan lat_r in
  let n2 = n *. n in
  let n3 = n2 *. n in
  let n4 = n3 *. n in
  let n5 = n4 *. n in
  let _n6 = n5 *. n in
  (* Meridional arc using series expansion *)
  let a0 = 1.0 +. n2 /. 4.0 +. n4 /. 64.0 in
  let a2 = 1.5 *. (n -. n3 /. 8.0 +. n5 /. 64.0) in
  let a4 = 15.0 /. 16.0 *. (n2 -. n4 /. 4.0) in
  let a6 = 35.0 /. 48.0 *. (n3 -. 5.0 *. n5 /. 16.0) in
  let a8 = 315.0 /. 512.0 *. n4 in
  let a10 = 693.0 /. 1280.0 *. n5 in
  let _ = a10 in
  let big_a = a /. (1.0 +. n) *. a0 in
  let s = big_a *. (lat_r
    -. a2 *. Float.sin (2.0 *. lat_r)
    +. a4 *. Float.sin (4.0 *. lat_r)
    -. a6 *. Float.sin (6.0 *. lat_r)
    +. a8 *. Float.sin (8.0 *. lat_r)) in
  let nu = a /. Float.sqrt (1.0 -. e2 *. sin_lat *. sin_lat) in
  let t = tan_lat in
  let t2 = t *. t in
  let t4 = t2 *. t2 in
  let c = e'2 *. cos_lat *. cos_lat in
  let dl = lon_r -. lon0_r in
  let dl2 = dl *. dl in
  let dl3 = dl2 *. dl in
  let dl4 = dl3 *. dl in
  let dl5 = dl4 *. dl in
  let dl6 = dl5 *. dl in
  let easting = 500000.0 +. k0 *. nu *. cos_lat *. (dl
    +. (1.0 -. t2 +. c) *. cos_lat *. cos_lat *. dl3 /. 6.0
    +. (5.0 -. 18.0 *. t2 +. t4 +. 72.0 *. c -. 58.0 *. e'2)
       *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl5 /. 120.0) in
  let northing = k0 *. (s +. nu *. tan_lat *. (
    cos_lat *. cos_lat *. dl2 /. 2.0
    +. (5.0 -. t2 +. 9.0 *. c +. 4.0 *. c *. c) *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl4 /. 24.0
    +. (61.0 -. 58.0 *. t2 +. t4 +. 600.0 *. c -. 330.0 *. e'2)
       *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. cos_lat *. dl6 /. 720.0)) in
  (easting, northing)

let utm_to_wgs84 ~zone easting northing =
  let lon0 = central_meridian zone in
  let lon0_r = lon0 *. Float.pi /. 180.0 in
  let e1 = (1.0 -. Float.sqrt (1.0 -. e2)) /. (1.0 +. Float.sqrt (1.0 -. e2)) in
  let e1_2 = e1 *. e1 in
  let e1_3 = e1_2 *. e1 in
  let e1_4 = e1_3 *. e1 in
  let x = easting -. 500000.0 in
  let m = northing /. k0 in
  let mu = m /. (a *. (1.0 -. e2 /. 4.0 -. 3.0 *. e2 *. e2 /. 64.0 -. 5.0 *. e2 *. e2 *. e2 /. 256.0)) in
  let phi1 = mu
    +. (3.0 *. e1 /. 2.0 -. 27.0 *. e1_3 /. 32.0) *. Float.sin (2.0 *. mu)
    +. (21.0 *. e1_2 /. 16.0 -. 55.0 *. e1_4 /. 32.0) *. Float.sin (4.0 *. mu)
    +. (151.0 *. e1_3 /. 96.0) *. Float.sin (6.0 *. mu)
    +. (1097.0 *. e1_4 /. 512.0) *. Float.sin (8.0 *. mu) in
  let sin_phi1 = Float.sin phi1 in
  let cos_phi1 = Float.cos phi1 in
  let tan_phi1 = Float.tan phi1 in
  let nu1 = a /. Float.sqrt (1.0 -. e2 *. sin_phi1 *. sin_phi1) in
  let rho1 = a *. (1.0 -. e2) /. (Float.pow (1.0 -. e2 *. sin_phi1 *. sin_phi1) 1.5) in
  let t1 = tan_phi1 in
  let t1_2 = t1 *. t1 in
  let c1 = e'2 *. cos_phi1 *. cos_phi1 in
  let d = x /. (nu1 *. k0) in
  let d2 = d *. d in
  let d3 = d2 *. d in
  let d4 = d3 *. d in
  let d5 = d4 *. d in
  let d6 = d5 *. d in
  let lat = phi1 -. (nu1 *. tan_phi1 /. rho1) *. (
    d2 /. 2.0
    -. (5.0 +. 3.0 *. t1_2 +. 10.0 *. c1 -. 4.0 *. c1 *. c1 -. 9.0 *. e'2) *. d4 /. 24.0
    +. (61.0 +. 90.0 *. t1_2 +. 298.0 *. c1 +. 45.0 *. t1_2 *. t1_2 -. 252.0 *. e'2 -. 3.0 *. c1 *. c1) *. d6 /. 720.0) in
  let lon = lon0_r +. (d
    -. (1.0 +. 2.0 *. t1_2 +. c1) *. d3 /. 6.0
    +. (5.0 -. 2.0 *. c1 +. 28.0 *. t1_2 -. 3.0 *. c1 *. c1 +. 8.0 *. e'2 +. 24.0 *. t1_2 *. t1_2) *. d5 /. 120.0)
    /. cos_phi1 in
  let lat_deg = lat *. 180.0 /. Float.pi in
  let lon_deg = lon *. 180.0 /. Float.pi in
  (lon_deg, lat_deg)