commit 4acf0814771f7d903928cddf10e1583f4bfdead5 · anil.recoil.org/slop

+18

stack/lonlat/dune-project

···

       1
       1
       +
       (lang dune 3.0)

     

       2
       2
       +
       (name lonlat)

     

       3
       3
       +
       

     

       4
       4
       +
       (generate_opam_files true)

     

       5
       5
       +
       

     

       6
       6
       +
       (package

     

       7
       7
       +
        (name lonlat)

     

       8
       8
       +
        (synopsis "Geographic coordinate manipulation library for OCaml")

     

       9
       9
       +
        (description "LonLat provides types and operations for working with \

     

       10
       10
       +
                     latitude/longitude coordinates, including distance calculations, \

     

       11
       11
       +
                     bearings, bounding boxes, and coordinate parsing.")

     

       12
       12
       +
        (depends

     

       13
       13
       +
         ocaml

     

       14
       14
       +
         dune))

     

       15
       15
       +
       

     

       16
       16
       +
       (authors "LonLat Authors")

     

       17
       17
       +
       (maintainers "maintainer@example.com")

     

       18
       18
       +
       (license MIT)

+159

stack/lonlat/lib/bbox.ml

···

       1
       1
       +
       (** Bounding box module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t = {

     

       4
       4
       +
         min_lat : float;

     

       5
       5
       +
         min_lon : float;

     

       6
       6
       +
         max_lat : float;

     

       7
       7
       +
         max_lon : float;

     

       8
       8
       +
       }

     

       9
       9
       +
       

     

       10
       10
       +
       let create ~min_lat ~min_lon ~max_lat ~max_lon =

     

       11
       11
       +
         if not (Lat.is_valid min_lat) then

     

       12
       12
       +
           raise (Error.Coordinate_error (Error.Invalid_latitude min_lat));

     

       13
       13
       +
         if not (Lat.is_valid max_lat) then

     

       14
       14
       +
           raise (Error.Coordinate_error (Error.Invalid_latitude max_lat));

     

       15
       15
       +
         if min_lat > max_lat then

     

       16
       16
       +
           raise (Error.Coordinate_error

     

       17
       17
       +
             (Error.Parse_error "min_lat must be less than or equal to max_lat"));

     

       18
       18
       +
         {

     

       19
       19
       +
           min_lat;

     

       20
       20
       +
           min_lon = Lon.normalize_180 min_lon;

     

       21
       21
       +
           max_lat;

     

       22
       22
       +
           max_lon = Lon.normalize_180 max_lon;

     

       23
       23
       +
         }

     

       24
       24
       +
       

     

       25
       25
       +
       let of_coords sw ne =

     

       26
       26
       +
         create

     

       27
       27
       +
           ~min_lat:(Coord.lat_float sw)

     

       28
       28
       +
           ~min_lon:(Coord.lon_float sw)

     

       29
       29
       +
           ~max_lat:(Coord.lat_float ne)

     

       30
       30
       +
           ~max_lon:(Coord.lon_float ne)

     

       31
       31
       +
       

     

       32
       32
       +
       let of_list coords =

     

       33
       33
       +
         match coords with

     

       34
       34
       +
         | [] -> raise (Error.Coordinate_error

     

       35
       35
       +
             (Error.Parse_error "Cannot create bbox from empty list"))

     

       36
       36
       +
         | first :: rest ->

     

       37
       37
       +
             let init_lat = Coord.lat_float first in

     

       38
       38
       +
             let init_lon = Coord.lon_float first in

     

       39
       39
       +
             let min_lat, min_lon, max_lat, max_lon =

     

       40
       40
       +
               List.fold_left (fun (min_lat, min_lon, max_lat, max_lon) coord ->

     

       41
       41
       +
                 let lat = Coord.lat_float coord in

     

       42
       42
       +
                 let lon = Coord.lon_float coord in

     

       43
       43
       +
                 (min min_lat lat, min min_lon lon, max max_lat lat, max max_lon lon)

     

       44
       44
       +
               ) (init_lat, init_lon, init_lat, init_lon) rest

     

       45
       45
       +
             in

     

       46
       46
       +
             create ~min_lat ~min_lon ~max_lat ~max_lon

     

       47
       47
       +
       

     

       48
       48
       +
       let min_lat t = t.min_lat

     

       49
       49
       +
       let min_lon t = t.min_lon

     

       50
       50
       +
       let max_lat t = t.max_lat

     

       51
       51
       +
       let max_lon t = t.max_lon

     

       52
       52
       +
       

     

       53
       53
       +
       let sw_corner t = Coord.create ~lat:t.min_lat ~lon:t.min_lon

     

       54
       54
       +
       let se_corner t = Coord.create ~lat:t.min_lat ~lon:t.max_lon

     

       55
       55
       +
       let ne_corner t = Coord.create ~lat:t.max_lat ~lon:t.max_lon

     

       56
       56
       +
       let nw_corner t = Coord.create ~lat:t.max_lat ~lon:t.min_lon

     

       57
       57
       +
       

     

       58
       58
       +
       let center t =

     

       59
       59
       +
         Coord.create

     

       60
       60
       +
           ~lat:((t.min_lat +. t.max_lat) /. 2.0)

     

       61
       61
       +
           ~lon:((t.min_lon +. t.max_lon) /. 2.0)

     

       62
       62
       +
       

     

       63
       63
       +
       let width t =

     

       64
       64
       +
         if t.max_lon >= t.min_lon then

     

       65
       65
       +
           t.max_lon -. t.min_lon

     

       66
       66
       +
         else

     

       67
       67
       +
           (* Handles wraparound case *)

     

       68
       68
       +
           (t.max_lon +. 360.0) -. t.min_lon

     

       69
       69
       +
       

     

       70
       70
       +
       let height t = t.max_lat -. t.min_lat

     

       71
       71
       +
       

     

       72
       72
       +
       let contains t coord =

     

       73
       73
       +
         let lat = Coord.lat_float coord in

     

       74
       74
       +
         let lon = Coord.lon_float coord in

     

       75
       75
       +
         lat >= t.min_lat && lat <= t.max_lat &&

     

       76
       76
       +
         if t.min_lon <= t.max_lon then

     

       77
       77
       +
           (* Normal case *)

     

       78
       78
       +
           lon >= t.min_lon && lon <= t.max_lon

     

       79
       79
       +
         else

     

       80
       80
       +
           (* Crosses antimeridian *)

     

       81
       81
       +
           lon >= t.min_lon || lon <= t.max_lon

     

       82
       82
       +
       

     

       83
       83
       +
       let intersects t1 t2 =

     

       84
       84
       +
         (* Check latitude overlap *)

     

       85
       85
       +
         let lat_overlaps =

     

       86
       86
       +
           not (t1.max_lat < t2.min_lat || t2.max_lat < t1.min_lat) in

     

       87
       87
       +
       

     

       88
       88
       +
         (* Check longitude overlap *)

     

       89
       89
       +
         let lon_overlaps =

     

       90
       90
       +
           if t1.min_lon <= t1.max_lon && t2.min_lon <= t2.max_lon then

     

       91
       91
       +
             (* Neither crosses antimeridian *)

     

       92
       92
       +
             not (t1.max_lon < t2.min_lon || t2.max_lon < t1.min_lon)

     

       93
       93
       +
           else if t1.min_lon > t1.max_lon && t2.min_lon > t2.max_lon then

     

       94
       94
       +
             (* Both cross antimeridian *)

     

       95
       95
       +
             true  (* They must overlap *)

     

       96
       96
       +
           else if t1.min_lon > t1.max_lon then

     

       97
       97
       +
             (* Only t1 crosses antimeridian *)

     

       98
       98
       +
             t2.max_lon >= t1.min_lon || t2.min_lon <= t1.max_lon

     

       99
       99
       +
           else

     

       100
       100
       +
             (* Only t2 crosses antimeridian *)

     

       101
       101
       +
             t1.max_lon >= t2.min_lon || t1.min_lon <= t2.max_lon

     

       102
       102
       +
         in

     

       103
       103
       +
       

     

       104
       104
       +
         lat_overlaps && lon_overlaps

     

       105
       105
       +
       

     

       106
       106
       +
       let union t1 t2 =

     

       107
       107
       +
         create

     

       108
       108
       +
           ~min_lat:(min t1.min_lat t2.min_lat)

     

       109
       109
       +
           ~min_lon:(min t1.min_lon t2.min_lon)

     

       110
       110
       +
           ~max_lat:(max t1.max_lat t2.max_lat)

     

       111
       111
       +
           ~max_lon:(max t1.max_lon t2.max_lon)

     

       112
       112
       +
       

     

       113
       113
       +
       let intersection t1 t2 =

     

       114
       114
       +
         if not (intersects t1 t2) then

     

       115
       115
       +
           None

     

       116
       116
       +
         else

     

       117
       117
       +
           Some (create

     

       118
       118
       +
             ~min_lat:(max t1.min_lat t2.min_lat)

     

       119
       119
       +
             ~min_lon:(max t1.min_lon t2.min_lon)

     

       120
       120
       +
             ~max_lat:(min t1.max_lat t2.max_lat)

     

       121
       121
       +
             ~max_lon:(min t1.max_lon t2.max_lon))

     

       122
       122
       +
       

     

       123
       123
       +
       let expand t degrees =

     

       124
       124
       +
         create

     

       125
       125
       +
           ~min_lat:(Lat.to_float (Lat.clamp (t.min_lat -. degrees)))

     

       126
       126
       +
           ~min_lon:(t.min_lon -. degrees)

     

       127
       127
       +
           ~max_lat:(Lat.to_float (Lat.clamp (t.max_lat +. degrees)))

     

       128
       128
       +
           ~max_lon:(t.max_lon +. degrees)

     

       129
       129
       +
       

     

       130
       130
       +
       let expand_km t km =

     

       131
       131
       +
         (* Approximate degrees per km at the center latitude *)

     

       132
       132
       +
         let center_lat = (t.min_lat +. t.max_lat) /. 2.0 in

     

       133
       133
       +
         let deg_lat_per_km = 1.0 /. 111.0 in  (* ~111 km per degree latitude *)

     

       134
       134
       +
         let deg_lon_per_km = 1.0 /. (111.0 *. cos (center_lat *. Float.pi /. 180.0)) in

     

       135
       135
       +
       

     

       136
       136
       +
         create

     

       137
       137
       +
           ~min_lat:(Lat.to_float (Lat.clamp (t.min_lat -. km *. deg_lat_per_km)))

     

       138
       138
       +
           ~min_lon:(t.min_lon -. km *. deg_lon_per_km)

     

       139
       139
       +
           ~max_lat:(Lat.to_float (Lat.clamp (t.max_lat +. km *. deg_lat_per_km)))

     

       140
       140
       +
           ~max_lon:(t.max_lon +. km *. deg_lon_per_km)

     

       141
       141
       +
       

     

       142
       142
       +
       let to_string t =

     

       143
       143
       +
         Printf.sprintf "%.6f,%.6f,%.6f,%.6f" t.min_lon t.min_lat t.max_lon t.max_lat

     

       144
       144
       +
       

     

       145
       145
       +
       let of_string s =

     

       146
       146
       +
         match String.split_on_char ',' s with

     

       147
       147
       +
         | [min_lon; min_lat; max_lon; max_lat] ->

     

       148
       148
       +
             (try

     

       149
       149
       +
               create

     

       150
       150
       +
                 ~min_lat:(float_of_string (String.trim min_lat))

     

       151
       151
       +
                 ~min_lon:(float_of_string (String.trim min_lon))

     

       152
       152
       +
                 ~max_lat:(float_of_string (String.trim max_lat))

     

       153
       153
       +
                 ~max_lon:(float_of_string (String.trim max_lon))

     

       154
       154
       +
             with Failure _ ->

     

       155
       155
       +
               raise (Error.Coordinate_error

     

       156
       156
       +
                 (Error.Parse_error "Invalid bbox format: expected numeric values")))

     

       157
       157
       +
         | _ ->

     

       158
       158
       +
             raise (Error.Coordinate_error

     

       159
       159
       +
               (Error.Parse_error "Invalid bbox format: expected min_lon,min_lat,max_lon,max_lat"))

+53

stack/lonlat/lib/bbox.mli

···

       1
       1
       +
       (** Bounding box module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t

     

       4
       4
       +
       

     

       5
       5
       +
       (** {2 Construction} *)

     

       6
       6
       +
       

     

       7
       7
       +
       val create : min_lat:float -> min_lon:float -> max_lat:float -> max_lon:float -> t

     

       8
       8
       +
       

     

       9
       9
       +
       val of_coords : Coord.t -> Coord.t -> t

     

       10
       10
       +
       (** From SW and NE corners *)

     

       11
       11
       +
       

     

       12
       12
       +
       val of_list : Coord.t list -> t

     

       13
       13
       +
       (** Compute bbox containing all points *)

     

       14
       14
       +
       

     

       15
       15
       +
       (** {2 Accessors} *)

     

       16
       16
       +
       

     

       17
       17
       +
       val min_lat : t -> float

     

       18
       18
       +
       val min_lon : t -> float

     

       19
       19
       +
       val max_lat : t -> float

     

       20
       20
       +
       val max_lon : t -> float

     

       21
       21
       +
       

     

       22
       22
       +
       val sw_corner : t -> Coord.t

     

       23
       23
       +
       val se_corner : t -> Coord.t

     

       24
       24
       +
       val ne_corner : t -> Coord.t

     

       25
       25
       +
       val nw_corner : t -> Coord.t

     

       26
       26
       +
       val center : t -> Coord.t

     

       27
       27
       +
       

     

       28
       28
       +
       val width : t -> float

     

       29
       29
       +
       (** Width in degrees *)

     

       30
       30
       +
       

     

       31
       31
       +
       val height : t -> float

     

       32
       32
       +
       (** Height in degrees *)

     

       33
       33
       +
       

     

       34
       34
       +
       (** {2 Operations} *)

     

       35
       35
       +
       

     

       36
       36
       +
       val contains : t -> Coord.t -> bool

     

       37
       37
       +
       val intersects : t -> t -> bool

     

       38
       38
       +
       val union : t -> t -> t

     

       39
       39
       +
       val intersection : t -> t -> t option

     

       40
       40
       +
       

     

       41
       41
       +
       val expand : t -> float -> t

     

       42
       42
       +
       (** Expand by degrees on all sides *)

     

       43
       43
       +
       

     

       44
       44
       +
       val expand_km : t -> float -> t

     

       45
       45
       +
       (** Expand by approximate km on all sides *)

     

       46
       46
       +
       

     

       47
       47
       +
       (** {2 Formatting} *)

     

       48
       48
       +
       

     

       49
       49
       +
       val to_string : t -> string

     

       50
       50
       +
       (** Format as "min_lon,min_lat,max_lon,max_lat" *)

     

       51
       51
       +
       

     

       52
       52
       +
       val of_string : string -> t

     

       53
       53
       +
       (** Parse from "min_lon,min_lat,max_lon,max_lat" *)

+195

stack/lonlat/lib/coord.ml

···

       1
       1
       +
       (** Coordinate pair module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t = {

     

       4
       4
       +
         lat : Lat.t;

     

       5
       5
       +
         lon : Lon.t;

     

       6
       6
       +
       }

     

       7
       7
       +
       

     

       8
       8
       +
       let create ~lat ~lon = {

     

       9
       9
       +
         lat = Lat.create lat;

     

       10
       10
       +
         lon = Lon.create lon;

     

       11
       11
       +
       }

     

       12
       12
       +
       

     

       13
       13
       +
       let of_lat_lon lat lon = { lat; lon }

     

       14
       14
       +
       

     

       15
       15
       +
       let lat t = t.lat

     

       16
       16
       +
       let lon t = t.lon

     

       17
       17
       +
       let lat_float t = Lat.to_float t.lat

     

       18
       18
       +
       let lon_float t = Lon.to_float t.lon

     

       19
       19
       +
       

     

       20
       20
       +
       (** Constants *)

     

       21
       21
       +
       let earth_radius_km = 6371.0

     

       22
       22
       +
       let deg_to_rad = Float.pi /. 180.0

     

       23
       23
       +
       let rad_to_deg = 180.0 /. Float.pi

     

       24
       24
       +
       

     

       25
       25
       +
       (** Distance calculations *)

     

       26
       26
       +
       

     

       27
       27
       +
       let distance_haversine p1 p2 =

     

       28
       28
       +
         let lat1 = lat_float p1 *. deg_to_rad in

     

       29
       29
       +
         let lat2 = lat_float p2 *. deg_to_rad in

     

       30
       30
       +
         let lon1 = lon_float p1 *. deg_to_rad in

     

       31
       31
       +
         let lon2 = lon_float p2 *. deg_to_rad in

     

       32
       32
       +
       

     

       33
       33
       +
         let dlat = lat2 -. lat1 in

     

       34
       34
       +
         let dlon = lon2 -. lon1 in

     

       35
       35
       +
       

     

       36
       36
       +
         let a =

     

       37
       37
       +
           (sin (dlat /. 2.0) ** 2.0) +.

     

       38
       38
       +
           cos lat1 *. cos lat2 *. (sin (dlon /. 2.0) ** 2.0) in

     

       39
       39
       +
       

     

       40
       40
       +
         let c = 2.0 *. atan2 (sqrt a) (sqrt (1.0 -. a)) in

     

       41
       41
       +
         earth_radius_km *. c

     

       42
       42
       +
       

     

       43
       43
       +
       let distance_vincenty p1 p2 =

     

       44
       44
       +
         (* WGS84 ellipsoid constants *)

     

       45
       45
       +
         let a = 6378.137 in  (* semi-major axis in km *)

     

       46
       46
       +
         let b = 6356.752 in  (* semi-minor axis in km *)

     

       47
       47
       +
         let f = 1.0 /. 298.257223563 in  (* flattening *)

     

       48
       48
       +
       

     

       49
       49
       +
         let lat1 = lat_float p1 *. deg_to_rad in

     

       50
       50
       +
         let lat2 = lat_float p2 *. deg_to_rad in

     

       51
       51
       +
         let lon1 = lon_float p1 *. deg_to_rad in

     

       52
       52
       +
         let lon2 = lon_float p2 *. deg_to_rad in

     

       53
       53
       +
       

     

       54
       54
       +
         let l = lon2 -. lon1 in

     

       55
       55
       +
         let u1 = atan ((1.0 -. f) *. tan lat1) in

     

       56
       56
       +
         let u2 = atan ((1.0 -. f) *. tan lat2) in

     

       57
       57
       +
         let sin_u1 = sin u1 in

     

       58
       58
       +
         let cos_u1 = cos u1 in

     

       59
       59
       +
         let sin_u2 = sin u2 in

     

       60
       60
       +
         let cos_u2 = cos u2 in

     

       61
       61
       +
       

     

       62
       62
       +
         let rec iterate lambda _prev_lambda iter =

     

       63
       63
       +
           if iter > 100 then

     

       64
       64
       +
             (* Failed to converge, fall back to haversine *)

     

       65
       65
       +
             distance_haversine p1 p2

     

       66
       66
       +
           else

     

       67
       67
       +
             let sin_lambda = sin lambda in

     

       68
       68
       +
             let cos_lambda = cos lambda in

     

       69
       69
       +
             let sin_sigma =

     

       70
       70
       +
               sqrt ((cos_u2 *. sin_lambda) ** 2.0 +.

     

       71
       71
       +
                     (cos_u1 *. sin_u2 -. sin_u1 *. cos_u2 *. cos_lambda) ** 2.0) in

     

       72
       72
       +
       

     

       73
       73
       +
             if sin_sigma = 0.0 then 0.0  (* coincident points *)

     

       74
       74
       +
             else

     

       75
       75
       +
               let cos_sigma = sin_u1 *. sin_u2 +. cos_u1 *. cos_u2 *. cos_lambda in

     

       76
       76
       +
               let sigma = atan2 sin_sigma cos_sigma in

     

       77
       77
       +
               let sin_alpha = cos_u1 *. cos_u2 *. sin_lambda /. sin_sigma in

     

       78
       78
       +
               let cos_sq_alpha = 1.0 -. sin_alpha ** 2.0 in

     

       79
       79
       +
               let cos_2sigma_m =

     

       80
       80
       +
                 if cos_sq_alpha = 0.0 then 0.0  (* equatorial line *)

     

       81
       81
       +
                 else cos_sigma -. 2.0 *. sin_u1 *. sin_u2 /. cos_sq_alpha in

     

       82
       82
       +
       

     

       83
       83
       +
               let c = f /. 16.0 *. cos_sq_alpha *. (4.0 +. f *. (4.0 -. 3.0 *. cos_sq_alpha)) in

     

       84
       84
       +
               let lambda' = l +. (1.0 -. c) *. f *. sin_alpha *.

     

       85
       85
       +
                             (sigma +. c *. sin_sigma *.

     

       86
       86
       +
                              (cos_2sigma_m +. c *. cos_sigma *.

     

       87
       87
       +
                               (-1.0 +. 2.0 *. cos_2sigma_m ** 2.0))) in

     

       88
       88
       +
       

     

       89
       89
       +
               if abs_float (lambda' -. lambda) < 1e-12 then

     

       90
       90
       +
                 (* Converged *)

     

       91
       91
       +
                 let u_sq = cos_sq_alpha *. (a ** 2.0 -. b ** 2.0) /. (b ** 2.0) in

     

       92
       92
       +
                 let cap_a = 1.0 +. u_sq /. 16384.0 *.

     

       93
       93
       +
                             (4096.0 +. u_sq *. (-768.0 +. u_sq *. (320.0 -. 175.0 *. u_sq))) in

     

       94
       94
       +
                 let cap_b = u_sq /. 1024.0 *. (256.0 +. u_sq *. (-128.0 +. u_sq *. (74.0 -. 47.0 *. u_sq))) in

     

       95
       95
       +
                 let delta_sigma = cap_b *. sin_sigma *.

     

       96
       96
       +
                                   (cos_2sigma_m +. cap_b /. 4.0 *.

     

       97
       97
       +
                                    (cos_sigma *. (-1.0 +. 2.0 *. cos_2sigma_m ** 2.0) -.

     

       98
       98
       +
                                     cap_b /. 6.0 *. cos_2sigma_m *.

     

       99
       99
       +
                                     (-3.0 +. 4.0 *. sin_sigma ** 2.0) *.

     

       100
       100
       +
                                     (-3.0 +. 4.0 *. cos_2sigma_m ** 2.0))) in

     

       101
       101
       +
       

     

       102
       102
       +
                 b *. cap_a *. (sigma -. delta_sigma)

     

       103
       103
       +
               else

     

       104
       104
       +
                 iterate lambda' lambda (iter + 1)

     

       105
       105
       +
         in

     

       106
       106
       +
       

     

       107
       107
       +
         iterate l 0.0 0

     

       108
       108
       +
       

     

       109
       109
       +
       (** Bearing calculations *)

     

       110
       110
       +
       

     

       111
       111
       +
       let bearing p1 p2 =

     

       112
       112
       +
         let lat1 = lat_float p1 *. deg_to_rad in

     

       113
       113
       +
         let lat2 = lat_float p2 *. deg_to_rad in

     

       114
       114
       +
         let lon1 = lon_float p1 *. deg_to_rad in

     

       115
       115
       +
         let lon2 = lon_float p2 *. deg_to_rad in

     

       116
       116
       +
       

     

       117
       117
       +
         let dlon = lon2 -. lon1 in

     

       118
       118
       +
         let y = sin dlon *. cos lat2 in

     

       119
       119
       +
         let x = cos lat1 *. sin lat2 -. sin lat1 *. cos lat2 *. cos dlon in

     

       120
       120
       +
       

     

       121
       121
       +
         let bearing = atan2 y x *. rad_to_deg in

     

       122
       122
       +
         mod_float (bearing +. 360.0) 360.0

     

       123
       123
       +
       

     

       124
       124
       +
       let bearing_final p1 p2 =

     

       125
       125
       +
         (* Final bearing is the reverse of the initial bearing from p2 to p1 *)

     

       126
       126
       +
         let reverse_bearing = bearing p2 p1 in

     

       127
       127
       +
         mod_float (reverse_bearing +. 180.0) 360.0

     

       128
       128
       +
       

     

       129
       129
       +
       (** Geographic operations *)

     

       130
       130
       +
       

     

       131
       131
       +
       let offset t ~heading ~distance =

     

       132
       132
       +
         let lat1 = lat_float t *. deg_to_rad in

     

       133
       133
       +
         let lon1 = lon_float t *. deg_to_rad in

     

       134
       134
       +
         let brng = heading *. deg_to_rad in

     

       135
       135
       +
         let d = distance /. earth_radius_km in

     

       136
       136
       +
       

     

       137
       137
       +
         let lat2 = asin (sin lat1 *. cos d +. cos lat1 *. sin d *. cos brng) in

     

       138
       138
       +
         let lon2 = lon1 +. atan2 (sin brng *. sin d *. cos lat1)

     

       139
       139
       +
                                  (cos d -. sin lat1 *. sin lat2) in

     

       140
       140
       +
       

     

       141
       141
       +
         create ~lat:(lat2 *. rad_to_deg) ~lon:(lon2 *. rad_to_deg)

     

       142
       142
       +
       

     

       143
       143
       +
       let midpoint p1 p2 =

     

       144
       144
       +
         let lat1 = lat_float p1 *. deg_to_rad in

     

       145
       145
       +
         let lat2 = lat_float p2 *. deg_to_rad in

     

       146
       146
       +
         let lon1 = lon_float p1 *. deg_to_rad in

     

       147
       147
       +
         let lon2 = lon_float p2 *. deg_to_rad in

     

       148
       148
       +
       

     

       149
       149
       +
         let dlon = lon2 -. lon1 in

     

       150
       150
       +
         let bx = cos lat2 *. cos dlon in

     

       151
       151
       +
         let by = cos lat2 *. sin dlon in

     

       152
       152
       +
       

     

       153
       153
       +
         let lat3 = atan2 (sin lat1 +. sin lat2)

     

       154
       154
       +
                          (sqrt ((cos lat1 +. bx) ** 2.0 +. by ** 2.0)) in

     

       155
       155
       +
         let lon3 = lon1 +. atan2 by (cos lat1 +. bx) in

     

       156
       156
       +
       

     

       157
       157
       +
         create ~lat:(lat3 *. rad_to_deg) ~lon:(lon3 *. rad_to_deg)

     

       158
       158
       +
       

     

       159
       159
       +
       let interpolate p1 p2 fraction =

     

       160
       160
       +
         if fraction <= 0.0 then p1

     

       161
       161
       +
         else if fraction >= 1.0 then p2

     

       162
       162
       +
         else

     

       163
       163
       +
           let lat1 = lat_float p1 in

     

       164
       164
       +
           let lat2 = lat_float p2 in

     

       165
       165
       +
           let lon1 = lon_float p1 in

     

       166
       166
       +
           let lon2 = lon_float p2 in

     

       167
       167
       +
       

     

       168
       168
       +
           (* Simple linear interpolation in degrees *)

     

       169
       169
       +
           let lat = lat1 +. (lat2 -. lat1) *. fraction in

     

       170
       170
       +
           let lon = lon1 +. (lon2 -. lon1) *. fraction in

     

       171
       171
       +
       

     

       172
       172
       +
           create ~lat ~lon

     

       173
       173
       +
       

     

       174
       174
       +
       (** Comparison *)

     

       175
       175
       +
       

     

       176
       176
       +
       let equal p1 p2 =

     

       177
       177
       +
         Lat.equal p1.lat p2.lat && Lon.equal p1.lon p2.lon

     

       178
       178
       +
       

     

       179
       179
       +
       let ( = ) = equal

     

       180
       180
       +
       

     

       181
       181
       +
       let almost_equal ?(epsilon = 0.000001) p1 p2 =

     

       182
       182
       +
         distance_haversine p1 p2 < epsilon

     

       183
       183
       +
       

     

       184
       184
       +
       (** Formatting *)

     

       185
       185
       +
       

     

       186
       186
       +
       let to_string t =

     

       187
       187
       +
         Printf.sprintf "%s, %s" (Lat.to_string t.lat) (Lon.to_string t.lon)

     

       188
       188
       +
       

     

       189
       189
       +
       let to_dms_string t =

     

       190
       190
       +
         Printf.sprintf "%s, %s" (Lat.to_dms_string t.lat) (Lon.to_dms_string t.lon)

     

       191
       191
       +
       

     

       192
       192
       +
       let format ?(lat_precision = 6) ?(lon_precision = 6) t =

     

       193
       193
       +
         Printf.sprintf "%s, %s"

     

       194
       194
       +
           (Lat.format ~precision:lat_precision t.lat)

     

       195
       195
       +
           (Lon.format ~precision:lon_precision t.lon)

+59

stack/lonlat/lib/coord.mli

···

       1
       1
       +
       (** Coordinate pair module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t

     

       4
       4
       +
       

     

       5
       5
       +
       (** {2 Construction} *)

     

       6
       6
       +
       

     

       7
       7
       +
       val create : lat:float -> lon:float -> t

     

       8
       8
       +
       

     

       9
       9
       +
       val of_lat_lon : Lat.t -> Lon.t -> t

     

       10
       10
       +
       

     

       11
       11
       +
       (** {2 Accessors} *)

     

       12
       12
       +
       

     

       13
       13
       +
       val lat : t -> Lat.t

     

       14
       14
       +
       val lon : t -> Lon.t

     

       15
       15
       +
       val lat_float : t -> float

     

       16
       16
       +
       val lon_float : t -> float

     

       17
       17
       +
       

     

       18
       18
       +
       (** {2 Distance calculations} *)

     

       19
       19
       +
       

     

       20
       20
       +
       val distance_haversine : t -> t -> float

     

       21
       21
       +
       (** Great circle distance in km using Earth radius 6371.0 *)

     

       22
       22
       +
       

     

       23
       23
       +
       val distance_vincenty : t -> t -> float

     

       24
       24
       +
       (** WGS84 ellipsoid distance in km *)

     

       25
       25
       +
       

     

       26
       26
       +
       (** {2 Bearing calculations} *)

     

       27
       27
       +
       

     

       28
       28
       +
       val bearing : t -> t -> float

     

       29
       29
       +
       (** Initial bearing in degrees [0, 360) *)

     

       30
       30
       +
       

     

       31
       31
       +
       val bearing_final : t -> t -> float

     

       32
       32
       +
       (** Final bearing in degrees *)

     

       33
       33
       +
       

     

       34
       34
       +
       (** {2 Geographic operations} *)

     

       35
       35
       +
       

     

       36
       36
       +
       val offset : t -> heading:float -> distance:float -> t

     

       37
       37
       +
       (** New point from heading (degrees) and distance (km) *)

     

       38
       38
       +
       

     

       39
       39
       +
       val midpoint : t -> t -> t

     

       40
       40
       +
       

     

       41
       41
       +
       val interpolate : t -> t -> float -> t

     

       42
       42
       +
       (** Linear interpolation: 0.0 = first, 1.0 = second *)

     

       43
       43
       +
       

     

       44
       44
       +
       (** {2 Comparison} *)

     

       45
       45
       +
       

     

       46
       46
       +
       val equal : t -> t -> bool

     

       47
       47
       +
       val ( = ) : t -> t -> bool

     

       48
       48
       +
       

     

       49
       49
       +
       val almost_equal : ?epsilon:float -> t -> t -> bool

     

       50
       50
       +
       (** Default epsilon = 0.000001 km *)

     

       51
       51
       +
       

     

       52
       52
       +
       (** {2 Formatting} *)

     

       53
       53
       +
       

     

       54
       54
       +
       val to_string : t -> string

     

       55
       55
       +
       (** "lat, lon" in decimal degrees *)

     

       56
       56
       +
       

     

       57
       57
       +
       val to_dms_string : t -> string

     

       58
       58
       +
       

     

       59
       59
       +
       val format : ?lat_precision:int -> ?lon_precision:int -> t -> string

stack/lonlat/lib/dune

···

       1
       1
       +
       (library

     

       2
       2
       +
        (public_name lonlat)

     

       3
       3
       +
        (name lonlat)

     

       4
       4
       +
        (modules error lat lon coord bbox vector parser lonlat)

     

       5
       5
       +
        (libraries str))

+33

stack/lonlat/lib/error.ml

···

       1
       1
       +
       (** Error handling for geographic coordinates *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t =

     

       4
       4
       +
         | Invalid_latitude of float

     

       5
       5
       +
         | Invalid_longitude of float

     

       6
       6
       +
         | Parse_error of string

     

       7
       7
       +
         | Invalid_format of string

     

       8
       8
       +
         | Invalid_hemisphere of string

     

       9
       9
       +
       

     

       10
       10
       +
       exception Coordinate_error of t

     

       11
       11
       +
       

     

       12
       12
       +
       let to_string = function

     

       13
       13
       +
         | Invalid_latitude f -> Printf.sprintf "Invalid latitude: %.6f (must be between -90 and 90)" f

     

       14
       14
       +
         | Invalid_longitude f -> Printf.sprintf "Invalid longitude: %.6f" f

     

       15
       15
       +
         | Parse_error s -> Printf.sprintf "Parse error: %s" s

     

       16
       16
       +
         | Invalid_format s -> Printf.sprintf "Invalid format: %s" s

     

       17
       17
       +
         | Invalid_hemisphere s -> Printf.sprintf "Invalid hemisphere: %s" s

     

       18
       18
       +
       

     

       19
       19
       +
       let catch f x =

     

       20
       20
       +
         try Ok (f x)

     

       21
       21
       +
         with Coordinate_error e -> Error e

     

       22
       22
       +
       

     

       23
       23
       +
       let catch_exn f x =

     

       24
       24
       +
         try Some (f x)

     

       25
       25
       +
         with Coordinate_error _ -> None

     

       26
       26
       +
       

     

       27
       27
       +
       let unwrap = function

     

       28
       28
       +
         | Ok x -> x

     

       29
       29
       +
         | Error e -> raise (Coordinate_error e)

     

       30
       30
       +
       

     

       31
       31
       +
       let unwrap_or default = function

     

       32
       32
       +
         | Ok x -> x

     

       33
       33
       +
         | Error _ -> default

+18

stack/lonlat/lib/error.mli

···

       1
       1
       +
       (** Error handling for geographic coordinates *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t =

     

       4
       4
       +
         | Invalid_latitude of float

     

       5
       5
       +
         | Invalid_longitude of float

     

       6
       6
       +
         | Parse_error of string

     

       7
       7
       +
         | Invalid_format of string

     

       8
       8
       +
         | Invalid_hemisphere of string

     

       9
       9
       +
       

     

       10
       10
       +
       exception Coordinate_error of t

     

       11
       11
       +
       

     

       12
       12
       +
       val to_string : t -> string

     

       13
       13
       +
       

     

       14
       14
       +
       (** Error handling combinators *)

     

       15
       15
       +
       val catch : ('a -> 'b) -> 'a -> ('b, t) result

     

       16
       16
       +
       val catch_exn : ('a -> 'b) -> 'a -> 'b option

     

       17
       17
       +
       val unwrap : ('a, t) result -> 'a

     

       18
       18
       +
       val unwrap_or : 'a -> ('a, t) result -> 'a

+66

stack/lonlat/lib/lat.ml

···

       1
       1
       +
       (** Latitude module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t = float

     

       4
       4
       +
       

     

       5
       5
       +
       let is_valid f = f >= -90.0 && f <= 90.0

     

       6
       6
       +
       

     

       7
       7
       +
       let create f =

     

       8
       8
       +
         if not (is_valid f) then

     

       9
       9
       +
           raise (Error.Coordinate_error (Error.Invalid_latitude f))

     

       10
       10
       +
         else f

     

       11
       11
       +
       

     

       12
       12
       +
       let unsafe_create f = f

     

       13
       13
       +
       

     

       14
       14
       +
       let clamp f =

     

       15
       15
       +
         if f < -90.0 then -90.0

     

       16
       16
       +
         else if f > 90.0 then 90.0

     

       17
       17
       +
         else f

     

       18
       18
       +
       

     

       19
       19
       +
       let of_dms ~degree ~minute ~second =

     

       20
       20
       +
         let sign = if degree < 0 || minute < 0 || second < 0.0 then -1.0 else 1.0 in

     

       21
       21
       +
         let d = abs degree in

     

       22
       22
       +
         let m = abs minute in

     

       23
       23
       +
         let s = abs_float second in

     

       24
       24
       +
         let decimal = float_of_int d +. (float_of_int m /. 60.0) +. (s /. 3600.0) in

     

       25
       25
       +
         create (sign *. decimal)

     

       26
       26
       +
       

     

       27
       27
       +
       let to_float t = t

     

       28
       28
       +
       

     

       29
       29
       +
       let degree t = int_of_float t

     

       30
       30
       +
       

     

       31
       31
       +
       let decimal_minute t =

     

       32
       32
       +
         let d = abs_float t in

     

       33
       33
       +
         (d -. floor d) *. 60.0

     

       34
       34
       +
       

     

       35
       35
       +
       let minute t = int_of_float (decimal_minute t)

     

       36
       36
       +
       

     

       37
       37
       +
       let second t =

     

       38
       38
       +
         let dm = decimal_minute t in

     

       39
       39
       +
         (dm -. floor dm) *. 60.0

     

       40
       40
       +
       

     

       41
       41
       +
       let hemisphere t = if t < 0.0 then `S else `N

     

       42
       42
       +
       

     

       43
       43
       +
       let add t delta = create (t +. delta)

     

       44
       44
       +
       let sub t delta = create (t -. delta)

     

       45
       45
       +
       let neg t = create (-.t)

     

       46
       46
       +
       let abs t = abs_float t

     

       47
       47
       +
       

     

       48
       48
       +
       let equal a b = a = b

     

       49
       49
       +
       let compare = Float.compare

     

       50
       50
       +
       let ( = ) = equal

     

       51
       51
       +
       let ( < ) a b = compare a b < 0

     

       52
       52
       +
       let ( > ) a b = compare a b > 0

     

       53
       53
       +
       let ( <= ) a b = compare a b <= 0

     

       54
       54
       +
       let ( >= ) a b = compare a b >= 0

     

       55
       55
       +
       

     

       56
       56
       +
       let to_string t = Printf.sprintf "%.6f" t

     

       57
       57
       +
       

     

       58
       58
       +
       let format ?(precision = 6) t =

     

       59
       59
       +
         Printf.sprintf "%.*f" precision t

     

       60
       60
       +
       

     

       61
       61
       +
       let to_dms_string t =

     

       62
       62
       +
         let h = match hemisphere t with `N -> "N" | `S -> "S" in

     

       63
       63
       +
         let d = abs_float (float_of_int (degree t)) |> int_of_float in

     

       64
       64
       +
         let m = minute (abs_float t) in

     

       65
       65
       +
         let s = second (abs_float t) in

     

       66
       66
       +
         Printf.sprintf "%d°%02d'%05.2f\"%s" d m s h

+57

stack/lonlat/lib/lat.mli

···

       1
       1
       +
       (** Latitude module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t

     

       4
       4
       +
       

     

       5
       5
       +
       (** {2 Construction} *)

     

       6
       6
       +
       

     

       7
       7
       +
       val create : float -> t

     

       8
       8
       +
       (** Raises Error.Coordinate_error if outside [-90, 90] *)

     

       9
       9
       +
       

     

       10
       10
       +
       val of_dms : degree:int -> minute:int -> second:float -> t

     

       11
       11
       +
       

     

       12
       12
       +
       val unsafe_create : float -> t

     

       13
       13
       +
       (** No validation - use with caution *)

     

       14
       14
       +
       

     

       15
       15
       +
       val clamp : float -> t

     

       16
       16
       +
       (** Clamp to valid range [-90, 90] *)

     

       17
       17
       +
       

     

       18
       18
       +
       (** {2 Accessors} *)

     

       19
       19
       +
       

     

       20
       20
       +
       val to_float : t -> float

     

       21
       21
       +
       val degree : t -> int

     

       22
       22
       +
       val minute : t -> int

     

       23
       23
       +
       val second : t -> float

     

       24
       24
       +
       val decimal_minute : t -> float

     

       25
       25
       +
       val hemisphere : t -> [`N | `S]

     

       26
       26
       +
       

     

       27
       27
       +
       (** {2 Operations} *)

     

       28
       28
       +
       

     

       29
       29
       +
       val add : t -> float -> t

     

       30
       30
       +
       val sub : t -> float -> t

     

       31
       31
       +
       val neg : t -> t

     

       32
       32
       +
       val abs : t -> float

     

       33
       33
       +
       

     

       34
       34
       +
       (** {2 Comparison} *)

     

       35
       35
       +
       

     

       36
       36
       +
       val equal : t -> t -> bool

     

       37
       37
       +
       val compare : t -> t -> int

     

       38
       38
       +
       val ( = ) : t -> t -> bool

     

       39
       39
       +
       val ( < ) : t -> t -> bool

     

       40
       40
       +
       val ( > ) : t -> t -> bool

     

       41
       41
       +
       val ( <= ) : t -> t -> bool

     

       42
       42
       +
       val ( >= ) : t -> t -> bool

     

       43
       43
       +
       

     

       44
       44
       +
       (** {2 Validation} *)

     

       45
       45
       +
       

     

       46
       46
       +
       val is_valid : float -> bool

     

       47
       47
       +
       

     

       48
       48
       +
       (** {2 Formatting} *)

     

       49
       49
       +
       

     

       50
       50
       +
       val to_string : t -> string

     

       51
       51
       +
       (** Decimal degrees *)

     

       52
       52
       +
       

     

       53
       53
       +
       val format : ?precision:int -> t -> string

     

       54
       54
       +
       (** With specified precision (default 6) *)

     

       55
       55
       +
       

     

       56
       56
       +
       val to_dms_string : t -> string

     

       57
       57
       +
       (** Format as "DD°MM'SS.S\"N" *)

+71

stack/lonlat/lib/lon.ml

···

       1
       1
       +
       (** Longitude module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t = float

     

       4
       4
       +
       

     

       5
       5
       +
       let normalize_180 f =

     

       6
       6
       +
         let x = mod_float (f +. 180.0) 360.0 in

     

       7
       7
       +
         if x < 0.0 then x +. 360.0 -. 180.0

     

       8
       8
       +
         else x -. 180.0

     

       9
       9
       +
       

     

       10
       10
       +
       let normalize_360 f =

     

       11
       11
       +
         let x = mod_float f 360.0 in

     

       12
       12
       +
         if x < 0.0 then x +. 360.0 else x

     

       13
       13
       +
       

     

       14
       14
       +
       let create f = normalize_180 f

     

       15
       15
       +
       

     

       16
       16
       +
       let unsafe_create f = f

     

       17
       17
       +
       

     

       18
       18
       +
       let of_dms ~degree ~minute ~second =

     

       19
       19
       +
         let sign = if degree < 0 || minute < 0 || second < 0.0 then -1.0 else 1.0 in

     

       20
       20
       +
         let d = abs degree in

     

       21
       21
       +
         let m = abs minute in

     

       22
       22
       +
         let s = abs_float second in

     

       23
       23
       +
         let decimal = float_of_int d +. (float_of_int m /. 60.0) +. (s /. 3600.0) in

     

       24
       24
       +
         create (sign *. decimal)

     

       25
       25
       +
       

     

       26
       26
       +
       let to_float t = t

     

       27
       27
       +
       

     

       28
       28
       +
       let to_range_180 t = normalize_180 t

     

       29
       29
       +
       let to_range_360 t = normalize_360 t

     

       30
       30
       +
       

     

       31
       31
       +
       let degree t = int_of_float t

     

       32
       32
       +
       

     

       33
       33
       +
       let decimal_minute t =

     

       34
       34
       +
         let d = abs_float t in

     

       35
       35
       +
         (d -. floor d) *. 60.0

     

       36
       36
       +
       

     

       37
       37
       +
       let minute t = int_of_float (decimal_minute t)

     

       38
       38
       +
       

     

       39
       39
       +
       let second t =

     

       40
       40
       +
         let dm = decimal_minute t in

     

       41
       41
       +
         (dm -. floor dm) *. 60.0

     

       42
       42
       +
       

     

       43
       43
       +
       let hemisphere t = if t < 0.0 then `W else `E

     

       44
       44
       +
       

     

       45
       45
       +
       let add t delta = create (t +. delta)

     

       46
       46
       +
       let sub t delta = create (t -. delta)

     

       47
       47
       +
       let neg t = create (-.t)

     

       48
       48
       +
       let abs t = abs_float t

     

       49
       49
       +
       

     

       50
       50
       +
       let equal a b =

     

       51
       51
       +
         (* Longitudes wrap around, so normalize before comparing *)

     

       52
       52
       +
         normalize_180 a = normalize_180 b

     

       53
       53
       +
       

     

       54
       54
       +
       let compare a b = Float.compare (normalize_180 a) (normalize_180 b)

     

       55
       55
       +
       let ( = ) = equal

     

       56
       56
       +
       let ( < ) a b = compare a b < 0

     

       57
       57
       +
       let ( > ) a b = compare a b > 0

     

       58
       58
       +
       let ( <= ) a b = compare a b <= 0

     

       59
       59
       +
       let ( >= ) a b = compare a b >= 0

     

       60
       60
       +
       

     

       61
       61
       +
       let to_string t = Printf.sprintf "%.6f" t

     

       62
       62
       +
       

     

       63
       63
       +
       let format ?(precision = 6) t =

     

       64
       64
       +
         Printf.sprintf "%.*f" precision t

     

       65
       65
       +
       

     

       66
       66
       +
       let to_dms_string t =

     

       67
       67
       +
         let h = match hemisphere t with `E -> "E" | `W -> "W" in

     

       68
       68
       +
         let d = abs_float (float_of_int (degree t)) |> int_of_float in

     

       69
       69
       +
         let m = minute (abs_float t) in

     

       70
       70
       +
         let s = second (abs_float t) in

     

       71
       71
       +
         Printf.sprintf "%03d°%02d'%05.2f\"%s" d m s h

+56

stack/lonlat/lib/lon.mli

···

       1
       1
       +
       (** Longitude module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t

     

       4
       4
       +
       

     

       5
       5
       +
       (** {2 Construction} *)

     

       6
       6
       +
       

     

       7
       7
       +
       val create : float -> t

     

       8
       8
       +
       (** Normalizes to [-180, 180] *)

     

       9
       9
       +
       

     

       10
       10
       +
       val of_dms : degree:int -> minute:int -> second:float -> t

     

       11
       11
       +
       

     

       12
       12
       +
       val unsafe_create : float -> t

     

       13
       13
       +
       

     

       14
       14
       +
       (** {2 Accessors} *)

     

       15
       15
       +
       

     

       16
       16
       +
       val to_float : t -> float

     

       17
       17
       +
       val degree : t -> int

     

       18
       18
       +
       val minute : t -> int

     

       19
       19
       +
       val second : t -> float

     

       20
       20
       +
       val decimal_minute : t -> float

     

       21
       21
       +
       val hemisphere : t -> [`E | `W]

     

       22
       22
       +
       

     

       23
       23
       +
       (** {2 Range operations} *)

     

       24
       24
       +
       

     

       25
       25
       +
       val normalize_180 : float -> float

     

       26
       26
       +
       (** Normalize any value to [-180, 180] *)

     

       27
       27
       +
       

     

       28
       28
       +
       val normalize_360 : float -> float

     

       29
       29
       +
       (** Normalize any value to [0, 360] *)

     

       30
       30
       +
       

     

       31
       31
       +
       val to_range_180 : t -> float

     

       32
       32
       +
       val to_range_360 : t -> float

     

       33
       33
       +
       

     

       34
       34
       +
       (** {2 Operations} *)

     

       35
       35
       +
       

     

       36
       36
       +
       val add : t -> float -> t

     

       37
       37
       +
       val sub : t -> float -> t

     

       38
       38
       +
       val neg : t -> t

     

       39
       39
       +
       val abs : t -> float

     

       40
       40
       +
       

     

       41
       41
       +
       (** {2 Comparison} *)

     

       42
       42
       +
       

     

       43
       43
       +
       val equal : t -> t -> bool

     

       44
       44
       +
       val compare : t -> t -> int

     

       45
       45
       +
       val ( = ) : t -> t -> bool

     

       46
       46
       +
       val ( < ) : t -> t -> bool

     

       47
       47
       +
       val ( > ) : t -> t -> bool

     

       48
       48
       +
       val ( <= ) : t -> t -> bool

     

       49
       49
       +
       val ( >= ) : t -> t -> bool

     

       50
       50
       +
       

     

       51
       51
       +
       (** {2 Formatting} *)

     

       52
       52
       +
       

     

       53
       53
       +
       val to_string : t -> string

     

       54
       54
       +
       val format : ?precision:int -> t -> string

     

       55
       55
       +
       val to_dms_string : t -> string

     

       56
       56
       +
       (** Format as "DDD°MM'SS.S\"E" *)

+15

stack/lonlat/lib/lonlat.ml

···

       1
       1
       +
       (** LonLat - Geographic coordinate manipulation library for OCaml *)

     

       2
       2
       +
       

     

       3
       3
       +
       module Error = Error

     

       4
       4
       +
       module Lat = Lat

     

       5
       5
       +
       module Lon = Lon

     

       6
       6
       +
       module Coord = Coord

     

       7
       7
       +
       module Bbox = Bbox

     

       8
       8
       +
       module Vector = Vector

     

       9
       9
       +
       module Parser = Parser

     

       10
       10
       +
       

     

       11
       11
       +
       (** Constants *)

     

       12
       12
       +
       let earth_radius_km = 6371.0

     

       13
       13
       +
       let wgs84_a = 6378.137

     

       14
       14
       +
       let wgs84_b = 6356.752

     

       15
       15
       +
       let wgs84_f = 1.0 /. 298.257223563

+23

stack/lonlat/lib/lonlat.mli

···

       1
       1
       +
       (** LonLat - Geographic coordinate manipulation library for OCaml *)

     

       2
       2
       +
       

     

       3
       3
       +
       module Error = Error

     

       4
       4
       +
       module Lat = Lat

     

       5
       5
       +
       module Lon = Lon

     

       6
       6
       +
       module Coord = Coord

     

       7
       7
       +
       module Bbox = Bbox

     

       8
       8
       +
       module Vector = Vector

     

       9
       9
       +
       module Parser = Parser

     

       10
       10
       +
       

     

       11
       11
       +
       (** {1 Constants} *)

     

       12
       12
       +
       

     

       13
       13
       +
       val earth_radius_km : float

     

       14
       14
       +
       (** Mean Earth radius: 6371.0 km *)

     

       15
       15
       +
       

     

       16
       16
       +
       val wgs84_a : float

     

       17
       17
       +
       (** WGS84 semi-major axis: 6378.137 km *)

     

       18
       18
       +
       

     

       19
       19
       +
       val wgs84_b : float

     

       20
       20
       +
       (** WGS84 semi-minor axis: 6356.752 km *)

     

       21
       21
       +
       

     

       22
       22
       +
       val wgs84_f : float

     

       23
       23
       +
       (** WGS84 flattening: 1/298.257223563 *)

+134

stack/lonlat/lib/parser.ml

···

       1
       1
       +
       (** Coordinate parsing module *)

     

       2
       2
       +
       

     

       3
       3
       +
       let parse_dms s =

     

       4
       4
       +
         let s = String.trim s in

     

       5
       5
       +
       

     

       6
       6
       +
         (* Check for hemisphere markers *)

     

       7
       7
       +
         let has_s = String.contains s 'S' || String.contains s 's' in

     

       8
       8
       +
         let has_w = String.contains s 'W' || String.contains s 'w' in

     

       9
       9
       +
       

     

       10
       10
       +
         let sign =

     

       11
       11
       +
           if has_s || has_w then -1.0

     

       12
       12
       +
           else 1.0 in

     

       13
       13
       +
       

     

       14
       14
       +
         (* Remove hemisphere markers and clean up *)

     

       15
       15
       +
         let cleaned =

     

       16
       16
       +
           s

     

       17
       17
       +
           |> Str.global_replace (Str.regexp "[NnSsEeWw]") ""

     

       18
       18
       +
           |> String.trim in

     

       19
       19
       +
       

     

       20
       20
       +
         (* Try to parse different formats *)

     

       21
       21
       +
         try

     

       22
       22
       +
           (* First try as simple decimal degrees *)

     

       23
       23
       +
           if not (String.contains cleaned (Char.chr 176) ||  (* degree symbol *)

     

       24
       24
       +
                   String.contains cleaned '\'' ||

     

       25
       25
       +
                   String.contains cleaned '"' ||

     

       26
       26
       +
                   String.contains cleaned ' ' ||

     

       27
       27
       +
                   String.contains cleaned ':') then

     

       28
       28
       +
             sign *. float_of_string cleaned

     

       29
       29
       +
           else

     

       30
       30
       +
             (* Parse DMS format *)

     

       31
       31
       +
             let deg_char = Char.chr 176 in  (* degree symbol *)

     

       32
       32
       +
             let parts =

     

       33
       33
       +
               if String.contains cleaned deg_char then

     

       34
       34
       +
                 (* Format: DD degree MM min SS sec or DD degree MM.MMM min *)

     

       35
       35
       +
                 let deg_parts = String.split_on_char deg_char cleaned in

     

       36
       36
       +
                 match deg_parts with

     

       37
       37
       +
                 | [deg_str; rest] ->

     

       38
       38
       +
                     let deg = float_of_string (String.trim deg_str) in

     

       39
       39
       +
                     if String.contains rest '\'' then

     

       40
       40
       +
                       let min_parts = String.split_on_char '\'' rest in

     

       41
       41
       +
                       match min_parts with

     

       42
       42
       +
                       | [min_str; sec_str] ->

     

       43
       43
       +
                           let min = float_of_string (String.trim min_str) in

     

       44
       44
       +
                           let sec_str = String.trim sec_str in

     

       45
       45
       +
                           let sec_str =

     

       46
       46
       +
                             if String.contains sec_str '"' then

     

       47
       47
       +
                               String.sub sec_str 0 (String.index sec_str '"')

     

       48
       48
       +
                             else sec_str in

     

       49
       49
       +
                           let sec =

     

       50
       50
       +
                             if sec_str = "" then 0.0

     

       51
       51
       +
                             else float_of_string sec_str in

     

       52
       52
       +
                           [deg; min; sec]

     

       53
       53
       +
                       | [min_str] ->

     

       54
       54
       +
                           let min = float_of_string (String.trim min_str) in

     

       55
       55
       +
                           [deg; min; 0.0]

     

       56
       56
       +
                       | _ -> raise (Error.Coordinate_error

     

       57
       57
       +
                           (Error.Parse_error "Invalid DMS format"))

     

       58
       58
       +
                     else

     

       59
       59
       +
                       [deg; 0.0; 0.0]

     

       60
       60
       +
                 | _ -> raise (Error.Coordinate_error

     

       61
       61
       +
                     (Error.Parse_error "Invalid degree format"))

     

       62
       62
       +
               else if String.contains cleaned ':' then

     

       63
       63
       +
                 (* Format: DD:MM:SS *)

     

       64
       64
       +
                 cleaned

     

       65
       65
       +
                 |> String.split_on_char ':'

     

       66
       66
       +
                 |> List.map String.trim

     

       67
       67
       +
                 |> List.map float_of_string

     

       68
       68
       +
               else if String.contains cleaned ' ' then

     

       69
       69
       +
                 (* Format: DD MM SS *)

     

       70
       70
       +
                 cleaned

     

       71
       71
       +
                 |> String.split_on_char ' '

     

       72
       72
       +
                 |> List.filter (fun s -> String.length s > 0)

     

       73
       73
       +
                 |> List.map String.trim

     

       74
       74
       +
                 |> List.map float_of_string

     

       75
       75
       +
               else

     

       76
       76
       +
                 raise (Error.Coordinate_error

     

       77
       77
       +
                   (Error.Parse_error "Unknown coordinate format"))

     

       78
       78
       +
             in

     

       79
       79
       +
       

     

       80
       80
       +
             (* Convert to decimal degrees *)

     

       81
       81
       +
             match parts with

     

       82
       82
       +
             | [deg] -> sign *. deg

     

       83
       83
       +
             | [deg; min] -> sign *. (abs_float deg +. min /. 60.0)

     

       84
       84
       +
             | [deg; min; sec] ->

     

       85
       85
       +
                 sign *. (abs_float deg +. min /. 60.0 +. sec /. 3600.0)

     

       86
       86
       +
             | _ -> raise (Error.Coordinate_error

     

       87
       87
       +
                 (Error.Parse_error "Invalid number of components"))

     

       88
       88
       +
         with

     

       89
       89
       +
         | Failure _ -> raise (Error.Coordinate_error

     

       90
       90
       +
             (Error.Parse_error ("Failed to parse coordinate: " ^ s)))

     

       91
       91
       +
       

     

       92
       92
       +
       let parse_lat s =

     

       93
       93
       +
         let value = parse_dms s in

     

       94
       94
       +
         Lat.create value

     

       95
       95
       +
       

     

       96
       96
       +
       let parse_lon s =

     

       97
       97
       +
         let value = parse_dms s in

     

       98
       98
       +
         Lon.create value

     

       99
       99
       +
       

     

       100
       100
       +
       let parse_coord s =

     

       101
       101
       +
         let s = String.trim s in

     

       102
       102
       +
       

     

       103
       103
       +
         (* Try comma-separated first *)

     

       104
       104
       +
         let parts =

     

       105
       105
       +
           if String.contains s ',' then

     

       106
       106
       +
             String.split_on_char ',' s

     

       107
       107
       +
             |> List.map String.trim

     

       108
       108
       +
           else

     

       109
       109
       +
             (* Try space-separated *)

     

       110
       110
       +
             s

     

       111
       111
       +
             |> String.split_on_char ' '

     

       112
       112
       +
             |> List.filter (fun s -> String.length s > 0)

     

       113
       113
       +
             |> List.map String.trim

     

       114
       114
       +
         in

     

       115
       115
       +
       

     

       116
       116
       +
         match parts with

     

       117
       117
       +
         | [lat_str; lon_str] ->

     

       118
       118
       +
             let lat = parse_dms lat_str in

     

       119
       119
       +
             let lon = parse_dms lon_str in

     

       120
       120
       +
             Coord.create ~lat ~lon

     

       121
       121
       +
         | _ -> raise (Error.Coordinate_error

     

       122
       122
       +
             (Error.Parse_error "Expected lat,lon or lat lon format"))

     

       123
       123
       +
       

     

       124
       124
       +
       let try_parse_lat s =

     

       125
       125
       +
         try Some (parse_lat s)

     

       126
       126
       +
         with Error.Coordinate_error _ -> None

     

       127
       127
       +
       

     

       128
       128
       +
       let try_parse_lon s =

     

       129
       129
       +
         try Some (parse_lon s)

     

       130
       130
       +
         with Error.Coordinate_error _ -> None

     

       131
       131
       +
       

     

       132
       132
       +
       let try_parse_coord s =

     

       133
       133
       +
         try Some (parse_coord s)

     

       134
       134
       +
         with Error.Coordinate_error _ -> None

+23

stack/lonlat/lib/parser.mli

···

       1
       1
       +
       (** Coordinate parsing module *)

     

       2
       2
       +
       

     

       3
       3
       +
       val parse_lat : string -> Lat.t

     

       4
       4
       +
       (** Parse decimal degrees or DMS formats *)

     

       5
       5
       +
       

     

       6
       6
       +
       val parse_lon : string -> Lon.t

     

       7
       7
       +
       

     

       8
       8
       +
       val parse_coord : string -> Coord.t

     

       9
       9
       +
       (** Parse "lat,lon" or "lat lon" formats *)

     

       10
       10
       +
       

     

       11
       11
       +
       val parse_dms : string -> float

     

       12
       12
       +
       (** Parse various DMS formats:

     

       13
       13
       +
           - "40.7128"

     

       14
       14
       +
           - "40°42'46\""

     

       15
       15
       +
           - "40°42.767'"

     

       16
       16
       +
           - "40 42 46"

     

       17
       17
       +
           - "40:42:46"

     

       18
       18
       +
           - "40.7128N" or "N40.7128" *)

     

       19
       19
       +
       

     

       20
       20
       +
       (** Try parsing - returns None on failure instead of raising *)

     

       21
       21
       +
       val try_parse_lat : string -> Lat.t option

     

       22
       22
       +
       val try_parse_lon : string -> Lon.t option

     

       23
       23
       +
       val try_parse_coord : string -> Coord.t option

+71

stack/lonlat/lib/vector.ml

···

       1
       1
       +
       (** Geographic vector module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t = {

     

       4
       4
       +
         dx : float;  (* km *)

     

       5
       5
       +
         dy : float;  (* km *)

     

       6
       6
       +
       }

     

       7
       7
       +
       

     

       8
       8
       +
       let deg_to_rad = Float.pi /. 180.0

     

       9
       9
       +
       let rad_to_deg = 180.0 /. Float.pi

     

       10
       10
       +
       

     

       11
       11
       +
       let create_cartesian ~dx ~dy = { dx; dy }

     

       12
       12
       +
       

     

       13
       13
       +
       let create_polar ~heading ~distance =

     

       14
       14
       +
         (* Convert heading to mathematical angle (east = 0, counterclockwise) *)

     

       15
       15
       +
         let angle = (90.0 -. heading) *. deg_to_rad in

     

       16
       16
       +
         {

     

       17
       17
       +
           dx = distance *. cos angle;

     

       18
       18
       +
           dy = distance *. sin angle;

     

       19
       19
       +
         }

     

       20
       20
       +
       

     

       21
       21
       +
       let between p1 p2 =

     

       22
       22
       +
         (* Calculate vector between two coordinates *)

     

       23
       23
       +
         let dist = Coord.distance_haversine p1 p2 in

     

       24
       24
       +
         let bearing = Coord.bearing p1 p2 in

     

       25
       25
       +
         create_polar ~heading:bearing ~distance:dist

     

       26
       26
       +
       

     

       27
       27
       +
       let dx t = t.dx

     

       28
       28
       +
       let dy t = t.dy

     

       29
       29
       +
       

     

       30
       30
       +
       let heading t =

     

       31
       31
       +
         (* Convert from mathematical angle to geographic heading *)

     

       32
       32
       +
         let angle = atan2 t.dy t.dx *. rad_to_deg in

     

       33
       33
       +
         mod_float (90.0 -. angle +. 360.0) 360.0

     

       34
       34
       +
       

     

       35
       35
       +
       let distance t =

     

       36
       36
       +
         sqrt (t.dx ** 2.0 +. t.dy ** 2.0)

     

       37
       37
       +
       

     

       38
       38
       +
       let add v1 v2 = {

     

       39
       39
       +
         dx = v1.dx +. v2.dx;

     

       40
       40
       +
         dy = v1.dy +. v2.dy;

     

       41
       41
       +
       }

     

       42
       42
       +
       

     

       43
       43
       +
       let sub v1 v2 = {

     

       44
       44
       +
         dx = v1.dx -. v2.dx;

     

       45
       45
       +
         dy = v1.dy -. v2.dy;

     

       46
       46
       +
       }

     

       47
       47
       +
       

     

       48
       48
       +
       let scale v factor = {

     

       49
       49
       +
         dx = v.dx *. factor;

     

       50
       50
       +
         dy = v.dy *. factor;

     

       51
       51
       +
       }

     

       52
       52
       +
       

     

       53
       53
       +
       let neg v = {

     

       54
       54
       +
         dx = -.v.dx;

     

       55
       55
       +
         dy = -.v.dy;

     

       56
       56
       +
       }

     

       57
       57
       +
       

     

       58
       58
       +
       let rotate v degrees =

     

       59
       59
       +
         let angle = degrees *. deg_to_rad in

     

       60
       60
       +
         let cos_a = cos angle in

     

       61
       61
       +
         let sin_a = sin angle in

     

       62
       62
       +
         {

     

       63
       63
       +
           dx = v.dx *. cos_a -. v.dy *. sin_a;

     

       64
       64
       +
           dy = v.dx *. sin_a +. v.dy *. cos_a;

     

       65
       65
       +
         }

     

       66
       66
       +
       

     

       67
       67
       +
       let apply coord v =

     

       68
       68
       +
         Coord.offset coord ~heading:(heading v) ~distance:(distance v)

     

       69
       69
       +
       

     

       70
       70
       +
       let to_string t =

     

       71
       71
       +
         Printf.sprintf "Vector(heading=%.2f°, distance=%.2fkm)" (heading t) (distance t)

+39

stack/lonlat/lib/vector.mli

···

       1
       1
       +
       (** Geographic vector module *)

     

       2
       2
       +
       

     

       3
       3
       +
       type t

     

       4
       4
       +
       

     

       5
       5
       +
       (** {2 Construction} *)

     

       6
       6
       +
       

     

       7
       7
       +
       val create_cartesian : dx:float -> dy:float -> t

     

       8
       8
       +
       (** dx, dy in kilometers *)

     

       9
       9
       +
       

     

       10
       10
       +
       val create_polar : heading:float -> distance:float -> t

     

       11
       11
       +
       (** heading in degrees, distance in km *)

     

       12
       12
       +
       

     

       13
       13
       +
       val between : Coord.t -> Coord.t -> t

     

       14
       14
       +
       (** Vector from first to second point *)

     

       15
       15
       +
       

     

       16
       16
       +
       (** {2 Accessors} *)

     

       17
       17
       +
       

     

       18
       18
       +
       val dx : t -> float

     

       19
       19
       +
       val dy : t -> float

     

       20
       20
       +
       val heading : t -> float

     

       21
       21
       +
       val distance : t -> float

     

       22
       22
       +
       

     

       23
       23
       +
       (** {2 Operations} *)

     

       24
       24
       +
       

     

       25
       25
       +
       val add : t -> t -> t

     

       26
       26
       +
       val sub : t -> t -> t

     

       27
       27
       +
       val scale : t -> float -> t

     

       28
       28
       +
       val neg : t -> t

     

       29
       29
       +
       val rotate : t -> float -> t

     

       30
       30
       +
       (** Rotate by degrees *)

     

       31
       31
       +
       

     

       32
       32
       +
       (** {2 Application} *)

     

       33
       33
       +
       

     

       34
       34
       +
       val apply : Coord.t -> t -> Coord.t

     

       35
       35
       +
       (** Apply vector to coordinate *)

     

       36
       36
       +
       

     

       37
       37
       +
       (** {2 Formatting} *)

     

       38
       38
       +
       

     

       39
       39
       +
       val to_string : t -> string

+27

stack/lonlat/lonlat.opam

···

       1
       1
       +
       # This file is generated by dune, edit dune-project instead

     

       2
       2
       +
       opam-version: "2.0"

     

       3
       3
       +
       synopsis: "Geographic coordinate manipulation library for OCaml"

     

       4
       4
       +
       description:

     

       5
       5
       +
         "LonLat provides types and operations for working with latitude/longitude coordinates, including distance calculations, bearings, bounding boxes, and coordinate parsing."

     

       6
       6
       +
       maintainer: ["maintainer@example.com"]

     

       7
       7
       +
       authors: ["LonLat Authors"]

     

       8
       8
       +
       license: "MIT"

     

       9
       9
       +
       depends: [

     

       10
       10
       +
         "ocaml"

     

       11
       11
       +
         "dune" {>= "3.0"}

     

       12
       12
       +
         "odoc" {with-doc}

     

       13
       13
       +
       ]

     

       14
       14
       +
       build: [

     

       15
       15
       +
         ["dune" "subst"] {dev}

     

       16
       16
       +
         [

     

       17
       17
       +
           "dune"

     

       18
       18
       +
           "build"

     

       19
       19
       +
           "-p"

     

       20
       20
       +
           name

     

       21
       21
       +
           "-j"

     

       22
       22
       +
           jobs

     

       23
       23
       +
           "@install"

     

       24
       24
       +
           "@runtest" {with-test}

     

       25
       25
       +
           "@doc" {with-doc}

     

       26
       26
       +
         ]

     

       27
       27
       +
       ]