WIP geodetic distance search

nitrite-spatial/pom.xml

@@ -42,6 +42,11 @@
             <groupId>org.locationtech.jts</groupId>
             <artifactId>jts-core</artifactId>
         </dependency>
+        <dependency>
+            <groupId>net.sf.geographiclib</groupId>
+            <artifactId>GeographicLib-Java</artifactId>
+            <version>2.0</version>
+        </dependency>
         <dependency>
             <groupId>org.projectlombok</groupId>
             <artifactId>lombok</artifactId>

nitrite-spatial/src/main/java/org/dizitart/no2/spatial/NearFilter.java

@@ -19,7 +19,7 @@
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.Point;
-import org.locationtech.jts.util.GeometricShapeFactory;
+import org.locationtech.jts.util.GeometricShapeFactoryExt;
 
 /**
  * @since 4.0
@@ -35,7 +35,7 @@ class NearFilter extends WithinFilter {
     }
 
     private static Geometry createCircle(Coordinate center, double radius) {
-        GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
+        GeometricShapeFactoryExt shapeFactory = new GeometricShapeFactoryExt();
         shapeFactory.setNumPoints(64);
         shapeFactory.setCentre(center);
         shapeFactory.setSize(radius * 2);

nitrite-spatial/src/main/java/org/dizitart/no2/spatial/SpatialIndex.java

@@ -138,6 +138,19 @@ public LinkedHashSet<NitriteId> findNitriteIds(FindPlan findPlan) {
         Geometry geometry = spatialFilter.getValue();
         BoundingBox boundingBox = fromGeometry(geometry);
 
+        // NOTE: It looks like we are painted into something of a corner here! The index only knows the bounding boxes,
+        //    because that's how an R-Tree is meant to work. And that's fine until we have to deal with the reality
+        //    of the fact that these are actually arbitrary geometries and the users are crafting queries that are
+        //    meant to test whether they are within each other or not.
+        //
+        //    The query *cannot* be satisfied by simply checking bounding boxes and calling it a day.
+        //
+        //    The closest thing I see so far to a solution is for us to *split the filter into two steps* and either:
+        //    (a) plumb the Map<NitriteId,Document> from ReadOperations down to here, so that we can look at the actual
+        //    geometry values; or (b) treat this as an initial "fast filtering" step and then plan to have something
+        //    in ReadOperations (or even further up the call stack) take the next step and perform the more
+        //    computationally intensive geometry check.
+
         if (filter instanceof WithinFilter) {
             keys = indexMap.findContainedKeys(boundingBox);
         } else if (filter instanceof IntersectsFilter) {

nitrite-spatial/src/main/java/org/locationtech/jts/util/GeometricShapeFactoryExt.java

@@ -0,0 +1,45 @@
+package org.locationtech.jts.util;
+
+import net.sf.geographiclib.Geodesic;
+import net.sf.geographiclib.GeodesicData;
+import net.sf.geographiclib.GeodesicMask;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.LinearRing;
+import org.locationtech.jts.geom.Polygon;
+
+/** Extends the JTS GeometricShapeFactory to add geodetic "small circle" geometry. */
+public class GeometricShapeFactoryExt extends GeometricShapeFactory {
+
+    /**
+     * Bitmask specifying which results should be returned from the "direct" calculation.
+     * We don't need to know the azimuth at s2, so we can save a couple of CPU cycles by not asking for it!
+     * <p/>
+     * See javadoc at {@link Geodesic#Direct(double, double, double, double, int)} for more details.
+     */
+    public static final int OUTMASK = GeodesicMask.STANDARD ^ GeodesicMask.AZIMUTH;
+
+    @Override
+    public Polygon createCircle()
+    {
+        Envelope env = dim.getEnvelope();
+        double radiusInMeters = env.getWidth() / 2.0;
+        double ctrX = dim.getCentre().x;
+        double ctrY = dim.getCentre().y;
+
+        Coordinate[] pts = new Coordinate[nPts + 1];
+        int i;
+        for (i = 0; i < nPts; i++) {
+            // because this is the "azimuth" value, it starts at "geodetic north" and proceeds clockwise
+            double azimuthInDegrees = i * (360.0 / nPts);
+            GeodesicData directResult = Geodesic.WGS84.Direct(ctrX, ctrY, azimuthInDegrees, radiusInMeters, OUTMASK);
+            double lat = directResult.lat2;
+            double lon = directResult.lon2;
+            pts[i] = coord(lat, lon);
+        }
+        pts[i] = new Coordinate(pts[0]);
+
+        LinearRing ring = geomFact.createLinearRing(pts);
+        return geomFact.createPolygon(ring);
+    }
+}

