Significantly improved performance of clustersDbscan

packages/turf-clusters-dbscan/index.ts

@@ -1,8 +1,8 @@
 import { GeoJsonProperties, FeatureCollection, Point } from "geojson";
 import { clone } from "@turf/clone";
-import { distance } from "@turf/distance";
-import { degreesToRadians, lengthToDegrees, Units } from "@turf/helpers";
-import { rbush as RBush } from "./lib/rbush-export.js";
+import { Units } from "@turf/helpers";
+import KDBush from "kdbush";
+import * as geokdbush from "geokdbush";
 
 /**
  * Point classification within the cluster.
@@ -66,6 +66,9 @@ function clustersDbscan(
   } = {}
 ): FeatureCollection<Point, DbscanProps> {
   // Input validation being handled by Typescript
+  // TODO oops! No it isn't. Typescript doesn't do runtime checking. We should
+  // re-enable these checks, though will have to wait for a major version bump
+  // as more restrictive checks could break currently working code.
   // collectionOf(points, 'Point', 'points must consist of a FeatureCollection of only Points');
   // if (maxDistance === null || maxDistance === undefined) throw new Error('maxDistance is required');
   // if (!(Math.sign(maxDistance) > 0)) throw new Error('maxDistance is invalid');
@@ -77,11 +80,13 @@ function clustersDbscan(
   // Defaults
   const minPoints = options.minPoints || 3;
 
-  // Calculate the distance in degrees for region queries
-  const latDistanceInDegrees = lengthToDegrees(maxDistance, options.units);
-
   // Create a spatial index
-  var tree = new RBush(points.features.length);
+  const kdIndex = new KDBush(points.features.length);
+  // Index each point for spatial queries
+  for (const point of points.features) {
+    kdIndex.add(point.geometry.coordinates[0], point.geometry.coordinates[1]);
+  }
+  kdIndex.finish();
 
   // Keeps track of whether a point has been visited or not.
   var visited = points.features.map((_) => false);
@@ -95,54 +100,22 @@ function clustersDbscan(
   // Keeps track of the clusterId for each point
   var clusterIds: number[] = points.features.map((_) => -1);
 
-  // Index each point for spatial queries
-  tree.load(
-    points.features.map((point, index) => {
-      var [x, y] = point.geometry.coordinates;
-      return {
-        minX: x,
-        minY: y,
-        maxX: x,
-        maxY: y,
-        index: index,
-      } as IndexedPoint;
-    })
-  );
-
   // Function to find neighbors of a point within a given distance
   const regionQuery = (index: number): IndexedPoint[] => {
     const point = points.features[index];
     const [x, y] = point.geometry.coordinates;
 
-    const minY = Math.max(y - latDistanceInDegrees, -90.0);
-    const maxY = Math.min(y + latDistanceInDegrees, 90.0);
-
-    const lonDistanceInDegrees = (function () {
-      // Handle the case where the bounding box crosses the poles
-      if (minY < 0 && maxY > 0) {
-        return latDistanceInDegrees;
-      }
-      if (Math.abs(minY) < Math.abs(maxY)) {
-        return latDistanceInDegrees / Math.cos(degreesToRadians(maxY));
-      } else {
-        return latDistanceInDegrees / Math.cos(degreesToRadians(minY));
-      }
-    })();
-
-    const minX = Math.max(x - lonDistanceInDegrees, -360.0);
-    const maxX = Math.min(x + lonDistanceInDegrees, 360.0);
-
-    // Calculate the bounding box for the region query
-    const bbox = { minX, minY, maxX, maxY };
-    return (tree.search(bbox) as ReadonlyArray<IndexedPoint>).filter(
-      (neighbor) => {
-        const neighborIndex = neighbor.index;
-        const neighborPoint = points.features[neighborIndex];
-        const distanceInKm = distance(point, neighborPoint, {
-          units: "kilometers",
-        });
-        return distanceInKm <= maxDistance;
-      }
+    return (
+      geokdbush
+        // @ts-expect-error - until https://github.com/mourner/geokdbush/issues/20 is resolved
+        .around<number>(kdIndex, x, y, undefined, maxDistance)
+        .map((id) => ({
+          minX: points.features[id].geometry.coordinates[0],
+          minY: points.features[id].geometry.coordinates[1],
+          maxX: points.features[id].geometry.coordinates[0],
+          maxY: points.features[id].geometry.coordinates[1],
+          index: id,
+        }))
     );
   };
 

packages/turf-clusters-dbscan/package.json

@@ -63,7 +63,6 @@
     "@turf/centroid": "workspace:*",
     "@turf/clusters": "workspace:*",
     "@types/benchmark": "^2.1.5",
-    "@types/rbush": "^3.0.2",
     "@types/tape": "^5.8.1",
     "benchmark": "^2.1.4",
     "chromatism": "^3.0.0",
@@ -78,11 +77,12 @@
   },
   "dependencies": {
     "@turf/clone": "workspace:*",
-    "@turf/distance": "workspace:*",
     "@turf/helpers": "workspace:*",
     "@turf/meta": "workspace:*",
     "@types/geojson": "^7946.0.10",
-    "rbush": "^3.0.1",
+    "@types/geokdbush": "^1.1.5",
+    "geokdbush": "^2.0.1",
+    "kdbush": "^4.0.2",
     "tslib": "^2.8.1"
   }
 }

packages/turf-clusters-dbscan/test/in/fiji.geojson

@@ -3,41 +3,47 @@
   "features": [
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [179.439697265625, -16.55196172197251]
       }
     },
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [179.01123046874997, -16.97274101999901]
       }
     },
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [179.505615234375, -17.035777250427184]
       }
     },
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [180.75805664062497, -16.41500926733237]
       }
     },
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [181.1865234375, -16.615137799987075]
       }
     },
     {
       "type": "Feature",
+      "properties": {},
       "geometry": {
         "type": "Point",
         "coordinates": [181.03271484375, -16.277960306212513]