Merge pull request #5656 from tchak/fix-geo-areas

Fix geo areas computations

app/models/geo_area.rb

@@ -78,14 +78,14 @@ class GeoArea < ApplicationRecord
   end
 
   def area
-    if polygon? && RGeo::Geos.supported?
-      rgeo_geometry&.area&.round(1)
+    if polygon?
+      GeojsonService.area(geometry.deep_symbolize_keys).round(1)
     end
   end
 
   def length
-    if line? && RGeo::Geos.supported?
-      rgeo_geometry.length.round(1)
+    if line?
+      GeojsonService.length(geometry.deep_symbolize_keys).round(1)
     end
   end
 

app/services/geojson_service.rb

@@ -14,4 +14,117 @@ class GeojsonService
 
     polygon.to_json
   end
+
+  # The following code is ported from turfjs
+  # https://github.com/Turfjs/turf/blob/master/packages/turf-area/index.ts
+
+  EQUATORIAL_RADIUS = 6378137
+
+  def self.area(geojson)
+    calculate_area(geojson)
+  end
+
+  def self.length(geojson)
+    segment_reduce(geojson, 0) do |previous_value, segment|
+      coordinates = segment[:geometry][:coordinates]
+      previous_value + distance(coordinates[0], coordinates[1])
+    end
+  end
+
+  def self.distance(from, to)
+    coordinates1 = from
+    coordinates2 = to
+    d_lat = degrees_to_radians(coordinates2[1] - coordinates1[1])
+    d_lon = degrees_to_radians(coordinates2[0] - coordinates1[0])
+    lat1 = degrees_to_radians(coordinates1[1])
+    lat2 = degrees_to_radians(coordinates2[1])
+
+    a = (Math.sin(d_lat / 2)**2) + (Math.sin(d_lon / 2)**2) * Math.cos(lat1) * Math.cos(lat2)
+
+    radians = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
+    radians * EQUATORIAL_RADIUS
+  end
+
+  def self.calculate_area(geom)
+    total = 0
+    case geom[:type]
+    when 'Polygon'
+      polygon_area(geom[:coordinates]);
+    when 'MultiPolygon'
+      geom[:coordinates].each do |coordinates|
+        total += polygon_area(coordinates)
+      end
+      total
+    else
+      total
+    end
+  end
+
+  def self.polygon_area(coordinates)
+    total = 0
+    if coordinates.present?
+      coordinates = coordinates.dup
+      total += ring_area(coordinates.shift).abs
+      coordinates.each do |coordinates|
+        total -= ring_area(coordinates).abs
+      end
+    end
+    total
+  end
+
+  def self.ring_area(coordinates)
+    total = 0
+    coords_length = coordinates.size
+
+    if coords_length > 2
+      coords_length.times do |i|
+        if i == coords_length - 2 # i = N-2
+          lower_index = coords_length - 2
+          middle_index = coords_length - 1
+          upper_index = 0
+        elsif i == coords_length - 1 # i = N-1
+          lower_index = coords_length - 1
+          middle_index = 0
+          upper_index = 1
+        else # i = 0 to N-3
+          lower_index = i
+          middle_index = i + 1
+          upper_index = i + 2
+        end
+        p1 = coordinates[lower_index]
+        p2 = coordinates[middle_index]
+        p3 = coordinates[upper_index]
+        total += (rad(p3[0]) - rad(p1[0])) * Math.sin(rad(p2[1]))
+      end
+      total = total * EQUATORIAL_RADIUS * EQUATORIAL_RADIUS / 2
+    end
+
+    total
+  end
+
+  def self.segment_reduce(geojson, initial_value)
+    previous_value = initial_value
+    started = false
+    coordinates = geojson[:coordinates].dup
+    from = coordinates.shift
+    coordinates.each do |to|
+      current_segment = { type: 'Feature', geometry: { type: 'LineString', coordinates: [from, to] } }
+      from = to
+      if started == false && initial_value.blank?
+        previous_value = current_segment
+      else
+        previous_value = yield previous_value, current_segment
+      end
+      started = true
+    end
+    previous_value
+  end
+
+  def self.rad(num)
+    num * Math::PI / 180
+  end
+
+  def self.degrees_to_radians(degrees)
+    rad(degrees % 360)
+  end
 end

spec/models/dossier_spec.rb

@@ -1305,7 +1305,7 @@ describe Dossier do
               'type' => 'Polygon'
             },
             properties: {
-              area: 219.0,
+              area: 103.6,
               champ_id: champ_carte.stable_id,
               dossier_id: dossier.id,
               id: geo_area.id,

spec/models/geo_area_spec.rb

@@ -2,21 +2,19 @@ RSpec.describe GeoArea, type: :model do
   describe '#area' do
     let(:geo_area) { build(:geo_area, :polygon) }
 
-    it { expect(geo_area.area).to eq(219.0) }
+    it { expect(geo_area.area).to eq(103.6) }
   end
 
   describe '#area (hourglass polygon)' do
     let(:geo_area) { build(:geo_area, :hourglass_polygon) }
 
-    # This test fails in my local environement end the problem exists in production.
-    # Must be some mismatch between CI/production. I still want this fix in production.
-    it.pending { expect(geo_area.area).to be_nil }
+    it { expect(geo_area.area).to eq(32.4) }
   end
 
   describe '#length' do
     let(:geo_area) { build(:geo_area, :line_string) }
 
-    it { expect(geo_area.length).to eq(30.8) }
+    it { expect(geo_area.length).to eq(21.2) }
   end
 
   describe '#location' do