Use same logic as Neo4j to calculate distances with 3 dimensions

polypheny · Oct 20, 2024 · b410a63 · b410a63
1 parent ee28cd6
commit b410a63
Show file tree

Hide file tree

Showing 2 changed files with 48 additions and 4 deletions.
diff --git a/core/src/main/java/org/polypheny/db/functions/CypherFunctions.java b/core/src/main/java/org/polypheny/db/functions/CypherFunctions.java
@@ -49,6 +49,8 @@
 import org.polypheny.db.type.entity.spatial.PolyGeometry;
 import org.polypheny.db.type.entity.spatial.PolyPoint;
 
+import static java.lang.Math.toRadians;
+import static org.polypheny.db.functions.spatial.GeoDistanceFunctions.EARTH_RADIUS_M;
 import static org.polypheny.db.type.entity.spatial.PolyGeometry.WGS_84;
 import static org.polypheny.db.type.entity.spatial.PolyGeometry.WGS_84_3D;
 
@@ -565,12 +567,39 @@ public static PolyDouble distance( PolyValue p1, PolyValue p2 ) {
                     return new PolyDouble(
                             Math.sqrt( Math.pow( g2.getX() - g1.getX(), 2 ) + Math.pow( g2.getY() - g1.getY(), 2 ) + Math.pow( g2.getZ() - g1.getZ(), 2 ) )
                     );
+                } else if ( srid == WGS_84_3D ) {
+                    // See https://github.com/neo4j/neo4j/blob/5.20/community/values/src/main/java/org/neo4j/values/storable/CRSCalculator.java
+                    double greatCircleDistance = getGreatCircleDistance( g1, g2 );
+                    double avgHeight = (g1.getZ() + g2.getZ()) / 2;
+                    // Note: Neo4j uses a different earth radius of 6378140.0, which is why the same calculation
+                    //       in Neo4j does not yield (exactly) the same results.
+                    double distance2D = (EARTH_RADIUS_M + avgHeight) * greatCircleDistance;
+                    double heightDifference = g1.getZ() - g2.getZ();
+                    return new PolyDouble( Math.sqrt( Math.pow( distance2D, 2 ) + Math.pow( heightDifference, 2 ) ) );
                 }
+            } else {
+                return new PolyDouble( g1.distance( g2 ) );
             }
-            return new PolyDouble( g1.distance( g2 ) );
         } catch ( GeometryTopologicalException e ) {
             throw new GenericRuntimeException( e );
         }
+
+        throw new GenericRuntimeException( "This should not be possible!" );
+    }
+
+
+    /**
+     * Use same logic as Neo4j to calculate the spherical distance.
+     * See: <a href="https://github.com/neo4j/neo4j/blob/5.20/community/values/src/main/java/org/neo4j/values/storable/CRSCalculator.java">GitHub</a>
+     */
+    private static double getGreatCircleDistance( PolyPoint g1, PolyPoint g2 ) {
+        double lat1 = toRadians( g1.getY() );
+        double lat2 = toRadians( g2.getY() );
+        double latDifference = lat2 - lat1;
+        double lonDifference = toRadians( g2.getX() - g1.getX() );
+        double alpha = Math.pow( Math.sin( latDifference / 2 ), 2 ) +
+                Math.cos( lat1 ) * Math.cos( lat2 ) * Math.pow( Math.sin( lonDifference / 2 ), 2 );
+        return 2.0 * Math.atan2( Math.sqrt( alpha ), Math.sqrt( 1 - alpha ) );
     }
 
 

diff --git a/dbms/src/test/java/org/polypheny/db/cypher/CypherGeoFunctionsTest.java b/dbms/src/test/java/org/polypheny/db/cypher/CypherGeoFunctionsTest.java
@@ -79,7 +79,7 @@ public void createPointFromNodeFields() {
 
     @Test
     public void distanceTest() {
-        // Compute distance in spherical coordinate system
+        // Compute distance in spherical coordinate system (2 dimensions)
         execute( """
                 CREATE (basel:City {name: 'Basel', latitude: 47.5595, longitude: 7.5885}),
                        (zurich:City {name: 'Zürich', latitude: 47.3770, longitude: 8.5416});
@@ -94,7 +94,22 @@ public void distanceTest() {
         assert res.data[0].length == 3;
         assert Math.abs( PolyValue.fromJson( res.data[0][2] ).asDocument().get( new PolyString( "value" ) ).asDouble().doubleValue() - 74460.31287583392 ) < 1e-9;
 
-        // Compute distance in euclidean coordinate system
+        // Compute distance in spherical coordinate system (3 dimensions)
+        execute( """
+                CREATE (basel:City {name: 'Basel', latitude: 47.5595, longitude: 7.5885}),
+                       (zurich:City {name: 'Zürich', latitude: 47.3770, longitude: 8.5416});
+                """ );
+        res = execute( """
+                MATCH (basel:City {name: 'Basel'}), (zurich:City {name: 'Zürich'})
+                WITH basel, zurich,
+                     point({latitude: basel.latitude, longitude: basel.longitude, height: 100}) AS point1,
+                     point({latitude: zurich.latitude, longitude: zurich.longitude, height: 200}) AS point2
+                RETURN basel.name, zurich.name, point.distance(point1, point2) AS distance;
+                """ );
+        assert res.data[0].length == 3;
+        assert Math.abs( PolyValue.fromJson( res.data[0][2] ).asDocument().get( new PolyString( "value" ) ).asDouble().doubleValue() - 74462.13313143898 ) < 1e-9;
+
+        // Compute distance in euclidean coordinate system (2 dimensions)
         execute( """
                 CREATE (a:Dot {x: 1, y: 1}),
                        (b:Dot {x: 2, y: 2}),
@@ -110,7 +125,7 @@ public void distanceTest() {
         assert res.data[0].length == 1;
         assert Math.abs( PolyValue.fromJson( res.data[0][0] ).asDocument().get( new PolyString( "value" ) ).asDouble().doubleValue() - Math.sqrt( 2 )) < 1e-9;
 
-        // Compute distance in euclidean coordinate system with 3 coordinates
+        // Compute distance in euclidean coordinate system (3 dimensions)
         execute( """
                 CREATE (a:Dot3D {x: 1, y: 1, z:1}),
                        (b:Dot3D {x: 2, y: 2, z:2}),