demarches-normaliennes/app/services/geojson_service.rb

class GeojsonService
  def self.to_json_polygon_for_cadastre(coordinates)
    polygon = {
      geom: {
        type: "Feature",
        geometry: {
          type: "Polygon",
          coordinates: [
            coordinates
          ]
        }
      }
    }

    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.bbox(geojson)
    result = [-Float::INFINITY, -Float::INFINITY, Float::INFINITY, Float::INFINITY]

    self.coord_each(geojson) do |coord|
      if result[3] > coord[1]
        result[3] = coord[1]
      end
      if result[2] > coord[0]
        result[2] = coord[0]
      end
      if result[1] < coord[1]
        result[1] = coord[1]
      end
      if result[0] < coord[0]
        result[0] = coord[0]
      end
    end

    result
  end

  def self.coord_each(geojson)
    geometries = if geojson.fetch(:type) == "FeatureCollection"
      geojson.fetch(:features).map { _1.fetch(:geometry) }
    else
      [geojson.fetch(:geometry)]
    end.compact

    geometries.each do |geometry|
      geometries = if geometry.fetch(:type) == "GeometryCollection"
        geometry.fetch(:geometries)
      else
        [geometry]
      end.compact

      geometries.each do |geometry|
        case geometry.fetch(:type)
        when "Point"
          yield geometry.fetch(:coordinates).map(&:to_f)
        when "LineString", "MultiPoint"
          geometry.fetch(:coordinates).each { yield _1.map(&:to_f) }
        when "Polygon", "MultiLineString"
          geometry.fetch(:coordinates).each do |shapes|
            shapes.each { yield _1.map(&:to_f) }
          end
        when "MultiPolygon"
          geometry.fetch(:coordinates).each do |polygons|
            polygons.each do |shapes|
              shapes.each { yield _1.map(&:to_f) }
            end
          end
        when "GeometryCollection"
          geometry.fetch(:geometries).each do |geometry|
            coord_each(geometry) { yield _1 }
          end
        end
      end
    end
  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