[turf-ellipse] fixes the wrong shape when drawing "big ellipses" near the poles

packages/turf-ellipse/README.md

@@ -17,6 +17,7 @@ Takes a [Point][1] and calculates the ellipse polygon given two semi-axes expres
     *   `options.pivot` **[Coord][2]** point around which any rotation will be performed (optional, default `center`)
     *   `options.steps` **[number][3]** number of steps (optional, default `64`)
     *   `options.units` **[string][5]** unit of measurement for axes (optional, default `'kilometers'`)
+    *   `options.accuracy` **[number][3]** level of precision used for the repartition of points along the curve (optional, default `3`)
     *   `options.properties` **[Object][4]** properties (optional, default `{}`)
 
 ### Examples

packages/turf-ellipse/index.ts

@@ -5,8 +5,10 @@ import {
   isNumber,
   Coord,
   Units,
+  point,
 } from "@turf/helpers";
-import { rhumbDestination } from "@turf/rhumb-destination";
+import { destination } from "@turf/destination";
+import { distance } from "@turf/distance";
 import { transformRotate } from "@turf/transform-rotate";
 import { getCoord } from "@turf/invariant";
 import { GeoJsonProperties, Feature, Polygon } from "geojson";
@@ -22,6 +24,7 @@ import { GeoJsonProperties, Feature, Polygon } from "geojson";
  * @param {Coord} [options.pivot=center] point around which any rotation will be performed
  * @param {number} [options.steps=64] number of steps
  * @param {string} [options.units='kilometers'] unit of measurement for axes
+ * @param {number} [options.accuracy=3] level of precision used for the repartition of points along the curve
  * @param {Object} [options.properties={}] properties
  * @returns {Feature<Polygon>} ellipse polygon
  * @example
