feat: Add geographic methods

This adds methods to compute latitude and longitude deltas from a
distance, which is typically useful to be able to query close
coordinates.
We also add a centroid computation in a sphere, that is uesful to
determine the center location of a set of geo points.

packages/cozy-client/src/models/geo.js

@@ -1,3 +1,6 @@
+const EARTH_RADIUS_M = 6378137
+const EARTH_CIRCUMFERENCE_M = 40075000
+
 /**
  * Compute the speed from distance and duration
  *
@@ -16,6 +19,15 @@ const degreesToRadians = degrees => {
   return (degrees * Math.PI) / 180
 }
 
+const radiansToDegrees = radians => {
+  return (radians * 180) / Math.PI
+}
+
+const roundToNDecimals = (value, n) => {
+  const multiplier = Math.pow(10, n)
+  return Math.round(value * multiplier) / multiplier
+}
+
 /**
  * Compute the distance between 2 geographic points, in meters.
  *
@@ -54,7 +66,91 @@ export const geodesicDistance = (point1, point2) => {
   const a = aLat + Math.cos(lat1) * Math.cos(lat2) * aLon
   const c = 2 * Math.asin(Math.sqrt(a))
 
-  // Radius of earth is 6371 km
-  const distance = 6371 * 1000 * c
-  return Math.round(distance * 100) / 100
+  const distance = EARTH_RADIUS_M * c
+  return roundToNDecimals(distance, 2)
+}
+
+/**
+ * Compute the geographical center of the given points
+ *
+ * This consists of finding the centroid of a set of points
+ * in a sphere.
+ * Note this assumes the Earth is a perfect sphere, which is not,
+ * but the approximation should be good enough.
+ *
+ * @param {Array<import("../types").Coordinates>} coordinates - The geo points
+ * @returns {import("../types").Coordinates} The center point
+ */
+export const computeSphericalCenter = coordinates => {
+  if (coordinates.length < 1) {
+    return null
+  }
+  if (coordinates.length === 1) {
+    return coordinates[0]
+  }
+  let totalX = 0
+  let totalY = 0
+  let totalZ = 0
+
+  for (const coord of coordinates) {
+    let lon = degreesToRadians(coord.lon)
+    let lat = degreesToRadians(coord.lat)
+
+    // Convert spherical coordinates to Cartesian coordinates
+    let x = Math.cos(lat) * Math.cos(lon)
+    let y = Math.cos(lat) * Math.sin(lon)
+    let z = Math.sin(lat)
+
+    totalX += x
+    totalY += y
+    totalZ += z
+  }
+
+  const avgX = totalX / coordinates.length
+  const avgY = totalY / coordinates.length
+  const avgZ = totalZ / coordinates.length
+
+  // Don't forget to convert Cartesian coordinates back to spherical
+  const centralLon = radiansToDegrees(Math.atan2(avgY, avgX))
+  const hyp = Math.sqrt(avgX * avgX + avgY * avgY)
+  const centralLat = radiansToDegrees(Math.atan2(avgZ, hyp))
+  return {
+    lat: roundToNDecimals(centralLat, 13),
+    lon: roundToNDecimals(centralLon, 13)
+  }
+}
+
+/**
+ * Compute the longitude delta from a distance, in meters.
+ *
+ * This requires the latitude: we want to compute the horizontal delta
+ * on the Earth surface. As it is a sphere (kind of), this delta won't be
+ * the same depending on whether it is on the equator (min variation)
+ * or on the poles (max variation), for instance.
+ *
+ * @param {number} latitude - The latitude
+ * @param {number} distance - The distance in meters
+ * @returns {number} the longitude delta degrees
+ */
+export const deltaLongitude = (latitude, distance) => {
+  const phi = degreesToRadians(latitude)
+  const deltaLambda = distance / (EARTH_RADIUS_M * Math.cos(phi))
+
+  return roundToNDecimals(radiansToDegrees(deltaLambda), 13)
+}
+
+/**
+ * Compute the latitude delta from a distance, in meters.
+ *
+ * The reasoning is rather simple: there are 360° of latitudes of same distance.
+ * Then, it consists of computing 1 degree distance, and divide the
+ * given distance by this value.
+ *
+ * @param {number} distance - The distance in meters
+ * @returns {number} The delta latitude degrees
+ */
+export const deltaLatitude = distance => {
+  const distOneLatDegree = EARTH_CIRCUMFERENCE_M / 360 // 111 319 meters per degree
+  const deltaLat = distance / distOneLatDegree
+  return roundToNDecimals(deltaLat, 13)
 }

