From b410a6325ebcbff460f4df395f18368a92946bf9 Mon Sep 17 00:00:00 2001 From: Rafael Biehler Date: Sun, 20 Oct 2024 23:33:22 +0200 Subject: [PATCH] Use same logic as Neo4j to calculate distances with 3 dimensions --- .../db/functions/CypherFunctions.java | 31 ++++++++++++++++++- .../db/cypher/CypherGeoFunctionsTest.java | 21 +++++++++++-- 2 files changed, 48 insertions(+), 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 index 82e3bce487..07fcb6ea5f 100644 --- 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: GitHub + */ + 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 index 1a1db10437..eff0d3424c 100644 --- 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}),