Two arities: + *
Haversine has roughly 0.5% error against full WGS84 ellipsoidal distance + * but is faster and avoids pulling more of GeoTools' referencing stack into + * the hot path. Users who need geodetic precision should transform to a + * projected CRS via {@code spatial:transform} before calling. + */ +public class FunSpatialDistance extends BasicFunction implements IndexUseReporter { + + private static final Logger LOG = LogManager.getLogger(FunSpatialDistance.class); + + private static final FunctionParameterSequenceType GEOMETRY_1 = + new FunctionParameterSequenceType("geometry1", Type.NODE, Cardinality.ZERO_OR_ONE, + "The first geometry"); + private static final FunctionParameterSequenceType GEOMETRY_2 = + new FunctionParameterSequenceType("geometry2", Type.NODE, Cardinality.ZERO_OR_ONE, + "The second geometry"); + private static final FunctionParameterSequenceType UNIT = + new FunctionParameterSequenceType("unit", Type.STRING, Cardinality.EXACTLY_ONE, + """ + Distance unit: "degree" (cartesian, source SRS), "meter", \ + "kilometer", "mile", or "nautical-mile"\ + """); + + /** Mean Earth radius (WGS84 authalic radius), in meters. Used for haversine. */ + private static final double EARTH_RADIUS_METERS = 6_371_007.2d; + + public static final FunctionSignature[] signatures = { + new FunctionSignature( + new QName("distance", SpatialModule.NAMESPACE_URI, SpatialModule.PREFIX), + """ + Returns the cartesian distance between $geometry1 and $geometry2 in \ + the units of their source spatial reference system. For lat/lon \ + geometries (EPSG:4326) this returns degrees; use the 3-arity form \ + with $unit = "meter" for great-circle distance.\ + """, + new SequenceType[]{GEOMETRY_1, GEOMETRY_2}, + new FunctionReturnSequenceType(Type.DOUBLE, Cardinality.ZERO_OR_ONE, + "the cartesian distance, or empty sequence if either argument is empty") + ), + new FunctionSignature( + new QName("distance", SpatialModule.NAMESPACE_URI, SpatialModule.PREFIX), + """ + Returns the distance between $geometry1 and $geometry2 in the \ + requested $unit. Supported units: "degree" (cartesian, source SRS), \ + "meter", "kilometer", "mile", "nautical-mile". Non-degree units \ + transform both geometries to EPSG:4326 and use haversine on the \ + closest-point pair (~0.5% error vs ellipsoidal).\ + """, + new SequenceType[]{GEOMETRY_1, GEOMETRY_2, UNIT}, + new FunctionReturnSequenceType(Type.DOUBLE, Cardinality.ZERO_OR_ONE, + """ + the distance in the requested unit, or empty sequence if \ + either geometry argument is empty\ + """) + ) + }; + + private boolean hasUsedIndex = false; + + public FunSpatialDistance(final XQueryContext context, final FunctionSignature signature) { + super(context, signature); + } + + @Override + public Sequence eval(final Sequence[] args, final Sequence contextSequence) throws XPathException { + if (args[0].isEmpty() || args[1].isEmpty()) { + return Sequence.EMPTY_SEQUENCE; + } + + final String unit = args.length == 3 + ? args[2].itemAt(0).getStringValue().trim().toLowerCase() + : "degree"; + if (!isSupportedUnit(unit)) { + throw new XPathException(this, ErrorCodes.FOER0000, """ + Unsupported distance unit '%s'. Expected one of: \ + degree, meter, kilometer, mile, nautical-mile.\ + """.formatted(unit)); + } + + final AbstractGMLJDBCIndexWorker indexWorker = (AbstractGMLJDBCIndexWorker) + context.getBroker().getIndexController().getWorkerByIndexId(AbstractGMLJDBCIndex.ID); + if (indexWorker == null) { + throw new XPathException(this, ErrorCodes.FOER0000, + "Unable to find a spatial index worker."); + } + + final boolean wantEpsg4326 = !"degree".equals(unit); + try { + final Geometry g1 = resolveGeometry(indexWorker, (NodeValue) args[0].itemAt(0), wantEpsg4326); + final Geometry g2 = resolveGeometry(indexWorker, (NodeValue) args[1].itemAt(0), wantEpsg4326); + if (g1 == null || g2 == null) { + throw new XPathException(this, ErrorCodes.FOER0000, + "Unable to resolve a geometry from the input nodes."); + } + final double result = "degree".equals(unit) + ? g1.distance(g2) + : haversineDistance(g1, g2, unit); + return new DoubleValue(this, result); + } catch (final SpatialIndexException e) { + LOG.error(e.getMessage(), e); + throw new XPathException(this, ErrorCodes.FOER0000, e); + } + } + + private Geometry resolveGeometry(final AbstractGMLJDBCIndexWorker indexWorker, + final NodeValue node, final boolean wantEpsg4326) + throws SpatialIndexException { + if (node.getImplementationType() == NodeValue.PERSISTENT_NODE) { + final Geometry indexed = indexWorker.getGeometryForNode( + context.getBroker(), (NodeProxy) node, wantEpsg4326); + if (indexed != null) { + hasUsedIndex = true; + return indexed; + } + } + // Fallback for in-memory geometries or unindexed persistent nodes. + final Geometry streamed = indexWorker.streamNodeToGeometry(context, node); + if (streamed == null) { + return null; + } + if (!wantEpsg4326) { + return streamed; + } + final String sourceCrs = ((Element) node.getNode()).getAttribute("srsName").trim(); + if (sourceCrs.isEmpty() || "EPSG:4326".equalsIgnoreCase(sourceCrs)) { + return streamed; + } + return indexWorker.transformGeometry(streamed, sourceCrs, "EPSG:4326"); + } + + /** + * Closest-point haversine. Uses JTS DistanceOp to find the nearest pair of + * coordinates on the two (already EPSG:4326-transformed) geometries, then + * computes great-circle distance via haversine and converts to the + * requested unit. + */ + private static double haversineDistance(final Geometry g1, final Geometry g2, final String unit) { + final Coordinate[] nearest = DistanceOp.nearestPoints(g1, g2); + final double meters = haversineMeters(nearest[0].y, nearest[0].x, nearest[1].y, nearest[1].x); + return switch (unit) { + case "meter" -> meters; + case "kilometer" -> meters / 1_000d; + case "mile" -> meters / 1_609.344d; + case "nautical-mile" -> meters / 1_852d; + // "degree" is handled by the caller via Geometry.distance(); not reached here. + default -> throw new IllegalArgumentException("Unsupported unit: " + unit); + }; + } + + private static double haversineMeters(final double lat1Deg, final double lon1Deg, + final double lat2Deg, final double lon2Deg) { + final double lat1 = Math.toRadians(lat1Deg); + final double lat2 = Math.toRadians(lat2Deg); + final double dLat = lat2 - lat1; + final double dLon = Math.toRadians(lon2Deg - lon1Deg); + final double sinDLat = Math.sin(dLat / 2d); + final double sinDLon = Math.sin(dLon / 2d); + final double a = sinDLat * sinDLat + + Math.cos(lat1) * Math.cos(lat2) * sinDLon * sinDLon; + final double c = 2d * Math.atan2(Math.sqrt(a), Math.sqrt(1d - a)); + return EARTH_RADIUS_METERS * c; + } + + private static boolean isSupportedUnit(final String unit) { + return switch (unit) { + case "degree", "meter", "kilometer", "mile", "nautical-mile" -> true; + default -> false; + }; + } + + @Override + public boolean hasUsedIndex() { + return hasUsedIndex; + } +} diff --git a/extensions/indexes/spatial/src/main/java/org/exist/xquery/modules/spatial/SpatialModule.java b/extensions/indexes/spatial/src/main/java/org/exist/xquery/modules/spatial/SpatialModule.java index 9a4869c094d..0b195068dbd 100644 --- a/extensions/indexes/spatial/src/main/java/org/exist/xquery/modules/spatial/SpatialModule.java +++ b/extensions/indexes/spatial/src/main/java/org/exist/xquery/modules/spatial/SpatialModule.java @@ -39,6 +39,10 @@ public class SpatialModule extends AbstractInternalModule { public final static String RELEASED_IN_VERSION = "eXist-1.2"; public static final FunctionDef[] functions = { + // --- Distance & proximity --- + new FunctionDef(FunSpatialDistance.signatures[0], FunSpatialDistance.class), + new FunctionDef(FunSpatialDistance.signatures[1], FunSpatialDistance.class), + // --- End Distance & proximity --- new FunctionDef(FunSpatialSearch.signatures[0], FunSpatialSearch.class), new FunctionDef(FunSpatialSearch.signatures[1], FunSpatialSearch.class), new FunctionDef(FunSpatialSearch.signatures[2], FunSpatialSearch.class), diff --git a/extensions/indexes/spatial/src/test/java/org/exist/indexing/spatial/FunSpatialDistanceTest.java b/extensions/indexes/spatial/src/test/java/org/exist/indexing/spatial/FunSpatialDistanceTest.java new file mode 100644 index 00000000000..6ad5b0f4a93 --- /dev/null +++ b/extensions/indexes/spatial/src/test/java/org/exist/indexing/spatial/FunSpatialDistanceTest.java @@ -0,0 +1,190 @@ +/* + * eXist-db Open Source Native XML Database + * Copyright (C) 2001 The eXist-db Authors + * + * info@exist-db.org + * http://www.exist-db.org + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +package org.exist.indexing.spatial; + +import org.exist.EXistException; +import org.exist.security.PermissionDeniedException; +import org.exist.storage.BrokerPool; +import org.exist.storage.DBBroker; +import org.exist.test.ExistEmbeddedServer; +import org.exist.xquery.XPathException; +import org.exist.xquery.XQuery; +import org.exist.xquery.value.DoubleValue; +import org.exist.xquery.value.Sequence; +import org.junit.ClassRule; +import org.junit.Test; + +import java.util.Optional; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertNotNull; +import static org.junit.Assert.assertTrue; +import static org.junit.Assert.fail; + +/** + * JUnit tests for {@link org.exist.xquery.modules.spatial.FunSpatialDistance}. + * + *
Uses in-memory GML geometries so the tests do not depend on the OS MasterMap
+ * fixture downloaded by the build's wget plugin. The spatial index is still
+ * configured (the embedded server picks up the test {@code conf.xml}); the
+ * geometries here flow through {@code streamNodeToGeometry} rather than
+ * {@code getGeometryForNode}, exercising the in-memory branch of
+ * {@link org.exist.xquery.modules.spatial.FunSpatialDistance#resolveGeometry}.
+ */
+public class FunSpatialDistanceTest {
+
+ @ClassRule
+ public static final ExistEmbeddedServer server = new ExistEmbeddedServer(true, true);
+
+ private static final String SPATIAL_PROLOG = """
+ import module namespace spatial='http://exist-db.org/xquery/spatial'
+ at 'java:org.exist.xquery.modules.spatial.SpatialModule';
+ declare namespace gml='http://www.opengis.net/gml';
+ """;
+
+ /** Times Square, NYC. */
+ private static final String POINT_NYC = """
+