nitrite-spatial/src/test/java/org/dizitart/no2/spatial/GeoSpatialFindNearTest.java

@@ -0,0 +1,180 @@
+/*
+ * Copyright (c) 2017-2021 Nitrite author or authors.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *       http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+package org.dizitart.no2.spatial;
+
+import net.sf.geographiclib.GeodesicData;
+import org.dizitart.no2.collection.DocumentCursor;
+import org.dizitart.no2.collection.FindOptions;
+import org.dizitart.no2.index.IndexOptions;
+import org.dizitart.no2.repository.Cursor;
+import org.junit.Test;
+import org.locationtech.jts.geom.Coordinate;
+
+import java.util.Random;
+import java.util.stream.Collectors;
+
+import static java.util.Collections.emptySet;
+import static net.sf.geographiclib.Geodesic.WGS84;
+import static org.dizitart.no2.collection.Document.createDocument;
+import static org.dizitart.no2.collection.FindOptions.orderBy;
+import static org.dizitart.no2.common.SortOrder.Ascending;
+import static org.dizitart.no2.common.util.Iterables.setOf;
+import static org.dizitart.no2.spatial.GeoSpatialTestBase.TestLocations.*;
+import static org.dizitart.no2.spatial.SpatialFluentFilter.where;
+import static org.dizitart.no2.spatial.SpatialIndexer.SPATIAL_INDEX;
+import static org.junit.Assert.assertEquals;
+
+/**
+ * Test cases to check whether the Spatial distance search is properly converting meters to/from degrees and
+ * accounting for the curvature of the Earth.
+ * <p>
+ * From the <a href="https://nitrite.dizitart.com/java-sdk/filter/index.html#near-filter">"near filter"
+ * documentation</a>, it is clear that the intended use is for the arguments to represent (1) a geo-coordinate
+ * as a lat/long pair; and (2) a distance in meters.
+ * <p>
+ * This appears to disagree with the example provided in <pre>nitrite-spatial/README.md</pre>, which makes no
+ * mention of distance units uses a coordinate value outside the space of possible lat/long coordinates.
+ **/
+public class GeoSpatialFindNearTest extends GeoSpatialTestBase {
+
+    /** 2.5 kilometers, in meters */
+    private static final double DIST_2500_M = 2_500.0d;
+    /** 20 millimeters, in meters */
+    private static final double DIST_20_MM = (20.0d / 1000.0d);
+
+    /** Just sorting by "id" field so that JUnit's "Expected" vs. "Actual" output is easier to visually inspect. */
+    public static final FindOptions ORDER_BY_ID = orderBy("id", Ascending);
+    public static final Random RANDOM = new Random();
+
+    @Test
+    public void testNearFilter_ObjectRepository_FindsCorrectSubset() {
+        repository.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "geometry");
+        Cursor<SpatialData> cursor = repository.find(
+            where("geometry").near(centerPt.getCoordinate(), DIST_2500_M),
+            ORDER_BY_ID);
+
+        assertEquals("Near filter should only find two locations within 2.5km",
+            setOf(obj2_950m_ESE, obj3_930m_W),
+            cursor.toSet());
+    }
+
+    @Test
+    public void testNearFilter_NitriteCollection_FindsCorrectSubset() {
+        collection.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "location");
+
+        DocumentCursor cursor = collection.find(
+            where("location").near(centerPt.getCoordinate(), DIST_2500_M),
+            ORDER_BY_ID);
+
+        assertEquals("Near filter should only find two locations within 2.5km",
+            setOf(doc2_950m_ESE, doc3_930m_W),
+            cursor.toList().stream().map(this::trimMeta).collect(Collectors.toSet()));
+    }
+
+    @Test
+    public void testNearFilter_ObjectRepository_SimpleEquatorDistances() {
+
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(0.0 1.00)"))); // 111km from (0,0)
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(0.0 0.01)"))); // 1.11km from (0,0)
+
+        repository.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "geometry");
+
+        Cursor<SpatialData> cursor = repository.find(where("geometry").near(new Coordinate(0.0, 0.0), 1120.0));
+
+        assertEquals("Near filter should find 1 location within 1.12km", 1, cursor.toList().size());
+    }
+
+    @Test
+    public void testNearFilter_NitriteCollection_SimpleEquatorDistances() {
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(0.0 1.00)")));
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(0.0 0.01)")));
+        collection.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "location");
+
+        DocumentCursor cursor = collection.find(where("location").near(new Coordinate(0.0, 0.0), 1120.0));
+
+        assertEquals("Near filter should find 1 location within 1.12km", 1, cursor.toList().size());
+    }
+
+    @Test
+    public void testNearFilter_ObjectRepository_Simple60NDistances() {
+
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(60.0 1.00)"))); // ~56km from (60,0)
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(60.0 2.00)"))); // ~111km from (60,0)
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(61.0 0.00)"))); // ~111km from (60,0)
+        repository.insert(new SpatialData(randLong(), readPoint("POINT(62.0 0.00)"))); // ~222km from (60,0)
+
+        repository.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "geometry");
+
+        Cursor<SpatialData> cursor = repository.find(
+            where("geometry").near(new Coordinate(60.0, 0.0), 112_000.0));
+
+        assertEquals("Near filter should find 3 locations within 1.12km", 3, cursor.toList().size());
+    }
+
+    @Test
+    public void testNearFilter_NitriteCollection_Simple60NDistances() {
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(60.0 1.00)")));  // ~56km from (60,0)
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(60.0 2.00)")));  // ~111km from (60,0)
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(61.0 0.00)")));  // ~111km from (60,0)
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(62.0 0.00)")));  // ~222km from (60,0)
+
+        // This point is outside the circle but inside the bounding box
+        collection.insert(createDocument("key", randLong()).put("location", readPoint("POINT(60.75 1.75)")));  // ~127km from (60,0)
+
+        collection.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "location");
+
+        DocumentCursor cursor = collection.find(where("location").near(new Coordinate(60.0, 0.0), 112_000.0));
+
+        assertEquals("Near filter should find 3 locations within 1.12km", 3, cursor.toList().size());
+    }
+
+    private static long randLong() {
+        return RANDOM.nextLong();
+    }
+
+    @Test
+    public void testNearFilter_ObjectRepository_TinyDistance_ReturnsNoResults() {
+        Coordinate coordinate = centerPt.getCoordinate();
+
+        Cursor<SpatialData> cursor = repository.find(where("geometry").near(coordinate, DIST_20_MM));
+        assertEquals("Near filter should return no results within 20mm distance", emptySet(), cursor.toSet());
+    }
+
+    @Test
+    public void testNearFilter_NitriteCollection_TinyDistance_ReturnsNoResults() {
+        Coordinate coordinate = centerPt.getCoordinate();
+
+        collection.createIndex(IndexOptions.indexOptions(SPATIAL_INDEX), "location");
+        DocumentCursor cursor = collection.find(where("location").near(coordinate, DIST_20_MM));
+        assertEquals("Near filter should return no results within 20mm distance", emptySet(), cursor.toSet());
+    }
+
+    @Test
+    public void testGeoLibrary_Distance_SanityCheck() {
+        final double EPSILON_10M = 10.0;
+        // All of these distances should be within 10 metres of what was estimated on Google Earth.
+
+        GeodesicData s12 = WGS84.Inverse(centerPt.getX(), centerPt.getY(), pt2_950m_ESE.getX(), pt2_950m_ESE.getY());
+        assertEquals(s12.s12, 950.0, EPSILON_10M);
+
+        GeodesicData s13 = WGS84.Inverse(centerPt.getX(), centerPt.getY(), pt3_930m_W.getX(), pt3_930m_W.getY());
+        assertEquals(s13.s12, 930.0, EPSILON_10M);
+
+        GeodesicData s14 = WGS84.Inverse(centerPt.getX(), centerPt.getY(), pt4_2750m_WSW.getX(), pt4_2750m_WSW.getY());
+        assertEquals(s14.s12, 2750.0, EPSILON_10M);
+    }
+}