packages/cozy-client/src/models/geo.spec.js

@@ -1,4 +1,10 @@
-const { geodesicDistance, computeSpeed } = require('./geo')
+const {
+  geodesicDistance,
+  computeSpeed,
+  computeSphericalCenter,
+  deltaLongitude,
+  deltaLatitude
+} = require('./geo')
 
 describe('geodesicDistance', () => {
   it('should return 0 for 2 points with same coordinates', () => {
@@ -15,7 +21,7 @@ describe('geodesicDistance', () => {
   it('should return the correct Paris-NY distance', () => {
     const paris = { lat: 48.8566, lon: 2.3522 }
     const NY = { lat: 40.7128, lon: -74.006 }
-    expect(geodesicDistance(paris, NY)).toEqual(5837240.9)
+    expect(geodesicDistance(paris, NY)).toEqual(5843779.97)
   })
 })
 
@@ -30,3 +36,57 @@ describe('speed', () => {
     expect(computeSpeed(10000, 1800)).toEqual(5.56)
   })
 })
+
+describe('computeSphericalCenter', () => {
+  it('should return null when coordinates are missing', () => {
+    expect(computeSphericalCenter([])).toEqual(null)
+  })
+  it('should return the correct coordinates when there are all equals', () => {
+    expect(computeSphericalCenter([{ lat: 10, lon: 10 }])).toEqual({
+      lat: 10,
+      lon: 10
+    })
+    expect(
+      computeSphericalCenter([
+        { lat: 0.0, lon: 0.0 },
+        { lat: 0.0, lon: 0.0 },
+        { lat: 0.0, lon: 0.0 }
+      ])
+    ).toEqual({ lat: 0, lon: 0 })
+  })
+  it('should return correct center for n points', () => {
+    expect(
+      computeSphericalCenter([{ lat: 10, lon: 10 }, { lat: 10, lon: 10 }])
+    ).toEqual({ lat: 10, lon: 10 })
+
+    expect(
+      computeSphericalCenter([{ lat: 48, lon: 2 }, { lat: 48, lon: 3 }])
+    ).toEqual({ lat: 48.0010848667078, lon: 2.5 })
+
+    expect(
+      computeSphericalCenter([
+        { lat: 10, lon: 10 },
+        { lat: 20, lon: 20 },
+        { lat: 30, lon: 30 }
+      ])
+    ).toEqual({ lat: 20.1866988557583, lon: 19.5721919393414 })
+  })
+})
+
+describe('delta longitude and latitude', () => {
+  it('null distance should give a null delta', () => {
+    expect(deltaLongitude(46.66, 0)).toEqual(0)
+    expect(deltaLatitude(0)).toEqual(0)
+  })
+
+  it('should give the correct longitude delta', () => {
+    expect(deltaLongitude(46.66, 200)).toEqual(0.002617749975)
+    expect(deltaLongitude(0, 200)).toEqual(0.0017966305682) // equator
+    expect(deltaLongitude(89.99, 200)).toEqual(10.2939349426746) // north pole
+    expect(deltaLongitude(-89.99, 200)).toEqual(10.2939349426746) // south pole
+  })
+  it('should give the correct latitude delta', () => {
+    expect(deltaLatitude(200)).toEqual(0.0017966313163)
+    expect(deltaLatitude(2000)).toEqual(0.0179663131628) // contrarily to the longitude, this is linear
+  })
+})

packages/cozy-client/types/models/geo.d.ts

@@ -1,2 +1,5 @@
 export function computeSpeed(distance: number, duration: number): number;
 export function geodesicDistance(point1: import("../types").Coordinates, point2: import("../types").Coordinates): number;
+export function computeSphericalCenter(coordinates: Array<import("../types").Coordinates>): import("../types").Coordinates;
+export function deltaLongitude(latitude: number, distance: number): number;
+export function deltaLatitude(distance: number): number;