demarches-normaliennes/app/services/geojson_service.rb

190 lines
4.9 KiB
Ruby

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