diff --git a/app/models/geo_area.rb b/app/models/geo_area.rb index c2693e80a..19b16ce21 100644 --- a/app/models/geo_area.rb +++ b/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 diff --git a/app/services/geojson_service.rb b/app/services/geojson_service.rb index 0e4c04914..d2a78d884 100644 --- a/app/services/geojson_service.rb +++ b/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 diff --git a/spec/models/dossier_spec.rb b/spec/models/dossier_spec.rb index b3f8966a8..0ab322673 100644 --- a/spec/models/dossier_spec.rb +++ b/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, diff --git a/spec/models/geo_area_spec.rb b/spec/models/geo_area_spec.rb index 4a42e972b..377b2bf00 100644 --- a/spec/models/geo_area_spec.rb +++ b/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