@@ -43,6 +46,7 @@ function ellipse(
     angle?: number;
     pivot?: Coord;
     properties?: GeoJsonProperties;
+    accuracy?: number;
   }
 ): Feature<Polygon> {
   // Optional params
@@ -52,71 +56,102 @@ function ellipse(
   const angle = options.angle || 0;
   const pivot = options.pivot || center;
   const properties = options.properties || {};
-
+  const accuracy = options.accuracy === undefined ? 0 : options.accuracy;
   // validation
   if (!center) throw new Error("center is required");
   if (!xSemiAxis) throw new Error("xSemiAxis is required");
   if (!ySemiAxis) throw new Error("ySemiAxis is required");
   if (!isObject(options)) throw new Error("options must be an object");
   if (!isNumber(steps)) throw new Error("steps must be a number");
   if (!isNumber(angle)) throw new Error("angle must be a number");
+  if (!isNumber(accuracy)) throw new Error("accuracy must be a number");
+  if (accuracy !== 0 && accuracy < 1)
+    throw new Error(
+      "accuracy must either be equal to -1 or greater or equal to 1"
+    );
 
-  const centerCoords = getCoord(center);
-  if (units !== "degrees") {
-    const xDest = rhumbDestination(center, xSemiAxis, 90, { units });
-    const yDest = rhumbDestination(center, ySemiAxis, 0, { units });
-    xSemiAxis = getCoord(xDest)[0] - centerCoords[0];
-    ySemiAxis = getCoord(yDest)[1] - centerCoords[1];
-  }
+  const internalSteps = Math.floor(Math.pow(accuracy, 2) * steps);
+
+  const centerCoords = getCoord(
+    transformRotate(point(getCoord(center)), angle, { pivot })
+  );
 
   const coordinates: number[][] = [];
-  for (let i = 0; i < steps; i += 1) {
-    const stepAngle = (i * -360) / steps;
-    let x =
-      (xSemiAxis * ySemiAxis) /
-      Math.sqrt(
-        Math.pow(ySemiAxis, 2) +
-          Math.pow(xSemiAxis, 2) * Math.pow(getTanDeg(stepAngle), 2)
+
+  let r = Math.sqrt(
+    (Math.pow(xSemiAxis, 2) * Math.pow(ySemiAxis, 2)) /
+      (Math.pow(xSemiAxis * Math.cos(degreesToRadians(-angle)), 2) +
+        Math.pow(ySemiAxis * Math.sin(degreesToRadians(-angle)), 2))
+  );
+  let currentCoords = getCoord(
+    destination(centerCoords, r, 0, { units: units })
+  );
+  coordinates.push(currentCoords);
+
+  let currentAngle = 0;
+  let currentArcLength = 0;
+  let previousCoords = currentCoords;
+  const cumulatedArcLength = [0];
+  const cumulatedAngle = [0];
+  let circumference = 0;
+  if (accuracy != 0) {
+    for (let i = 1; i < internalSteps + 1; i += 1) {
+      previousCoords = currentCoords;
+      currentAngle = (360 * i) / internalSteps;
+      r = Math.sqrt(
+        (Math.pow(xSemiAxis, 2) * Math.pow(ySemiAxis, 2)) /
+          (Math.pow(
+            xSemiAxis * Math.cos(degreesToRadians(currentAngle - angle)),
+            2
+          ) +
+            Math.pow(
+              ySemiAxis * Math.sin(degreesToRadians(currentAngle - angle)),
+              2
+            ))
       );
-    let y =
-      (xSemiAxis * ySemiAxis) /
-      Math.sqrt(
-        Math.pow(xSemiAxis, 2) +
-          Math.pow(ySemiAxis, 2) / Math.pow(getTanDeg(stepAngle), 2)
+      currentCoords = getCoord(
+        destination(centerCoords, r, currentAngle, { units: units })
       );
-
-    if (stepAngle < -90 && stepAngle >= -270) x = -x;
-    if (stepAngle < -180 && stepAngle >= -360) y = -y;
-    if (units === "degrees") {
-      const angleRad = degreesToRadians(angle);
-      const newx = x * Math.cos(angleRad) + y * Math.sin(angleRad);
-      const newy = y * Math.cos(angleRad) - x * Math.sin(angleRad);
-      x = newx;
-      y = newy;
+      currentArcLength += distance(previousCoords, currentCoords);
+      cumulatedAngle.push(currentAngle);
+      cumulatedArcLength.push(currentArcLength);
     }
+    circumference = cumulatedArcLength[cumulatedArcLength.length - 1];
+  }
 
-    coordinates.push([x + centerCoords[0], y + centerCoords[1]]);
+  let j = 0;
+  for (let i = 1; i < steps; i += 1) {
+    let angleNewPoint = (360 * i) / steps;
+    if (accuracy != 0) {
+      const targetArcLength = (i * circumference) / steps;
+      while (cumulatedArcLength[j] < targetArcLength) {
+        j += 1;
+      }
+      const ratio =
+        (targetArcLength - cumulatedArcLength[j - 1]) /
+        (cumulatedArcLength[j] - cumulatedArcLength[j - 1]);
+      angleNewPoint =
+        cumulatedAngle[j - 1] +
+        ratio * (cumulatedAngle[j] - cumulatedAngle[j - 1]);
+    }
+    r = Math.sqrt(
+      (Math.pow(xSemiAxis, 2) * Math.pow(ySemiAxis, 2)) /
+        (Math.pow(
+          xSemiAxis * Math.cos(degreesToRadians(angleNewPoint - angle)),
+          2
+        ) +
+          Math.pow(
+            ySemiAxis * Math.sin(degreesToRadians(angleNewPoint - angle)),
+            2
+          ))
+    );
+    currentCoords = getCoord(
+      destination(centerCoords, r, angleNewPoint, { units: units })
+    );
+    coordinates.push(currentCoords);
   }
   coordinates.push(coordinates[0]);
-  if (units === "degrees") {
-    return polygon([coordinates], properties);
-  } else {
-    return transformRotate(polygon([coordinates], properties), angle, {
-      pivot,
-    });
-  }
-}
-
-/**
- * Get Tan Degrees
- *
- * @private
- * @param {number} deg Degrees
- * @returns {number} Tan Degrees
- */
-function getTanDeg(deg: number) {
-  const rad = (deg * Math.PI) / 180;
-  return Math.tan(rad);
+  return polygon([coordinates], properties);
 }
 
 export { ellipse };

packages/turf-ellipse/package.json

@@ -53,9 +53,10 @@
   },
   "devDependencies": {
     "@placemarkio/check-geojson": "^0.1.12",
+    "@turf/area": "workspace:^",
     "@turf/bbox-polygon": "workspace:^",
     "@turf/circle": "workspace:^",
-    "@turf/destination": "workspace:^",
+    "@turf/intersect": "workspace:^",
     "@turf/truncate": "workspace:^",
     "@types/benchmark": "^2.1.5",
     "@types/tape": "^4.2.32",
@@ -70,9 +71,10 @@
     "write-json-file": "^5.0.0"
   },
   "dependencies": {
+    "@turf/destination": "workspace:^",
+    "@turf/distance": "workspace:^",
     "@turf/helpers": "workspace:^",
     "@turf/invariant": "workspace:^",
-    "@turf/rhumb-destination": "workspace:^",
     "@turf/transform-rotate": "workspace:^",
     "@types/geojson": "^7946.0.10",
     "tslib": "^2.6.2"