Merge pull request #10 from geocrystal/bearing

add a method to calculate the initial bearing between two points

spec/geo/bearing_spec.cr

@@ -0,0 +1,18 @@
+require "../spec_helper"
+require "../../src/geo/bearing.cr"
+
+describe Geo::Coord do
+  point_a = Geo::Coord.new(48.8566, 2.3522)   # Paris
+  point_b = Geo::Coord.new(40.7128, -74.0060) # New York
+
+  context "bearing" do
+    # Checked on https://www.movable-type.co.uk/scripts/latlong.html
+    context "calculates initial bearing" do
+      it { point_a.bearing(point_b).should eq(291.7938627483058) }
+    end
+
+    context "calculates final bearing" do
+      it { point_a.bearing(point_b, true).should eq(233.70448129781204) }
+    end
+  end
+end

src/geo/bearing.cr

@@ -0,0 +1,35 @@
+module Geo
+  struct Coord
+    # The formula used for calculating the initial bearing between two points on the Earth's surfaceis
+    # is derived from the broader concepts of great-circle distance and navigation on a spherical Earth model.
+    # https://en.wikipedia.org/wiki/Great-circle_distance
+    #
+    # https://github.com/Turfjs/turf/blob/master/packages/turf-bearing/index.ts
+    def bearing(to : Geo::Coord, final = false) : Float64
+      if final
+        # Calculate the bearing from the destination point back to the original point
+        reverse_bearing = calculate_bearing(to.lat, to.lng, lat, lng)
+
+        # Adjust by 180 degrees to get the final bearing in the correct direction
+        (reverse_bearing + 180) % 360
+      else
+        calculate_bearing(lat, lng, to.lat, to.lng)
+      end
+    end
+
+    private def calculate_bearing(lat1, lng1, lat2, lng2) : Float64
+      rad_lat1 = Geo::Utils.degrees_to_radians(lat1)
+      rad_lng1 = Geo::Utils.degrees_to_radians(lng1)
+      rad_lat2 = Geo::Utils.degrees_to_radians(lat2)
+      rad_lng2 = Geo::Utils.degrees_to_radians(lng2)
+
+      delta_lng = rad_lng2 - rad_lng1
+
+      a = Math.sin(delta_lng) * Math.cos(rad_lat2)
+      b = Math.cos(rad_lat1) * Math.sin(rad_lat2) -
+          Math.sin(rad_lat1) * Math.cos(rad_lat2) * Math.cos(delta_lng)
+
+      Geo::Utils.radians_to_degrees(Math.atan2(a, b)) % 360
+    end
+  end
+end

src/geo/polygon.cr

@@ -27,7 +27,7 @@ module Geo
 
     # Evaluate area of a polygon using shoelace formula
     def area
-      @coords.each_cons(2).sum { |p| p[0].shoelace(p[1]) }.abs.fdiv(2)
+      @coords.each_cons_pair.sum { |lat, lng| lat.shoelace(lng) }.abs.fdiv(2)
     end
 
     # Order coords in lexicographical order.
@@ -61,20 +61,18 @@ module Geo
     def contains?(coord : Geo::Coord) : Bool
       last_coord = @coords.last
       odd_nodes = false
-      x = coord.lng
-      y = coord.lat
-
-      @coords.each do |p|
-        yi = p.lat
-        xi = p.lng
-        yj = last_coord.lat
-        xj = last_coord.lng
+      y, x = coord.ll
+
+      @coords.each do |iter_coord|
+        yi, xi = iter_coord.ll
+        yj, xj = last_coord.ll
+
         if yi < y && yj >= y ||
            yj < y && yi >= y
           odd_nodes = !odd_nodes if xi + (y - yi) / (yj - yi) * (xj - xi) < x
         end
 
-        last_coord = p
+        last_coord = iter_coord
       end
 
       odd_nodes
@@ -92,6 +90,7 @@ module Geo
 
     def ==(other : Geo::Polygon) : Bool
       min_size = Math.min(size, other.size)
+
       0.upto(min_size - 1) do |i|
         return false unless self[i] == other[i]
       end

src/geo/utils.cr

@@ -21,5 +21,13 @@ module Geo
       return 0 if val == 0 # colinear
       val > 0 ? 1 : 2      # clockwise or counterclockwise
     end
+
+    def degrees_to_radians(degrees : Number) : Float64
+      degrees * Math::PI / 180.0
+    end
+
+    def radians_to_degrees(radians : Number) : Float64
+      radians * 180.0 / Math::PI
+    end
   end
 end