diff --git a/spec/geo/bearing_spec.cr b/spec/geo/bearing_spec.cr new file mode 100644 index 0000000..c7615f4 --- /dev/null +++ b/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 diff --git a/src/geo/bearing.cr b/src/geo/bearing.cr new file mode 100644 index 0000000..3798e5c --- /dev/null +++ b/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 diff --git a/src/geo/polygon.cr b/src/geo/polygon.cr index c711ff4..39f15ec 100644 --- a/src/geo/polygon.cr +++ b/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 diff --git a/src/geo/utils.cr b/src/geo/utils.cr index 0f27a57..2d7b5be 100644 --- a/src/geo/utils.cr +++ b/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