nitrite-spatial/src/test/java/org/dizitart/no2/spatial/GeoSpatialTestBase.java

@@ -0,0 +1,120 @@
+/*
+ * Copyright (c) 2017-2021 Nitrite author or authors.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *       http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package org.dizitart.no2.spatial;
+
+import lombok.SneakyThrows;
+import org.dizitart.no2.Nitrite;
+import org.dizitart.no2.collection.Document;
+import org.dizitart.no2.collection.NitriteCollection;
+import org.dizitart.no2.repository.ObjectRepository;
+import org.junit.After;
+import org.junit.Before;
+import org.junit.Rule;
+import org.locationtech.jts.geom.Point;
+import org.locationtech.jts.io.ParseException;
+import org.locationtech.jts.io.WKTReader;
+
+import static org.dizitart.no2.collection.Document.createDocument;
+import static org.dizitart.no2.common.Constants.DOC_ID;
+import static org.dizitart.no2.common.Constants.DOC_MODIFIED;
+import static org.dizitart.no2.common.Constants.DOC_REVISION;
+import static org.dizitart.no2.common.Constants.DOC_SOURCE;
+import static org.dizitart.no2.spatial.TestUtil.createDb;
+import static org.dizitart.no2.spatial.TestUtil.deleteDb;
+import static org.dizitart.no2.spatial.TestUtil.getRandomTempDbFile;
+
+public abstract class GeoSpatialTestBase {
+    private String fileName;
+    protected Nitrite db;
+    protected NitriteCollection collection;
+    protected ObjectRepository<SpatialData> repository;
+
+    static WKTReader reader = new WKTReader();
+    protected SpatialData obj2_950m_ESE, obj3_930m_W, obj4_2750m_WSW;
+    protected Document doc2_950m_ESE, doc3_930m_W, doc4_2750m_WSW;
+
+    @Rule
+    public Retry retry = new Retry(3);
+
+    /**
+     * I chose a set of simple landmarks in a major city at high latitude, near 60°N,
+     * such that the separation between them is primarily east-west.
+     * <p>
+     * At the equator, 1 degree of either latitude or longitude measures approx. 111km wide.
+     * However, at 60°N, 1 degree of longitude is only half as wide. (cf. cos(60°) == 0.5)
+     * <p>
+     * This math is not exact enough for the needs of a geographer, but it's close enough to create
+     * simple test cases that can distinguish whether we are properly converting meters to/from degrees,
+     * including accounting for the curvature of the Earth.
+     */
+    public static class TestLocations {
+        static Point centerPt = readPoint("POINT(59.91437 10.73402)");   // National Theater (Oslo)
+        static Point pt2_950m_ESE = readPoint("POINT(59.9115306 10.7501574)"); // Olso Central Station
+        static Point pt3_930m_W = readPoint("POINT(59.91433 10.71730)");  // National Library of Norway
+        static Point pt4_2750m_WSW = readPoint("POINT(59.90749 10.68670)");  //  Norwegian Museum of Cultural Hist.
+    }
+
+    @SneakyThrows
+    protected static Point readPoint(String wkt) {
+        return (Point) reader.read(wkt);
+    }
+
+    @Before
+    public void before() throws ParseException {
+        fileName = getRandomTempDbFile();
+        db = createDb(fileName);
+
+        collection = db.getCollection("test");
+        repository = db.getRepository(SpatialData.class);
+
+        insertObjects();
+        insertDocuments();
+    }
+
+    protected void insertObjects() throws ParseException {
+        obj2_950m_ESE = new SpatialData(2L, TestLocations.pt2_950m_ESE);
+        obj3_930m_W = new SpatialData(3L, TestLocations.pt3_930m_W);
+        obj4_2750m_WSW = new SpatialData(4L, TestLocations.pt4_2750m_WSW);
+        repository.insert(obj2_950m_ESE, obj3_930m_W, obj4_2750m_WSW);
+    }
+
+    protected void insertDocuments() throws ParseException {
+        doc2_950m_ESE = createDocument("key", 2L).put("location", TestLocations.pt2_950m_ESE);
+        doc3_930m_W = createDocument("key", 3L).put("location", TestLocations.pt3_930m_W);
+        doc4_2750m_WSW = createDocument("key", 4L).put("location", TestLocations.pt4_2750m_WSW);
+
+        collection.insert(doc2_950m_ESE, doc3_930m_W, doc4_2750m_WSW);
+    }
+
+    @After
+    public void after() {
+        if (db != null && !db.isClosed()) {
+            db.close();
+        }
+
+        deleteDb(fileName);
+    }
+
+    protected Document trimMeta(Document document) {
+        document.remove(DOC_ID);
+        document.remove(DOC_REVISION);
+        document.remove(DOC_MODIFIED);
+        document.remove(DOC_SOURCE);
+        return document;
+    }
+}