diff --git a/common/src/main/java/org/apache/sedona/common/Functions.java b/common/src/main/java/org/apache/sedona/common/Functions.java
index e0932a01d45..fa8a7da9e25 100644
--- a/common/src/main/java/org/apache/sedona/common/Functions.java
+++ b/common/src/main/java/org/apache/sedona/common/Functions.java
@@ -1362,6 +1362,10 @@ public static Geometry intersection(Geometry leftGeometry, Geometry rightGeometr
return leftGeometry.intersection(rightGeometry);
}
+ public static Geometry sharedPaths(Geometry leftGeometry, Geometry rightGeometry) {
+ return SharedPaths.compute(leftGeometry, rightGeometry);
+ }
+
public static Geometry makeValid(Geometry geometry, boolean keepCollapsed) {
GeometryFixer fixer = new GeometryFixer(geometry);
fixer.setKeepCollapsed(keepCollapsed);
diff --git a/common/src/main/java/org/apache/sedona/common/utils/GeomUtils.java b/common/src/main/java/org/apache/sedona/common/utils/GeomUtils.java
index eec89f9d8a8..7648c87a6ac 100644
--- a/common/src/main/java/org/apache/sedona/common/utils/GeomUtils.java
+++ b/common/src/main/java/org/apache/sedona/common/utils/GeomUtils.java
@@ -175,14 +175,23 @@ public static String getEWKT(Geometry geometry) {
if (srid != 0) {
sridString = "SRID=" + String.valueOf(srid) + ";";
}
- return sridString + new WKTWriter(4).write(geometry);
+ return sridString + writeWKT(geometry);
}
public static String getWKT(Geometry geometry) {
if (geometry == null) {
return null;
}
- return new WKTWriter(4).write(geometry);
+ return writeWKT(geometry);
+ }
+
+ private static String writeWKT(Geometry geometry) {
+ // JTS 1.20 omits the separator between a dimensional marker and EMPTY for nested empty
+ // geometries (for example, "MULTILINESTRING ZEMPTY"), which its own WKTReader rejects.
+ String wkt = new WKTWriter(4).write(geometry);
+ return wkt.replace("ZMEMPTY", "ZM EMPTY")
+ .replace("ZEMPTY", "Z EMPTY")
+ .replace("MEMPTY", "M EMPTY");
}
public static String getHexEWKB(Geometry geometry, int endian) {
diff --git a/common/src/main/java/org/apache/sedona/common/utils/SharedPaths.java b/common/src/main/java/org/apache/sedona/common/utils/SharedPaths.java
new file mode 100644
index 00000000000..d577a96ab9e
--- /dev/null
+++ b/common/src/main/java/org/apache/sedona/common/utils/SharedPaths.java
@@ -0,0 +1,764 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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.apache.sedona.common.utils;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.locationtech.jts.algorithm.Distance;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateSequence;
+import org.locationtech.jts.geom.CoordinateSequenceFactory;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryCollection;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LineString;
+import org.locationtech.jts.geom.MultiLineString;
+import org.locationtech.jts.geom.PrecisionModel;
+import org.locationtech.jts.geom.util.LineStringExtracter;
+import org.locationtech.jts.index.strtree.ItemDistance;
+import org.locationtech.jts.index.strtree.STRtree;
+import org.locationtech.jts.linearref.LinearLocation;
+import org.locationtech.jts.operation.overlayng.OverlayNG;
+import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
+
+/**
+ * Computes the portions of two lineal geometries which follow the same or opposite direction.
+ *
+ * For valid, non-null inputs, the result always has this shape:
+ *
+ *
+ * GeometryCollection
+ * +-- element 0: MultiLineString of same-direction paths
+ * `-- element 1: MultiLineString of opposite-direction paths
+ *
+ *
+ * Every returned path is oriented like the first input. For example:
+ *
+ *
+ * left : A --------------------> B
+ * right: A --------------------> B => element 0 (same)
+ * right: A <-------------------- B => element 1 (opposite)
+ *
+ *
+ * The implementation follows these stages:
+ *
+ *
+ * validate inputs
+ * |
+ * Stage 1: robust 2D intersection and shared-line extraction
+ * |
+ * Stage 2: build an index over each source geometry
+ * |
+ * Stage 3: locate and compare each fragment in both sources
+ * |
+ * Stage 4: orient like the first input and restore source-derived Z values
+ * |
+ * Stage 5: assemble the same-direction and opposite-direction buckets
+ *
+ *
+ * The spatial indices are important: an overlay can emit one fragment per source segment, so
+ * rescanning both inputs for every fragment would make identical long lines quadratic.
+ */
+public final class SharedPaths {
+ // Interior samples avoid the endpoint ambiguity of closed lines and shared source vertices.
+ // The same pair of fractions is used consistently when classifying traversal direction and
+ // when disambiguating the source segment used for Z interpolation.
+ private static final double START_SAMPLE_FRACTION = 0.1;
+ private static final double END_SAMPLE_FRACTION = 0.9;
+
+ private SharedPaths() {}
+
+ /**
+ * Returns a collection containing same-direction paths followed by opposite-direction paths.
+ * Paths are oriented in the direction of {@code left}. Matching SRIDs are required. When either
+ * input has Z, the result retains Z only if every returned coordinate resolves to a finite source
+ * Z; otherwise the whole result is XY. Like PostGIS, the result does not retain M.
+ */
+ public static Geometry compute(Geometry left, Geometry right) {
+ // Match null-propagating SQL behavior before attempting type or SRID validation.
+ if (left == null || right == null) {
+ return null;
+ }
+ validateSRIDs(left, right);
+ validateLineal(left);
+ validateLineal(right);
+
+ // The first input owns the result's precision model, coordinate sequence implementation, and
+ // (already validated) SRID. This also keeps every nested result geometry on the same factory.
+ GeometryFactory resultFactory = resultFactory(left);
+ if (left.isEmpty() || right.isEmpty()) {
+ return emptyResult(resultFactory);
+ }
+
+ // A shared path requires the envelopes to overlap at overlay precision. Avoid overlay setup,
+ // dimensional scanning, and result assembly when the precision-adjusted envelopes prove the
+ // result is empty.
+ if (envelopesDisjoint(left, right)) {
+ return emptyResult(resultFactory);
+ }
+
+ // Stage 1: robust overlay finds the physical shared coverage. Point-only intersections are
+ // intentionally ignored later because LineStringExtracter returns only linear components.
+ Geometry intersection = OverlayNGRobust.overlay(left, right, OverlayNG.INTERSECTION);
+ boolean retainZ = hasZ(left) || hasZ(right);
+
+ @SuppressWarnings("unchecked")
+ List paths = LineStringExtracter.getLines(intersection);
+
+ // Stage 2: build each source index once. Besides avoiding repeated scans, the index preserves
+ // source traversal order so repeated or self-overlapping paths have deterministic semantics.
+ double overlayTolerance = overlayTolerance(left.getPrecisionModel());
+ SourceSegmentIndex leftSegments =
+ !paths.isEmpty() ? new SourceSegmentIndex(left, overlayTolerance) : null;
+ SourceSegmentIndex rightSegments =
+ !paths.isEmpty() ? new SourceSegmentIndex(right, overlayTolerance) : null;
+ List preparedPaths = new ArrayList<>();
+ for (LineString path : paths) {
+ if (path.isEmpty()) {
+ continue;
+ }
+ // Stage 3: determine how this overlay fragment is traversed in each original input.
+ Direction direction = direction(path, leftSegments, rightSegments);
+
+ // Stage 4: overlay output direction is not a contract. Normalize it to the direction of the
+ // first input before reconstructing ordinates or adding the path to a result bucket.
+ LineString orientedPath = direction.forwardOnLeft ? path : (LineString) path.reverse();
+ double[] zValues =
+ retainZ ? resolvePostGISZ(orientedPath, leftSegments, rightSegments) : null;
+ if (retainZ && zValues == null) {
+ // WKT and WKB assign one dimensionality to the whole result. If any returned coordinate
+ // has no finite source Z, use XY throughout rather than fabricating a value or emitting
+ // NaN.
+ retainZ = false;
+ }
+ preparedPaths.add(new PreparedPath(orientedPath, direction.same, zValues));
+ }
+
+ // Stage 5: even when one bucket is empty, retain both typed MultiLineString children.
+ List sameDirection = new ArrayList<>();
+ List oppositeDirection = new ArrayList<>();
+ for (PreparedPath preparedPath : preparedPaths) {
+ LineString resultPath =
+ copyWithOrdinates(
+ preparedPath.path, resultFactory, retainZ ? preparedPath.zValues : null);
+ (preparedPath.same ? sameDirection : oppositeDirection).add(resultPath);
+ }
+ MultiLineString same =
+ resultFactory.createMultiLineString(sameDirection.toArray(new LineString[0]));
+ MultiLineString opposite =
+ resultFactory.createMultiLineString(oppositeDirection.toArray(new LineString[0]));
+ return resultFactory.createGeometryCollection(new Geometry[] {same, opposite});
+ }
+
+ private static void validateLineal(Geometry geometry) {
+ if (!(geometry instanceof LineString) && !(geometry instanceof MultiLineString)) {
+ throw new IllegalArgumentException("Geometry is not lineal");
+ }
+ }
+
+ private static void validateSRIDs(Geometry left, Geometry right) {
+ if (left.getSRID() != right.getSRID()) {
+ throw new IllegalArgumentException(
+ String.format(
+ "Operation on mixed SRID geometries (%d != %d)", left.getSRID(), right.getSRID()));
+ }
+ }
+
+ private static GeometryFactory resultFactory(Geometry geometry) {
+ return new GeometryFactory(
+ new PrecisionModel(geometry.getPrecisionModel()),
+ geometry.getSRID(),
+ geometry.getFactory().getCoordinateSequenceFactory());
+ }
+
+ private static GeometryCollection emptyResult(GeometryFactory factory) {
+ // PostGIS returns two empty MultiLineStrings rather than an empty GeometryCollection or null.
+ MultiLineString same = factory.createMultiLineString();
+ MultiLineString opposite = factory.createMultiLineString();
+ return factory.createGeometryCollection(new Geometry[] {same, opposite});
+ }
+
+ private static double overlayTolerance(PrecisionModel precisionModel) {
+ if (precisionModel.isFloating()) {
+ return 0;
+ }
+ // A fixed-precision overlay may round both ordinates by half a grid cell. Use the diagonal
+ // displacement so the snapped result can still be located on the original source linework.
+ return Math.sqrt(0.5) * precisionModel.gridSize();
+ }
+
+ private static boolean envelopesDisjoint(Geometry left, Geometry right) {
+ Envelope leftEnvelope = left.getEnvelopeInternal();
+ Envelope rightEnvelope = right.getEnvelopeInternal();
+ PrecisionModel precisionModel = left.getPrecisionModel();
+ if (precisionModel.isFloating()) {
+ return leftEnvelope.disjoint(rightEnvelope);
+ }
+
+ // Fixed-precision overlay rounds coordinates to the first input's grid. Compare rounded
+ // envelope bounds so nearby raw envelopes are not rejected when rounding makes them overlap.
+ return precisionModel.makePrecise(rightEnvelope.getMinX())
+ > precisionModel.makePrecise(leftEnvelope.getMaxX())
+ || precisionModel.makePrecise(rightEnvelope.getMaxX())
+ < precisionModel.makePrecise(leftEnvelope.getMinX())
+ || precisionModel.makePrecise(rightEnvelope.getMinY())
+ > precisionModel.makePrecise(leftEnvelope.getMaxY())
+ || precisionModel.makePrecise(rightEnvelope.getMaxY())
+ < precisionModel.makePrecise(leftEnvelope.getMinY());
+ }
+
+ private static Direction direction(
+ LineString path, SourceSegmentIndex leftSegments, SourceSegmentIndex rightSegments) {
+ CoordinateSequence sequence = path.getCoordinateSequence();
+
+ // Each boolean says whether the overlay fragment's current coordinate order is forward in the
+ // corresponding source. Equal booleans mean both sources traverse the fragment the same way.
+ boolean forwardOnLeft = leftSegments.isForward(sequence);
+ boolean forwardOnRight = rightSegments.isForward(sequence);
+ return new Direction(forwardOnLeft, forwardOnLeft == forwardOnRight);
+ }
+
+ private static int firstNonZeroSegment(CoordinateSequence sequence) {
+ // Overlay output should be non-degenerate, but repeated coordinates are legal in a line. Skip
+ // them so the direction samples always span a segment with a meaningful orientation.
+ for (int i = 0; i < sequence.size() - 1; i++) {
+ if (sequence.getX(i) != sequence.getX(i + 1) || sequence.getY(i) != sequence.getY(i + 1)) {
+ return i;
+ }
+ }
+ throw new IllegalArgumentException("Shared path does not contain a non-zero segment");
+ }
+
+ private static LineString copyWithOrdinates(
+ LineString path, GeometryFactory factory, double[] zValues) {
+ CoordinateSequence source = path.getCoordinateSequence();
+ CoordinateSequenceFactory sequenceFactory = factory.getCoordinateSequenceFactory();
+
+ // Shared-path topology is two-dimensional. Build a fresh XY or XYZ sequence so M is always
+ // removed. A non-null zValues array has already been checked to contain only finite values.
+ CoordinateSequence target = sequenceFactory.create(source.size(), zValues != null ? 3 : 2, 0);
+
+ for (int i = 0; i < source.size(); i++) {
+ double x = source.getX(i);
+ double y = source.getY(i);
+ target.setOrdinate(i, CoordinateSequence.X, x);
+ target.setOrdinate(i, CoordinateSequence.Y, y);
+ if (zValues != null) {
+ target.setOrdinate(i, CoordinateSequence.Z, zValues[i]);
+ }
+ }
+ return factory.createLineString(target);
+ }
+
+ private static double[] resolvePostGISZ(
+ LineString path, SourceSegmentIndex leftSegments, SourceSegmentIndex rightSegments) {
+ CoordinateSequence sequence = path.getCoordinateSequence();
+ double[] zValues = new double[sequence.size()];
+ for (int i = 0; i < sequence.size(); i++) {
+ zValues[i] = postGISZ(sequence, i, leftSegments, rightSegments);
+ if (!Double.isFinite(zValues[i])) {
+ return null;
+ }
+ }
+ return zValues;
+ }
+
+ private static double postGISZ(
+ CoordinateSequence path,
+ int coordinateIndex,
+ SourceSegmentIndex leftSegments,
+ SourceSegmentIndex rightSegments) {
+ // Z precedence captures the stable PostGIS/GEOS behavior used by shared paths:
+ //
+ // exact vertex in left
+ // -> exact vertex in right
+ // -> interpolation on the matching left traversal
+ // -> interpolation on the matching right traversal
+ //
+ // Exact vertices are indexed globally because a noded overlay coordinate can coincide with a
+ // vertex on another traversal of the same non-simple input. Interpolation, by contrast, must
+ // use the traversal containing the shared path or it can pick the wrong Z at a crossing.
+ double x = path.getX(coordinateIndex);
+ double y = path.getY(coordinateIndex);
+ double exactZ = leftSegments.exactZ(x, y);
+ if (Double.isFinite(exactZ)) {
+ return exactZ;
+ }
+ exactZ = rightSegments.exactZ(x, y);
+ if (Double.isFinite(exactZ)) {
+ return exactZ;
+ }
+ // Segment searches are deliberately lazy. Exact vertices are the common case, and a source
+ // with no finite-Z segment can never contribute an interpolated value.
+ if (leftSegments.hasInterpolatableZ()) {
+ double leftZ = leftSegments.findForPathCoordinate(path, coordinateIndex).zAt(x, y);
+ if (Double.isFinite(leftZ)) {
+ return leftZ;
+ }
+ }
+ if (rightSegments.hasInterpolatableZ()) {
+ return rightSegments.findForPathCoordinate(path, coordinateIndex).zAt(x, y);
+ }
+ return Double.NaN;
+ }
+
+ private static boolean hasZ(Geometry geometry) {
+ // LineString and MultiLineString are the only accepted inputs, so their component coordinate
+ // sequences fully describe the dimensionality relevant to the output.
+ for (int component = 0; component < geometry.getNumGeometries(); component++) {
+ Geometry child = geometry.getGeometryN(component);
+ if (child instanceof LineString && ((LineString) child).getCoordinateSequence().hasZ()) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ private static final class Direction {
+ // Whether the overlay fragment's current order agrees with the first input. This controls
+ // whether the fragment must be reversed before it is returned.
+ private final boolean forwardOnLeft;
+
+ // Whether both inputs traverse the physical path in the same direction. This selects the
+ // first or second MultiLineString in the result GeometryCollection.
+ private final boolean same;
+
+ private Direction(boolean forwardOnLeft, boolean same) {
+ this.forwardOnLeft = forwardOnLeft;
+ this.same = same;
+ }
+ }
+
+ private static final class PreparedPath {
+ private final LineString path;
+ private final boolean same;
+ private final double[] zValues;
+
+ private PreparedPath(LineString path, boolean same, double[] zValues) {
+ this.path = path;
+ this.same = same;
+ this.zValues = zValues;
+ }
+ }
+
+ /**
+ * Indexes source locations for path direction and segment-aware Z interpolation.
+ *
+ * Two complementary indices are built together during construction:
+ *
+ *
+ * source vertices -----------------> exactZByVertex (O(1) exact-Z lookup)
+ * source segment envelopes --------> STRtree (spatial candidate lookup)
+ * source traversal order ----------> segment ordinal (stable tie-breaking)
+ *
+ *
+ * A point may belong to several segments at a self-intersection or repeated path. Candidate
+ * selection therefore uses geometric distance first and source traversal order second.
+ */
+ private static final class SourceSegmentIndex {
+ private static final ItemDistance SEGMENT_TO_POINT_DISTANCE =
+ (first, second) -> {
+ Object firstItem = first.getItem();
+ SourceSegment segment =
+ firstItem instanceof SourceSegment
+ ? (SourceSegment) firstItem
+ : (SourceSegment) second.getItem();
+ Coordinate coordinate =
+ firstItem instanceof Coordinate
+ ? (Coordinate) firstItem
+ : (Coordinate) second.getItem();
+ return segment.distanceTo(coordinate);
+ };
+
+ private final STRtree index = new STRtree();
+ private final Map exactZByVertex = new HashMap<>();
+ private boolean hasInterpolatableZ;
+
+ // Overlay interpolation can move a theoretically-on-segment sample by a few ULPs, while a
+ // fixed precision model can move it to the nearest grid point. Track both displacements so the
+ // STRtree query envelope still reaches the true segment before the distance check below.
+ private double queryTolerance;
+
+ private SourceSegmentIndex(Geometry geometry, double overlayTolerance) {
+ queryTolerance = Math.max(16 * Math.ulp(1.0), overlayTolerance);
+ // Ordinals flatten component/segment positions into source traversal order. Zero-length
+ // segments are omitted because no interior sample can be located on them; omitting them does
+ // not change the relative order of usable segments.
+ int ordinal = 0;
+ for (int component = 0; component < geometry.getNumGeometries(); component++) {
+ LineString line = (LineString) geometry.getGeometryN(component);
+ CoordinateSequence sequence = line.getCoordinateSequence();
+
+ // Use a source-aware query expansion at large coordinate magnitudes while retaining a
+ // small absolute floor near zero.
+ double maxCoordinateMagnitude = 1.0;
+ for (int vertex = 0; vertex < sequence.size(); vertex++) {
+ maxCoordinateMagnitude =
+ Math.max(
+ maxCoordinateMagnitude,
+ Math.max(Math.abs(sequence.getX(vertex)), Math.abs(sequence.getY(vertex))));
+ }
+ queryTolerance = Math.max(queryTolerance, 16 * Math.ulp(maxCoordinateMagnitude));
+ if (sequence.hasZ()) {
+ // putIfAbsent preserves the first finite Z encountered along the source traversal. This
+ // makes an exact vertex deterministic when a non-simple line visits the same XY more
+ // than once with different Z values.
+ for (int vertex = 0; vertex < sequence.size(); vertex++) {
+ double z = sequence.getZ(vertex);
+ if (Double.isFinite(z)) {
+ exactZByVertex.putIfAbsent(
+ new XYKey(sequence.getX(vertex), sequence.getY(vertex)), z);
+ }
+ }
+ }
+ for (int segment = 0; segment < sequence.size() - 1; segment++) {
+ Coordinate start = sequence.getCoordinateCopy(segment);
+ Coordinate end = sequence.getCoordinateCopy(segment + 1);
+ if (start.equals2D(end)) {
+ continue;
+ }
+
+ // The envelope is only a candidate filter. Actual distance and traversal order rank the
+ // segments after querying the tree.
+ SourceSegment sourceSegment = new SourceSegment(start, end, sequence.hasZ(), ordinal++);
+ hasInterpolatableZ |= sourceSegment.hasInterpolatableZ();
+ index.insert(new Envelope(start, end), sourceSegment);
+ }
+ }
+ index.build();
+ }
+
+ private double exactZ(double x, double y) {
+ // Exact means bitwise-equal XY after normalizing signed zero; interpolated coordinates fall
+ // through to the traversal-aware segment logic.
+ Double z = exactZByVertex.get(new XYKey(x, y));
+ return z == null ? Double.NaN : z;
+ }
+
+ private boolean hasInterpolatableZ() {
+ return hasInterpolatableZ;
+ }
+
+ private boolean isForward(CoordinateSequence path) {
+ int segmentIndex = firstNonZeroSegment(path);
+ Coordinate start = path.getCoordinateCopy(segmentIndex);
+ Coordinate end = path.getCoordinateCopy(segmentIndex + 1);
+
+ // Sample inside the first non-zero output segment. Sampling away from endpoints avoids
+ // choosing the preceding segment at a source vertex and avoids the coincident start/end of a
+ // closed line. The samples are located independently because an overlay edge may span one or
+ // more source vertices after collinear nodes are collapsed.
+ Coordinate startSample =
+ LinearLocation.pointAlongSegmentByFraction(start, end, START_SAMPLE_FRACTION);
+ Coordinate endSample =
+ LinearLocation.pointAlongSegmentByFraction(start, end, END_SAMPLE_FRACTION);
+ return locate(startSample).compareTo(locate(endSample)) < 0;
+ }
+
+ private SourceLocation locate(Coordinate coordinate) {
+ SourceSegment closest = closestSegment(coordinate, candidates(coordinate), queryTolerance);
+ if (closest == null) {
+ // Robust overlay snapping can move a valid result farther than a predictable number of
+ // ULPs. An indexed nearest-neighbour fallback avoids failing the whole query while keeping
+ // the normal on-source path bounded by the small tolerance above.
+ closest =
+ closestSegment(coordinate, nearestCandidates(coordinate), Double.POSITIVE_INFINITY);
+ }
+ if (closest == null) {
+ throw new IllegalStateException("Shared path point is not present in source geometry");
+ }
+ return new SourceLocation(closest, coordinate);
+ }
+
+ private SourceSegment findForPathCoordinate(CoordinateSequence path, int coordinateIndex) {
+ int adjacentIndex = adjacentNonZeroCoordinate(path, coordinateIndex);
+ Coordinate coordinate = path.getCoordinateCopy(coordinateIndex);
+ Coordinate adjacent = path.getCoordinateCopy(adjacentIndex);
+ SourceSegment closest =
+ closestPathSegment(coordinate, adjacent, candidates(coordinate), true);
+ if (closest == null) {
+ closest = closestPathSegment(coordinate, adjacent, nearestCandidates(coordinate), false);
+ }
+ if (closest == null) {
+ throw new IllegalStateException("Shared path segment is not present in source geometry");
+ }
+ return closest;
+ }
+
+ private SourceSegment closestPathSegment(
+ Coordinate coordinate,
+ Coordinate adjacent,
+ List candidates,
+ boolean requireTolerance) {
+ SourceSegment closest = null;
+ double closestDistance = Double.POSITIVE_INFINITY;
+ for (SourceSegment candidate : candidates) {
+ // Overlay can collapse many collinear source vertices into one long output edge. Limit the
+ // test samples to the portion of that edge covered by this candidate source segment.
+ double sharedFraction = candidate.sharedFractionFrom(coordinate, adjacent);
+ if (!(sharedFraction > 0)) {
+ continue;
+ }
+ Coordinate startSample =
+ LinearLocation.pointAlongSegmentByFraction(
+ coordinate, adjacent, sharedFraction * START_SAMPLE_FRACTION);
+ Coordinate endSample =
+ LinearLocation.pointAlongSegmentByFraction(
+ coordinate, adjacent, sharedFraction * END_SAMPLE_FRACTION);
+ // Requiring two interior samples, rather than one midpoint, distinguishes the intended
+ // traversal when the midpoint itself is a self-intersection with another source segment.
+ double startDistance = candidate.distanceTo(startSample);
+ double endDistance = candidate.distanceTo(endSample);
+ if (requireTolerance && (startDistance > queryTolerance || endDistance > queryTolerance)) {
+ continue;
+ }
+ double distance = Math.max(startDistance, endDistance);
+
+ // As in locate(), distance protects against tolerance-admitted nearby segments and ordinal
+ // makes repeated coincident traversals deterministic.
+ if (isBetter(candidate, distance, closest, closestDistance)) {
+ closest = candidate;
+ closestDistance = distance;
+ }
+ }
+ return closest;
+ }
+
+ private static SourceSegment closestSegment(
+ Coordinate coordinate, List candidates, double maximumDistance) {
+ SourceSegment closest = null;
+ double closestDistance = Double.POSITIVE_INFINITY;
+ for (SourceSegment candidate : candidates) {
+ double distance = candidate.distanceTo(coordinate);
+ if (distance <= maximumDistance
+ && isBetter(candidate, distance, closest, closestDistance)) {
+ closest = candidate;
+ closestDistance = distance;
+ }
+ }
+ return closest;
+ }
+
+ private static boolean isBetter(
+ SourceSegment candidate, double distance, SourceSegment closest, double closestDistance) {
+ int distanceComparison = Double.compare(distance, closestDistance);
+ return closest == null
+ || distanceComparison < 0
+ || (distanceComparison == 0 && candidate.ordinal < closest.ordinal);
+ }
+
+ private List candidates(Coordinate coordinate) {
+ // Expand before querying; otherwise an interpolated or grid-snapped sample just outside a
+ // segment envelope would never reach the tolerance-aware distance check.
+ Envelope searchEnvelope = new Envelope(coordinate);
+ searchEnvelope.expandBy(queryTolerance);
+ @SuppressWarnings("unchecked")
+ List candidates = index.query(searchEnvelope);
+ return candidates;
+ }
+
+ private List nearestCandidates(Coordinate coordinate) {
+ SourceSegment nearest =
+ (SourceSegment)
+ index.nearestNeighbour(
+ new Envelope(coordinate), coordinate, SEGMENT_TO_POINT_DISTANCE);
+ if (nearest == null) {
+ return new ArrayList<>();
+ }
+ double distance = nearest.distanceTo(coordinate);
+ Envelope searchEnvelope = new Envelope(coordinate);
+ searchEnvelope.expandBy(Math.nextUp(distance + queryTolerance));
+ @SuppressWarnings("unchecked")
+ List candidates = index.query(searchEnvelope);
+ return candidates;
+ }
+
+ private static int adjacentNonZeroCoordinate(CoordinateSequence sequence, int coordinateIndex) {
+ // Prefer the following path edge so an interior vertex uses the traversal leaving it. For the
+ // final coordinate, fall back to the preceding edge. Repeated XY coordinates are skipped.
+ double x = sequence.getX(coordinateIndex);
+ double y = sequence.getY(coordinateIndex);
+ for (int i = coordinateIndex + 1; i < sequence.size(); i++) {
+ if (sequence.getX(i) != x || sequence.getY(i) != y) {
+ return i;
+ }
+ }
+ for (int i = coordinateIndex - 1; i >= 0; i--) {
+ if (sequence.getX(i) != x || sequence.getY(i) != y) {
+ return i;
+ }
+ }
+ throw new IllegalArgumentException("Shared path does not contain a non-zero segment");
+ }
+ }
+
+ /** A non-zero source segment plus its position in the flattened source traversal. */
+ private static final class SourceSegment {
+ private final Coordinate start;
+ private final Coordinate end;
+ private final boolean hasZ;
+ private final int ordinal;
+
+ private SourceSegment(Coordinate start, Coordinate end, boolean hasZ, int ordinal) {
+ this.start = start;
+ this.end = end;
+ this.hasZ = hasZ;
+ this.ordinal = ordinal;
+ }
+
+ private double zAt(double x, double y) {
+ boolean atStart = start.getX() == x && start.getY() == y;
+ boolean atEnd = end.getX() == x && end.getY() == y;
+ if (!hasZ) {
+ return Double.NaN;
+ }
+
+ // Preserve exact endpoint ordinates before considering interpolation. The global exact-Z map
+ // normally handles these cases; keeping this local rule also makes the segment primitive
+ // correct when used independently.
+ if (atStart) {
+ return start.getZ();
+ }
+ if (atEnd) {
+ return end.getZ();
+ }
+ double startZ = start.getZ();
+ double endZ = end.getZ();
+ if (!Double.isFinite(startZ) || !Double.isFinite(endZ)) {
+ return Double.NaN;
+ }
+
+ // Interpolate on the dominant XY axis. For an exactly on-segment point the two axes are
+ // mathematically equivalent; the dominant axis reduces roundoff and avoids division by a
+ // tiny delta.
+ double dx = end.getX() - start.getX();
+ double dy = end.getY() - start.getY();
+ double fraction =
+ Math.abs(dx) >= Math.abs(dy) ? (x - start.getX()) / dx : (y - start.getY()) / dy;
+ return startZ + fraction * (endZ - startZ);
+ }
+
+ private boolean hasInterpolatableZ() {
+ return hasZ && Double.isFinite(start.getZ()) && Double.isFinite(end.getZ());
+ }
+
+ private double sharedFractionFrom(Coordinate coordinate, Coordinate adjacent) {
+ // Project both source endpoints onto the outgoing output edge. The largest positive projected
+ // fraction bounds candidate-local samples; the later distance checks reject candidates
+ // which are not actually collinear and overlapping. Clamp to one because only the current
+ // output edge is relevant.
+ double dx = adjacent.getX() - coordinate.getX();
+ double dy = adjacent.getY() - coordinate.getY();
+ double startFraction;
+ double endFraction;
+ if (Math.abs(dx) >= Math.abs(dy)) {
+ startFraction = (start.getX() - coordinate.getX()) / dx;
+ endFraction = (end.getX() - coordinate.getX()) / dx;
+ } else {
+ startFraction = (start.getY() - coordinate.getY()) / dy;
+ endFraction = (end.getY() - coordinate.getY()) / dy;
+ }
+ return Math.min(1.0, Math.max(startFraction, endFraction));
+ }
+
+ private double distanceTo(Coordinate coordinate) {
+ return Distance.pointToSegment(coordinate, start, end);
+ }
+ }
+
+ /** Hash key for an exact XY vertex lookup; Z and M deliberately do not participate. */
+ private static final class XYKey {
+ private final long x;
+ private final long y;
+
+ private XYKey(double x, double y) {
+ this.x = normalizedBits(x);
+ this.y = normalizedBits(y);
+ }
+
+ private static long normalizedBits(double value) {
+ // Primitive double == treats -0.0 and +0.0 as equal, but their raw bit patterns differ.
+ // Canonicalize both so the hash key has the same semantics as an XY coordinate comparison.
+ return Double.doubleToLongBits(value == 0.0 ? 0.0 : value);
+ }
+
+ @Override
+ public boolean equals(Object other) {
+ if (!(other instanceof XYKey)) {
+ return false;
+ }
+ XYKey key = (XYKey) other;
+ return x == key.x && y == key.y;
+ }
+
+ @Override
+ public int hashCode() {
+ int result = Long.hashCode(x);
+ return 31 * result + Long.hashCode(y);
+ }
+ }
+
+ /**
+ * Linear-reference position in a source geometry.
+ *
+ * The segment ordinal is compared first, then the displacement between two samples is
+ * projected onto that segment:
+ *
+ *
+ * segment 4 @ 90 < segment 5 @ 10
+ * segment 5 moving with the source < segment 5 farther along the source
+ *
+ *
+ * This ordering lets two interior samples reveal whether a source traverses a shared fragment
+ * forward or backward without rescanning the original geometry.
+ */
+ private static final class SourceLocation implements Comparable {
+ private final SourceSegment segment;
+ private final Coordinate coordinate;
+
+ private SourceLocation(SourceSegment segment, Coordinate coordinate) {
+ this.segment = segment;
+ this.coordinate = coordinate;
+ }
+
+ @Override
+ public int compareTo(SourceLocation other) {
+ // Segment order dominates position because the flattened ordinal follows source traversal.
+ int segmentComparison = Integer.compare(segment.ordinal, other.segment.ordinal);
+ if (segmentComparison != 0) {
+ return segmentComparison;
+ }
+
+ // Compare the two samples directly rather than normalizing each against a possibly distant
+ // segment endpoint. Scaling the source vector avoids overflow without changing the sign of
+ // the projection. Using both axes also handles fixed-precision overlays which change which
+ // axis is dominant after snapping.
+ double sourceDx = segment.end.getX() - segment.start.getX();
+ double sourceDy = segment.end.getY() - segment.start.getY();
+ double scale = Math.max(Math.abs(sourceDx), Math.abs(sourceDy));
+ double sampleDx = coordinate.getX() - other.coordinate.getX();
+ double sampleDy = coordinate.getY() - other.coordinate.getY();
+ double projection = Math.fma(sampleDx, sourceDx / scale, sampleDy * (sourceDy / scale));
+ return projection == 0.0 ? 0 : Double.compare(projection, 0.0);
+ }
+ }
+}
diff --git a/common/src/test/java/org/apache/sedona/common/GeometryUtilTest.java b/common/src/test/java/org/apache/sedona/common/GeometryUtilTest.java
index 6e67df75cf5..7a9847a8ca2 100644
--- a/common/src/test/java/org/apache/sedona/common/GeometryUtilTest.java
+++ b/common/src/test/java/org/apache/sedona/common/GeometryUtilTest.java
@@ -21,8 +21,10 @@
import java.util.List;
import java.util.Objects;
import org.apache.sedona.common.utils.GeomUtils;
+import org.junit.Assert;
import org.junit.Test;
import org.locationtech.jts.geom.*;
+import org.locationtech.jts.io.WKTReader;
public class GeometryUtilTest {
public static final GeometryFactory GEOMETRY_FACTORY = new GeometryFactory();
@@ -54,4 +56,23 @@ public void extractGeometryCollection() {
GEOMETRY_FACTORY.createGeometryCollection(geoms.toArray(new Geometry[geoms.size()]))),
"GEOMETRYCOLLECTION (POLYGON ((0 1, 3 0, 4 3, 0 4, 0 1)), POLYGON ((3 4, 6 3, 5 5, 3 4)), POINT (5 8), POLYGON ((0 1, 3 0, 4 3, 0 4, 0 1)), POLYGON ((3 4, 6 3, 5 5, 3 4)))"));
}
+
+ @Test
+ public void dimensionalNestedEmptyWKTRoundTrips() throws Exception {
+ WKTReader reader = new WKTReader();
+ for (String wkt :
+ new String[] {
+ "GEOMETRYCOLLECTION Z (LINESTRING Z (0 0 1, 1 1 2), MULTILINESTRING Z EMPTY)",
+ "GEOMETRYCOLLECTION M (LINESTRING M (0 0 1, 1 1 2), MULTILINESTRING M EMPTY)",
+ "GEOMETRYCOLLECTION ZM (LINESTRING ZM (0 0 1 2, 1 1 3 4), " + "MULTILINESTRING ZM EMPTY)"
+ }) {
+ Geometry geometry = reader.read(wkt);
+ String serialized = GeomUtils.getWKT(geometry);
+
+ Assert.assertFalse(serialized.contains("ZEMPTY"));
+ Assert.assertFalse(serialized.contains("MEMPTY"));
+ Assert.assertFalse(serialized.contains("ZMEMPTY"));
+ Assert.assertTrue(geometry.equalsExact(reader.read(serialized)));
+ }
+ }
}
diff --git a/common/src/test/java/org/apache/sedona/common/utils/SharedPathsTest.java b/common/src/test/java/org/apache/sedona/common/utils/SharedPathsTest.java
new file mode 100644
index 00000000000..46835a4d26a
--- /dev/null
+++ b/common/src/test/java/org/apache/sedona/common/utils/SharedPathsTest.java
@@ -0,0 +1,528 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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.apache.sedona.common.utils;
+
+import static org.junit.Assert.*;
+
+import org.apache.sedona.common.Functions;
+import org.apache.sedona.common.geometrySerde.GeometrySerializer;
+import org.junit.Test;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateSequence;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryCollection;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LineString;
+import org.locationtech.jts.geom.MultiLineString;
+import org.locationtech.jts.geom.PrecisionModel;
+import org.locationtech.jts.geom.impl.CoordinateArraySequence;
+import org.locationtech.jts.io.ParseException;
+import org.locationtech.jts.io.WKBReader;
+import org.locationtech.jts.io.WKTReader;
+
+public class SharedPathsTest {
+ private final WKTReader reader = new WKTReader();
+
+ @Test
+ public void matchesPostGISRegressionCases() throws ParseException {
+ assertSharedPaths(
+ "LINESTRING (0 0, 10 0)",
+ "LINESTRING (10 0, 0 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((0 0, 10 0)))");
+ assertSharedPaths(
+ "LINESTRING (0 0, 10 0)",
+ "LINESTRING (20 0, 30 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)");
+ assertSharedPaths(
+ "LINESTRING (0 0, 100 0)",
+ "LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((20 0, 30 0)), MULTILINESTRING ((70 0, 80 0)))");
+ assertSharedPaths(
+ "MULTILINESTRING ((1 3, 4 2, 7 2, 7 5), (13 10, 14 7, 11 6, 15 5))",
+ "LINESTRING (2 1, 4 2, 7 2, 8 3, 10 6, 11 6, 14 7, 16 9)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((4 2, 7 2)), MULTILINESTRING ((14 7, 11 6)))");
+ }
+
+ @Test
+ public void orientsPathsLikeFirstInput() throws ParseException {
+ assertSharedPaths(
+ "LINESTRING (15 0, 5 0)",
+ "LINESTRING (0 0, 10 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((10 0, 5 0)))");
+ }
+
+ @Test
+ public void classifiesShortSharedPathOnLongSourceSegment() throws ParseException {
+ assertSharedPaths(
+ "LINESTRING (-1000000000 0, 1000000000 0)",
+ "LINESTRING (0 0, 0.0000000001 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((0 0, 0.0000000001 0)), " + "MULTILINESTRING EMPTY)");
+ }
+
+ @Test
+ public void locatesFixedPrecisionOverlayOnOriginalLinework() throws ParseException {
+ GeometryFactory fixedFactory = new GeometryFactory(new PrecisionModel(1.0));
+ LineString left =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(0.24, 0.24), new Coordinate(10.24, 0.24)});
+ LineString right =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(5.24, 0.24), new Coordinate(15.24, 0.24)});
+
+ Geometry result = Functions.sharedPaths(left, right);
+
+ Geometry expected =
+ reader.read("GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)");
+ assertTrue("Expected " + expected + " but found " + result, expected.equalsExact(result));
+ assertEquals(fixedFactory.getPrecisionModel(), result.getPrecisionModel());
+ }
+
+ @Test
+ public void classifiesDirectionWhenFixedPrecisionChangesTheDominantAxis() throws ParseException {
+ GeometryFactory fixedFactory = new GeometryFactory(new PrecisionModel(1.0));
+ LineString left =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(0.49, 0.1), new Coordinate(0.51, 0.4)});
+ LineString sameDirection =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(0.49, 0.1), new Coordinate(0.51, 0.4)});
+ LineString oppositeDirection = (LineString) sameDirection.reverse();
+
+ Geometry sameResult = Functions.sharedPaths(left, sameDirection);
+ Geometry oppositeResult = Functions.sharedPaths(left, oppositeDirection);
+
+ Geometry expectedSame =
+ reader.read("GEOMETRYCOLLECTION (MULTILINESTRING ((0 0, 1 0)), MULTILINESTRING EMPTY)");
+ Geometry expectedOpposite =
+ reader.read("GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((0 0, 1 0)))");
+ assertTrue(
+ "Expected " + expectedSame + " but found " + sameResult,
+ expectedSame.equalsExact(sameResult));
+ assertTrue(
+ "Expected " + expectedOpposite + " but found " + oppositeResult,
+ expectedOpposite.equalsExact(oppositeResult));
+ }
+
+ @Test
+ public void keepsFixedPrecisionPathsThatOverlapAfterRounding() throws ParseException {
+ GeometryFactory fixedFactory = new GeometryFactory(new PrecisionModel(1.0));
+ LineString left =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(0, 0.24), new Coordinate(10, 0.24)});
+ LineString right =
+ fixedFactory.createLineString(
+ new Coordinate[] {new Coordinate(5, 0.26), new Coordinate(15, 0.26)});
+
+ assertTrue(left.getEnvelopeInternal().disjoint(right.getEnvelopeInternal()));
+ Geometry result = Functions.sharedPaths(left, right);
+
+ Geometry expected =
+ reader.read("GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)");
+ assertTrue("Expected " + expected + " but found " + result, expected.equalsExact(result));
+ }
+
+ @Test
+ public void choosesTheClosestSegmentBeforeTheEarliestTraversal() throws ParseException {
+ assertSharedPaths(
+ "MULTILINESTRING ((1000000000000010 1, 1000000000000000 1), "
+ + "(1000000000000000 0, 1000000000000010 0))",
+ "LINESTRING (1000000000000000 0, 1000000000000010 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((1000000000000000 0, "
+ + "1000000000000010 0)), MULTILINESTRING EMPTY)");
+ }
+
+ @Test
+ public void handlesClosedAndNonSimpleLines() throws ParseException {
+ assertSharedPaths(
+ "LINESTRING (0 0, 10 0)",
+ "LINESTRING (0 0, 0 10, 10 10, 10 0, 0 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((0 0, 10 0)))");
+ assertSharedPaths(
+ "LINESTRING (0 0, 2 0, 0 0)",
+ "LINESTRING (0 0, 2 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((0 0, 2 0)), MULTILINESTRING EMPTY)");
+ assertSharedPaths(
+ "LINESTRING (0 0, 2 2, 0 2, 2 0)",
+ "LINESTRING (0 0, 2 2)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((0 0, 1 1), (1 1, 2 2)), MULTILINESTRING EMPTY)");
+ }
+
+ @Test
+ public void matchesPostGISDirectionBucketsForAmbiguousNonSimpleTraversals()
+ throws ParseException {
+ String left = "LINESTRING (2.5 0.5, 1 2, 2.5 0.5, 0 0.5, 2 1.5)";
+ String right = "LINESTRING (1 2, 2.5 0.5, 2.000678017814481 2.5, 2.5 0.5)";
+
+ // Both lines revisit the same path, so a shared coordinate can belong to more than one source
+ // traversal. PostGIS does not promise operand-order symmetry for this ambiguous case: swapping
+ // the operands changes the direction buckets. Compare bucket topology because a repeated path
+ // also has no unique output orientation.
+ assertSharedPathBuckets(
+ left,
+ right,
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((2.5 0.5, 1.6666666666666667 "
+ + "1.3333333333333333)), MULTILINESTRING ((1.6666666666666667 "
+ + "1.3333333333333333, 1 2)))");
+ assertSharedPathBuckets(
+ right,
+ left,
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((1 2, "
+ + "1.6666666666666667 1.3333333333333333), (1.6666666666666667 "
+ + "1.3333333333333333, 2.5 0.5)))");
+ }
+
+ @Test
+ public void ignoresPointIntersections() throws ParseException {
+ assertSharedPaths(
+ "LINESTRING (0 0, 10 0)",
+ "LINESTRING (5 -5, 5 5)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)");
+ assertSharedPaths(
+ "LINESTRING (0 0, 10 0)",
+ "LINESTRING (10 0, 20 0)",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)");
+ }
+
+ @Test
+ public void handlesEmptyLinealInputsAndNulls() throws ParseException {
+ Geometry emptyResult =
+ Functions.sharedPaths(
+ reader.read("LINESTRING EMPTY"), reader.read("MULTILINESTRING EMPTY"));
+ assertBucketsEmpty(emptyResult);
+
+ emptyResult =
+ Functions.sharedPaths(
+ reader.read("LINESTRING EMPTY"), reader.read("LINESTRING (0 0, 1 0)"));
+ assertBucketsEmpty(emptyResult);
+
+ Geometry disjointLeft = reader.read("LINESTRING Z (0 0 1, 1 0 2)");
+ Geometry disjointRight = reader.read("LINESTRING Z (2 0 3, 3 0 4)");
+ disjointLeft.setSRID(4326);
+ disjointRight.setSRID(4326);
+ emptyResult = Functions.sharedPaths(disjointLeft, disjointRight);
+ assertEquals(
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)",
+ Functions.asWKT(emptyResult));
+ assertBucketsEmpty(emptyResult);
+ assertEquals(4326, emptyResult.getSRID());
+ assertEquals(4326, emptyResult.getGeometryN(0).getSRID());
+ assertEquals(4326, emptyResult.getGeometryN(1).getSRID());
+ assertNull(Functions.sharedPaths(null, reader.read("LINESTRING (0 0, 1 0)")));
+ assertNull(Functions.sharedPaths(reader.read("LINESTRING (0 0, 1 0)"), null));
+ }
+
+ @Test
+ public void rejectsNonLinealInputs() throws ParseException {
+ Geometry line = reader.read("LINESTRING (0 0, 1 0)");
+ for (String invalidWkt :
+ new String[] {
+ "POINT (100 100)",
+ "POLYGON ((0 0, 1 0, 1 1, 0 0))",
+ "GEOMETRYCOLLECTION (LINESTRING (0 0, 1 0))",
+ "GEOMETRYCOLLECTION EMPTY"
+ }) {
+ IllegalArgumentException error =
+ assertThrows(
+ IllegalArgumentException.class,
+ () -> Functions.sharedPaths(readUnchecked(invalidWkt), line));
+ assertEquals("Geometry is not lineal", error.getMessage());
+
+ error =
+ assertThrows(
+ IllegalArgumentException.class,
+ () -> Functions.sharedPaths(line, readUnchecked(invalidWkt)));
+ assertEquals("Geometry is not lineal", error.getMessage());
+ }
+ }
+
+ @Test
+ public void rejectsMixedSRIDsAndRetainsMatchingSRID() throws ParseException {
+ Geometry left = reader.read("LINESTRING (0 0, 10 0)");
+ Geometry right = reader.read("LINESTRING (20 0, 30 0)");
+ left.setSRID(10);
+ right.setSRID(5);
+
+ IllegalArgumentException error =
+ assertThrows(IllegalArgumentException.class, () -> Functions.sharedPaths(left, right));
+ assertEquals("Operation on mixed SRID geometries (10 != 5)", error.getMessage());
+
+ right.setSRID(10);
+ Geometry result = Functions.sharedPaths(left, right);
+ assertEquals(10, result.getSRID());
+ assertEquals(10, result.getGeometryN(0).getSRID());
+ assertEquals(10, result.getGeometryN(1).getSRID());
+ }
+
+ @Test
+ public void matchesPostGISZAndMBehavior() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING ZM (0 1 5 4, 0 0 6 5, 1 0 7 6, 1 1 8 7)"),
+ reader.read("LINESTRING ZM (0 -1 3 8, 0 0 2 9, 1 0 1 10, 1 -1 0 11)"));
+ CoordinateSequence sequence = firstPath(result).getCoordinateSequence();
+ assertTrue(sequence.hasZ());
+ assertFalse(sequence.hasM());
+ assertArrayEquals(new double[] {6, 7}, zValues(sequence), 0.0);
+
+ result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING (0 1, 0 0, 1 0, 1 1)"),
+ reader.read("LINESTRING Z (0 -1 3, 0 0 2, 1 0 1, 1 -1 0)"));
+ sequence = firstPath(result).getCoordinateSequence();
+ assertTrue(sequence.hasZ());
+ assertArrayEquals(new double[] {2, 1}, zValues(sequence), 0.0);
+
+ result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING M (0 0 1, 10 0 2)"),
+ reader.read("LINESTRING M (0 0 3, 10 0 4)"));
+ sequence = firstPath(result).getCoordinateSequence();
+ assertFalse(sequence.hasZ());
+ assertFalse(sequence.hasM());
+ assertEquals(2, sequence.getDimension());
+ }
+
+ @Test
+ public void writesParseableWKTFor3DResultWithEmptyBucket() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 0 6, 1 0 7)"), reader.read("LINESTRING Z (0 0 1, 1 0 2)"));
+
+ String wkt = Functions.asWKT(result);
+
+ assertEquals(
+ "GEOMETRYCOLLECTION Z(MULTILINESTRING Z((0 0 6, 1 0 7)), " + "MULTILINESTRING Z EMPTY)",
+ wkt);
+ assertTrue(result.equalsExact(reader.read(wkt)));
+ }
+
+ @Test
+ public void downgradesWholeResultWhenAnySharedPathHasNoSourceZ() throws Exception {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("MULTILINESTRING ((0 0 0, 10 0 10), (0 1, 10 1))"),
+ reader.read("MULTILINESTRING ((2 0, 8 0), (2 1, 8 1))"));
+ String expectedWkt =
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((2 0, 8 0), (2 1, 8 1)), " + "MULTILINESTRING EMPTY)";
+ Geometry expected = reader.read(expectedWkt);
+
+ assertEquals(expectedWkt, Functions.asWKT(result));
+ assertFalse(firstPath(result).getCoordinateSequence().hasZ());
+ assertFalse(
+ ((LineString) result.getGeometryN(0).getGeometryN(1)).getCoordinateSequence().hasZ());
+
+ Geometry wktRoundTrip = reader.read(Functions.asWKT(result));
+ assertTrue(expected.equalsExact(wktRoundTrip));
+ Geometry wkbRoundTrip = new WKBReader().read(Functions.asWKB(result));
+ assertTrue(expected.equalsExact(wkbRoundTrip));
+ assertFalse(firstPath(wkbRoundTrip).getCoordinateSequence().hasZ());
+ Geometry serdeRoundTrip = GeometrySerializer.deserialize(GeometrySerializer.serialize(result));
+ assertTrue(expected.equalsExact(serdeRoundTrip));
+ assertFalse(firstPath(serdeRoundTrip).getCoordinateSequence().hasZ());
+ }
+
+ @Test
+ public void selectsAndInterpolatesZFromTheSourceLinework() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (2 0 12, 8 0 18)"),
+ reader.read("LINESTRING Z (0 0 100, 10 0 110)"));
+ assertArrayEquals(
+ new double[] {12, 18}, zValues(firstPath(result).getCoordinateSequence()), 0.0);
+
+ result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 0 0, 10 0 10)"),
+ reader.read("LINESTRING Z (2 0 102, 8 0 108)"));
+ assertArrayEquals(
+ new double[] {102, 108}, zValues(firstPath(result).getCoordinateSequence()), 0.0);
+
+ result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 0 0, 10 0 10)"), reader.read("LINESTRING (2 0, 8 0)"));
+ assertArrayEquals(new double[] {2, 8}, zValues(firstPath(result).getCoordinateSequence()), 0.0);
+ }
+
+ @Test
+ public void selectsZFromTheTraversalContainingTheSharedPath() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 0 0, 2 2 2, 0 2 20, 2 0 22)"),
+ reader.read("LINESTRING (1 1, 2 0)"));
+
+ assertArrayEquals(
+ new double[] {21, 22}, zValues(firstPath(result).getCoordinateSequence()), 0.0);
+ }
+
+ @Test
+ public void disambiguatesAPathWhoseMidpointIsASelfIntersection() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 0 0, 2 2 2, 0 2 20, 2 0 22)"),
+ reader.read("LINESTRING (0 2, 2 0)"));
+
+ MultiLineString sameDirection = (MultiLineString) result.getGeometryN(0);
+ assertEquals(2, sameDirection.getNumGeometries());
+ assertArrayEquals(
+ new double[] {20, 21},
+ zValues(((LineString) sameDirection.getGeometryN(0)).getCoordinateSequence()),
+ 0.0);
+ assertArrayEquals(
+ new double[] {21, 22},
+ zValues(((LineString) sameDirection.getGeometryN(1)).getCoordinateSequence()),
+ 0.0);
+ }
+
+ @Test
+ public void givesAnExactVertexPriorityOverSegmentInterpolation() throws ParseException {
+ Geometry result =
+ Functions.sharedPaths(
+ reader.read("LINESTRING Z (0 1 100, 1 1 999, 2 1 102, 0 0 0, 2 2 2)"),
+ reader.read("LINESTRING (0 0, 2 2)"));
+
+ MultiLineString sameDirection = (MultiLineString) result.getGeometryN(0);
+ assertEquals(2, sameDirection.getNumGeometries());
+ assertArrayEquals(
+ new double[] {0, 999},
+ zValues(((LineString) sameDirection.getGeometryN(0)).getCoordinateSequence()),
+ 0.0);
+ assertArrayEquals(
+ new double[] {999, 2},
+ zValues(((LineString) sameDirection.getGeometryN(1)).getCoordinateSequence()),
+ 0.0);
+ }
+
+ @Test
+ public void indexesZSourcesInsteadOfRescanningThemForEveryOutputCoordinate() {
+ int pointCount = 1000;
+ Coordinate[] coordinates = new Coordinate[pointCount];
+ for (int i = 0; i < pointCount; i++) {
+ coordinates[i] = new Coordinate(i, i % 2, i);
+ }
+ CountingCoordinateSequence leftSequence = new CountingCoordinateSequence(coordinates);
+ CountingCoordinateSequence rightSequence = new CountingCoordinateSequence(coordinates);
+ GeometryFactory factory = new GeometryFactory();
+ LineString left = factory.createLineString(leftSequence);
+ LineString right = factory.createLineString(rightSequence);
+
+ Geometry result = Functions.sharedPaths(left, right);
+
+ MultiLineString sameDirection = (MultiLineString) result.getGeometryN(0);
+ assertEquals(pointCount - 1, sameDirection.getNumGeometries());
+ assertArrayEquals(
+ new double[] {0, 1},
+ zValues(((LineString) sameDirection.getGeometryN(0)).getCoordinateSequence()),
+ 0.0);
+ assertArrayEquals(
+ new double[] {pointCount - 2, pointCount - 1},
+ zValues(((LineString) sameDirection.getGeometryN(pointCount - 2)).getCoordinateSequence()),
+ 0.0);
+ long xyReads = leftSequence.xyReads + rightSequence.xyReads;
+ assertTrue(
+ "Source coordinate reads should scale linearly, but found " + xyReads, xyReads < 100000);
+ long coordinateReads = leftSequence.coordinateReads + rightSequence.coordinateReads;
+ assertTrue(
+ "Source segment reads should scale linearly, but found " + coordinateReads,
+ coordinateReads < 100000);
+ }
+
+ private void assertSharedPaths(String leftWkt, String rightWkt, String expectedWkt)
+ throws ParseException {
+ Geometry actual = Functions.sharedPaths(reader.read(leftWkt), reader.read(rightWkt));
+ Geometry expected = reader.read(expectedWkt);
+ assertTrue("Expected " + expected + " but found " + actual, expected.equalsExact(actual));
+ }
+
+ private void assertSharedPathBuckets(String leftWkt, String rightWkt, String expectedWkt)
+ throws ParseException {
+ Geometry actual = Functions.sharedPaths(reader.read(leftWkt), reader.read(rightWkt));
+ Geometry expected = reader.read(expectedWkt);
+ for (int bucket = 0; bucket < 2; bucket++) {
+ Geometry actualBucket = actual.getGeometryN(bucket);
+ Geometry expectedBucket = expected.getGeometryN(bucket);
+ boolean matches =
+ expectedBucket.isEmpty()
+ ? actualBucket.isEmpty()
+ : expectedBucket.equalsTopo(actualBucket);
+ assertTrue("Expected bucket " + expectedBucket + " but found " + actualBucket, matches);
+ }
+ }
+
+ private void assertBucketsEmpty(Geometry result) {
+ assertTrue(result instanceof GeometryCollection);
+ assertEquals(2, result.getNumGeometries());
+ assertTrue(result.getGeometryN(0) instanceof MultiLineString);
+ assertTrue(result.getGeometryN(1) instanceof MultiLineString);
+ assertTrue(result.getGeometryN(0).isEmpty());
+ assertTrue(result.getGeometryN(1).isEmpty());
+ }
+
+ private LineString firstPath(Geometry result) {
+ return (LineString) result.getGeometryN(0).getGeometryN(0);
+ }
+
+ private double[] zValues(CoordinateSequence sequence) {
+ double[] values = new double[sequence.size()];
+ for (int i = 0; i < sequence.size(); i++) {
+ values[i] = sequence.getZ(i);
+ }
+ return values;
+ }
+
+ private Geometry readUnchecked(String wkt) {
+ try {
+ return reader.read(wkt);
+ } catch (ParseException e) {
+ throw new AssertionError(e);
+ }
+ }
+
+ private static final class CountingCoordinateSequence extends CoordinateArraySequence {
+ private long xyReads;
+ private long coordinateReads;
+
+ private CountingCoordinateSequence(Coordinate[] coordinates) {
+ super(coordinates, 3, 0);
+ }
+
+ @Override
+ public double getX(int index) {
+ xyReads++;
+ return super.getX(index);
+ }
+
+ @Override
+ public double getY(int index) {
+ xyReads++;
+ return super.getY(index);
+ }
+
+ @Override
+ public Coordinate getCoordinate(int index) {
+ coordinateReads++;
+ return super.getCoordinate(index);
+ }
+
+ @Override
+ public Coordinate getCoordinateCopy(int index) {
+ coordinateReads++;
+ return super.getCoordinateCopy(index);
+ }
+ }
+}
diff --git a/docs/api/flink/Geometry-Functions.md b/docs/api/flink/Geometry-Functions.md
index 1dee6619226..452fa05f75d 100644
--- a/docs/api/flink/Geometry-Functions.md
+++ b/docs/api/flink/Geometry-Functions.md
@@ -245,6 +245,7 @@ These functions compute results arising from the overlay of two geometries. Thes
| :--- | :--- | :--- | :--- |
| [ST_Difference](Overlay-Functions/ST_Difference.md) | Geometry | Return the difference between geometry A and B (return part of geometry A that does not intersect geometry B) | v1.5.0 |
| [ST_Intersection](Overlay-Functions/ST_Intersection.md) | Geometry | Return the intersection geometry of A and B | v1.5.0 |
+| [ST_SharedPaths](Overlay-Functions/ST_SharedPaths.md) | Geometry | Returns the paths shared by two lineal geometries, grouped by traversal direction. | v2.0.0 |
| [ST_SubDivide](Overlay-Functions/ST_SubDivide.md) | `Array` | Returns list of geometries divided based of given maximum number of vertices. | v1.5.0 |
| [ST_SymDifference](Overlay-Functions/ST_SymDifference.md) | Geometry | Return the symmetrical difference between geometry A and B (return parts of geometries which are in either of the sets, but not in their intersection) | v1.5.0 |
| [ST_UnaryUnion](Overlay-Functions/ST_UnaryUnion.md) | Geometry | This variant of [ST_Union](Overlay-Functions/ST_Union.md) operates on a single geometry input. The input geometry can be a simple Geometry type, a MultiGeometry, or a GeometryCollection. The function calculates the ge... | v1.6.1 |
diff --git a/docs/api/flink/Overlay-Functions/ST_SharedPaths.md b/docs/api/flink/Overlay-Functions/ST_SharedPaths.md
new file mode 100644
index 00000000000..cc4382c7647
--- /dev/null
+++ b/docs/api/flink/Overlay-Functions/ST_SharedPaths.md
@@ -0,0 +1,197 @@
+
+
+# ST_SharedPaths
+
+Introduction: Returns the paths shared by two lineal geometries, grouped by traversal direction.
+
+Format: `ST_SharedPaths(A: Geometry, B: Geometry)`
+
+Return type: `Geometry`
+
+Since: `v2.0.0`
+
+Both inputs must be a `LineString` or `MultiLineString` and must have the same SRID. The result is
+a `GeometryCollection` containing exactly two `MultiLineString` elements. The first contains paths
+traversed in the same direction by both inputs; the second contains paths traversed in opposite
+directions. Coordinates in both elements follow the direction of `A`.
+
+For non-simple inputs that traverse the same path more than once, direction follows the first
+matching source traversal. Consequently, swapping `A` and `B` can change which result element
+contains a path.
+
+The matching input SRID is retained. A non-empty result retains Z only when at least one input has
+Z and every coordinate of every returned shared path resolves to a finite Z from the input
+geometries. If any returned coordinate does not resolve to a finite source Z, the entire result is
+two-dimensional. M values are always dropped. A result with no shared paths is also
+two-dimensional.
+
+When the inputs do not share a path, including when they intersect only at points, both elements
+are empty `MultiLineString` geometries.
+
+## Visual examples
+
+### Same and opposite directions
+
+Input `B` first follows `A`, leaves it, and later traverses another part of `A` backwards. The
+single result therefore populates both direction buckets. Notice that the path in element 1 is
+reoriented to follow `A`, even though `B` traverses it in the opposite direction.
+
+
+
+### Multipart input and noded output
+
+Input `A` can be a `MultiLineString`. In this example, `B` shares part of two different components
+of `A`. The vertex at `(90 161)` belongs to `B` and splits the shared diagonal into two paths in
+element 0.
+
+
+
+### Point-only intersections
+
+Crossing at an interior point and touching at an endpoint are both zero-dimensional contacts.
+Neither is a shared path, so both `MultiLineString` elements are empty.
+
+
+
+## SQL examples
+
+Wrap the result in `ST_AsText` to display its two direction buckets as WKT.
+
+### Same direction
+
+The shared interval is returned in element 0 because both inputs traverse it from left to right.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (5 0, 15 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)
+```
+
+### Opposite direction
+
+The shared interval is returned in element 1 because `B` traverses it from right to left. Its
+coordinates still follow the left-to-right direction of `A`.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (15 0, 5 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))
+```
+
+### Same and opposite paths in one result
+
+One input pair can populate both elements. Here, `B` first follows `A`, leaves it, and later
+returns along `A` in the opposite direction.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 100 0)'),
+ ST_GeomFromWKT('LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((20 0, 30 0)), MULTILINESTRING ((70 0, 80 0)))
+```
+
+### Multipart input and noded output
+
+Shared paths can come from different components of a `MultiLineString`. The overlay is noded at
+vertices from either input, so a shared path may be shorter than its containing input segment.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT(
+ 'MULTILINESTRING ((1 3, 4 2, 7 2, 7 5), (13 10, 14 7, 11 6, 15 5))'
+ ),
+ ST_GeomFromWKT(
+ 'LINESTRING (2 1, 4 2, 7 2, 8 3, 10 6, 11 6, 14 7, 16 9)'
+ )
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((4 2, 7 2)), MULTILINESTRING ((14 7, 11 6)))
+```
+
+### No lineal overlap
+
+Lines that cross only at a point do not share a path, so both elements are empty.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (5 -5, 5 5)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)
+```
+
+## Practical use: comparing linear networks
+
+`ST_SharedPaths` is useful when two road, rail, pipeline, or utility-network datasets contain
+coincident centerlines and direction matters. For example, it can distinguish a road segment
+digitized consistently in two datasets from a one-way segment digitized in reverse. It can also
+measure how much of an inspected utility route follows its reference alignment in either
+direction.
+
+The following query assumes that a spatial join has already paired candidate segments. Element 0
+measures consistently directed overlap, while element 1 measures reversed overlap:
+
+```sql
+WITH compared AS (
+ SELECT
+ segment_id,
+ ST_SharedPaths(reference_geom, candidate_geom) AS shared
+ FROM network_segment_pairs
+)
+SELECT
+ segment_id,
+ ST_Length(ST_GeometryN(shared, 0)) AS same_direction_length,
+ ST_Length(ST_GeometryN(shared, 1)) AS opposite_direction_length
+FROM compared;
+```
+
+Rows with a positive `opposite_direction_length` are useful candidates for direction or one-way
+attribute review. Apply `ST_SharedPaths` after candidate matching rather than to every possible
+pair in two large networks. The function finds exact lineal overlap; snap or otherwise normalize
+nearly coincident datasets first when their coordinates differ within an accepted tolerance.
diff --git a/docs/api/snowflake/vector-data/Geometry-Functions.md b/docs/api/snowflake/vector-data/Geometry-Functions.md
index cd888855499..de388312887 100644
--- a/docs/api/snowflake/vector-data/Geometry-Functions.md
+++ b/docs/api/snowflake/vector-data/Geometry-Functions.md
@@ -238,6 +238,7 @@ These functions compute results arising from the overlay of two geometries. Thes
| :--- | :--- |
| [ST_Difference](Overlay-Functions/ST_Difference.md) | Return the difference between geometry A and B (return part of geometry A that does not intersect geometry B) |
| [ST_Intersection](Overlay-Functions/ST_Intersection.md) | Return the intersection geometry of A and B |
+| [ST_SharedPaths](Overlay-Functions/ST_SharedPaths.md) | Returns the paths shared by two lineal geometries, grouped by traversal direction. |
| [ST_Split](Overlay-Functions/ST_Split.md) | Split an input geometry by another geometry (called the blade). Linear (LineString or MultiLineString) geometry can be split by a Point, MultiPoint, LineString, MultiLineString, Polygon, or MultiPo... |
| [ST_SubDivide](Overlay-Functions/ST_SubDivide.md) | Returns a multi-geometry divided based of given maximum number of vertices. |
| [ST_SubDivideExplode](Overlay-Functions/ST_SubDivideExplode.md) | It works the same as ST_SubDivide but returns new rows with geometries instead of a multi-geometry. |
diff --git a/docs/api/snowflake/vector-data/Overlay-Functions/ST_SharedPaths.md b/docs/api/snowflake/vector-data/Overlay-Functions/ST_SharedPaths.md
new file mode 100644
index 00000000000..ee42132609b
--- /dev/null
+++ b/docs/api/snowflake/vector-data/Overlay-Functions/ST_SharedPaths.md
@@ -0,0 +1,198 @@
+
+
+# ST_SharedPaths
+
+Introduction: Returns the paths shared by two lineal geometries, grouped by traversal direction.
+
+Format: `ST_SharedPaths(A: geometry, B: geometry)`
+
+Return type: `Geometry`
+
+This function supports Snowflake `GEOMETRY` values. It is not available for `GEOGRAPHY` values.
+Both inputs must be a `LineString` or `MultiLineString`. The result is a `GeometryCollection`
+containing exactly two `MultiLineString` elements. The first contains paths traversed in the same
+direction by both inputs; the second contains paths traversed in opposite directions. Coordinates
+in both elements follow the direction of `A`.
+
+For non-simple inputs that traverse the same path more than once, direction follows the first
+matching source traversal. Consequently, swapping `A` and `B` can change which result element
+contains a path.
+
+Snowflake's native `GEOMETRY` bridge passes values to Sedona as GeoJSON, which does not carry SRID
+metadata. Consequently, this variant cannot reject mixed input SRIDs or preserve a matching SRID;
+its result has SRID 0. The legacy binary geometry interface preserves matching SRIDs. A non-empty
+result retains Z only when at least one input has Z and every coordinate of every returned shared
+path resolves to a finite Z from the input geometries. If any returned coordinate does not resolve
+to a finite source Z, the entire result is two-dimensional. M values are always dropped. A result
+with no shared paths is also two-dimensional.
+
+When the inputs do not share a path, including when they intersect only at points, both elements
+are empty `MultiLineString` geometries.
+
+## Visual examples
+
+### Same and opposite directions
+
+Input `B` first follows `A`, leaves it, and later traverses another part of `A` backwards. The
+single result therefore populates both direction buckets. Notice that the path in element 1 is
+reoriented to follow `A`, even though `B` traverses it in the opposite direction.
+
+
+
+### Multipart input and noded output
+
+Input `A` can be a `MultiLineString`. In this example, `B` shares part of two different components
+of `A`. The vertex at `(90 161)` belongs to `B` and splits the shared diagonal into two paths in
+element 0.
+
+
+
+### Point-only intersections
+
+Crossing at an interior point and touching at an endpoint are both zero-dimensional contacts.
+Neither is a shared path, so both `MultiLineString` elements are empty.
+
+
+
+## SQL examples
+
+Wrap the result in `ST_AsText` to display its two direction buckets as WKT.
+
+### Same direction
+
+The shared interval is returned in element 0 because both inputs traverse it from left to right.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeometryFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeometryFromWKT('LINESTRING (5 0, 15 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)
+```
+
+### Opposite direction
+
+The shared interval is returned in element 1 because `B` traverses it from right to left. Its
+coordinates still follow the left-to-right direction of `A`.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeometryFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeometryFromWKT('LINESTRING (15 0, 5 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))
+```
+
+### Same and opposite paths in one result
+
+One input pair can populate both elements. Here, `B` first follows `A`, leaves it, and later
+returns along `A` in the opposite direction.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeometryFromWKT('LINESTRING (0 0, 100 0)'),
+ ST_GeometryFromWKT('LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((20 0, 30 0)), MULTILINESTRING ((70 0, 80 0)))
+```
+
+### Multipart input and noded output
+
+Shared paths can come from different components of a `MultiLineString`. The overlay is noded at
+vertices from either input, so a shared path may be shorter than its containing input segment.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeometryFromWKT(
+ 'MULTILINESTRING ((1 3, 4 2, 7 2, 7 5), (13 10, 14 7, 11 6, 15 5))'
+ ),
+ ST_GeometryFromWKT(
+ 'LINESTRING (2 1, 4 2, 7 2, 8 3, 10 6, 11 6, 14 7, 16 9)'
+ )
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((4 2, 7 2)), MULTILINESTRING ((14 7, 11 6)))
+```
+
+### No lineal overlap
+
+Lines that cross only at a point do not share a path, so both elements are empty.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeometryFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeometryFromWKT('LINESTRING (5 -5, 5 5)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)
+```
+
+## Practical use: comparing linear networks
+
+`ST_SharedPaths` is useful when two road, rail, pipeline, or utility-network datasets contain
+coincident centerlines and direction matters. For example, it can distinguish a road segment
+digitized consistently in two datasets from a one-way segment digitized in reverse. It can also
+measure how much of an inspected utility route follows its reference alignment in either
+direction.
+
+The following query assumes that a spatial join has already paired candidate segments. Element 0
+measures consistently directed overlap, while element 1 measures reversed overlap:
+
+```sql
+WITH compared AS (
+ SELECT
+ segment_id,
+ ST_SharedPaths(reference_geom, candidate_geom) AS shared
+ FROM network_segment_pairs
+)
+SELECT
+ segment_id,
+ ST_Length(ST_GeometryN(shared, 0)) AS same_direction_length,
+ ST_Length(ST_GeometryN(shared, 1)) AS opposite_direction_length
+FROM compared;
+```
+
+Rows with a positive `opposite_direction_length` are useful candidates for direction or one-way
+attribute review. Apply `ST_SharedPaths` after candidate matching rather than to every possible
+pair in two large networks. The function finds exact lineal overlap; snap or otherwise normalize
+nearly coincident datasets first when their coordinates differ within an accepted tolerance.
diff --git a/docs/api/sql/Geometry-Functions.md b/docs/api/sql/Geometry-Functions.md
index c29b9eefe33..84cd8a8f3f0 100644
--- a/docs/api/sql/Geometry-Functions.md
+++ b/docs/api/sql/Geometry-Functions.md
@@ -248,6 +248,7 @@ These functions compute results arising from the overlay of two geometries. Thes
| :--- | :--- | :--- | :--- |
| [ST_Difference](Overlay-Functions/ST_Difference.md) | Geometry | Return the difference between geometry A and B (return part of geometry A that does not intersect geometry B) | v1.2.0 |
| [ST_Intersection](Overlay-Functions/ST_Intersection.md) | Geometry | Return the intersection geometry of A and B | v1.0.0 |
+| [ST_SharedPaths](Overlay-Functions/ST_SharedPaths.md) | Geometry | Returns the paths shared by two lineal geometries, grouped by traversal direction. | v2.0.0 |
| [ST_Split](Overlay-Functions/ST_Split.md) | Geometry | Split an input geometry by another geometry (called the blade). Linear (LineString or MultiLineString) geometry can be split by a Point, MultiPoint, LineString, MultiLineString, Polygon, or MultiPo... | v1.4.0 |
| [ST_SubDivide](Overlay-Functions/ST_SubDivide.md) | `Array` | Returns list of geometries divided based of given maximum number of vertices. | v1.1.0 |
| [ST_SubDivideExplode](Overlay-Functions/ST_SubDivideExplode.md) | Geometry | It works the same as ST_SubDivide but returns new rows with geometries instead of list. | v1.1.0 |
diff --git a/docs/api/sql/Overlay-Functions/ST_SharedPaths.md b/docs/api/sql/Overlay-Functions/ST_SharedPaths.md
new file mode 100644
index 00000000000..cc4382c7647
--- /dev/null
+++ b/docs/api/sql/Overlay-Functions/ST_SharedPaths.md
@@ -0,0 +1,197 @@
+
+
+# ST_SharedPaths
+
+Introduction: Returns the paths shared by two lineal geometries, grouped by traversal direction.
+
+Format: `ST_SharedPaths(A: Geometry, B: Geometry)`
+
+Return type: `Geometry`
+
+Since: `v2.0.0`
+
+Both inputs must be a `LineString` or `MultiLineString` and must have the same SRID. The result is
+a `GeometryCollection` containing exactly two `MultiLineString` elements. The first contains paths
+traversed in the same direction by both inputs; the second contains paths traversed in opposite
+directions. Coordinates in both elements follow the direction of `A`.
+
+For non-simple inputs that traverse the same path more than once, direction follows the first
+matching source traversal. Consequently, swapping `A` and `B` can change which result element
+contains a path.
+
+The matching input SRID is retained. A non-empty result retains Z only when at least one input has
+Z and every coordinate of every returned shared path resolves to a finite Z from the input
+geometries. If any returned coordinate does not resolve to a finite source Z, the entire result is
+two-dimensional. M values are always dropped. A result with no shared paths is also
+two-dimensional.
+
+When the inputs do not share a path, including when they intersect only at points, both elements
+are empty `MultiLineString` geometries.
+
+## Visual examples
+
+### Same and opposite directions
+
+Input `B` first follows `A`, leaves it, and later traverses another part of `A` backwards. The
+single result therefore populates both direction buckets. Notice that the path in element 1 is
+reoriented to follow `A`, even though `B` traverses it in the opposite direction.
+
+
+
+### Multipart input and noded output
+
+Input `A` can be a `MultiLineString`. In this example, `B` shares part of two different components
+of `A`. The vertex at `(90 161)` belongs to `B` and splits the shared diagonal into two paths in
+element 0.
+
+
+
+### Point-only intersections
+
+Crossing at an interior point and touching at an endpoint are both zero-dimensional contacts.
+Neither is a shared path, so both `MultiLineString` elements are empty.
+
+
+
+## SQL examples
+
+Wrap the result in `ST_AsText` to display its two direction buckets as WKT.
+
+### Same direction
+
+The shared interval is returned in element 0 because both inputs traverse it from left to right.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (5 0, 15 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)
+```
+
+### Opposite direction
+
+The shared interval is returned in element 1 because `B` traverses it from right to left. Its
+coordinates still follow the left-to-right direction of `A`.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (15 0, 5 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))
+```
+
+### Same and opposite paths in one result
+
+One input pair can populate both elements. Here, `B` first follows `A`, leaves it, and later
+returns along `A` in the opposite direction.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 100 0)'),
+ ST_GeomFromWKT('LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((20 0, 30 0)), MULTILINESTRING ((70 0, 80 0)))
+```
+
+### Multipart input and noded output
+
+Shared paths can come from different components of a `MultiLineString`. The overlay is noded at
+vertices from either input, so a shared path may be shorter than its containing input segment.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT(
+ 'MULTILINESTRING ((1 3, 4 2, 7 2, 7 5), (13 10, 14 7, 11 6, 15 5))'
+ ),
+ ST_GeomFromWKT(
+ 'LINESTRING (2 1, 4 2, 7 2, 8 3, 10 6, 11 6, 14 7, 16 9)'
+ )
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING ((4 2, 7 2)), MULTILINESTRING ((14 7, 11 6)))
+```
+
+### No lineal overlap
+
+Lines that cross only at a point do not share a path, so both elements are empty.
+
+```sql
+SELECT ST_AsText(ST_SharedPaths(
+ ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ ST_GeomFromWKT('LINESTRING (5 -5, 5 5)')
+));
+```
+
+Output:
+
+```text
+GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)
+```
+
+## Practical use: comparing linear networks
+
+`ST_SharedPaths` is useful when two road, rail, pipeline, or utility-network datasets contain
+coincident centerlines and direction matters. For example, it can distinguish a road segment
+digitized consistently in two datasets from a one-way segment digitized in reverse. It can also
+measure how much of an inspected utility route follows its reference alignment in either
+direction.
+
+The following query assumes that a spatial join has already paired candidate segments. Element 0
+measures consistently directed overlap, while element 1 measures reversed overlap:
+
+```sql
+WITH compared AS (
+ SELECT
+ segment_id,
+ ST_SharedPaths(reference_geom, candidate_geom) AS shared
+ FROM network_segment_pairs
+)
+SELECT
+ segment_id,
+ ST_Length(ST_GeometryN(shared, 0)) AS same_direction_length,
+ ST_Length(ST_GeometryN(shared, 1)) AS opposite_direction_length
+FROM compared;
+```
+
+Rows with a positive `opposite_direction_length` are useful candidates for direction or one-way
+attribute review. Apply `ST_SharedPaths` after candidate matching rather than to every possible
+pair in two large networks. The function finds exact lineal overlap; snap or otherwise normalize
+nearly coincident datasets first when their coordinates differ within an accepted tolerance.
diff --git a/docs/image/ST_SharedPaths/ST_SharedPaths.svg b/docs/image/ST_SharedPaths/ST_SharedPaths.svg
new file mode 100644
index 00000000000..366eaa48059
--- /dev/null
+++ b/docs/image/ST_SharedPaths/ST_SharedPaths.svg
@@ -0,0 +1,89 @@
+
+ ST_SharedPaths with same-direction and opposite-direction results
+ Four-stage diagram. Input A is shown alone as LINESTRING (0 0, 100 0). Input B is shown alone as LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0). The overlay shows B following A from 20 0 to 30 0 and returning over A from 80 0 to 70 0. The output is GEOMETRYCOLLECTION (MULTILINESTRING ((20 0, 30 0)), MULTILINESTRING ((70 0, 80 0))). Both output paths are oriented like A.
+
+
+
+
+ Same and opposite shared paths
+ See each input separately, then compare the overlay and the two ordered output buckets
+
+ 1 · Input A
+ Reference direction: left → right
+ (0 0)
+ (100 0)
+ A = LINESTRING (0 0, 100 0)
+
+ 2 · Input B
+ Follows A, leaves it, then returns backwards
+ start (20 0)
+ end (70 0)
+ off A
+ B = LINESTRING (20 0, 30 0, 30 50, 80 0, 70 0)
+
+ 3 · Overlay and classify
+ Blue A remains visible between red B's dashes; tinted bands mark the shared portions
+ same direction
+ B goes ←; result follows A →
+ 0
+ 20
+ 30
+ 70
+ 80
+ 100
+ A
+ B
+ same bucket
+ opposite bucket
+
+ 4 · GeometryCollection output
+ Element 0 · same direction
+ MULTILINESTRING ((20 0, 30 0))
+ Element 1 · opposite direction
+ MULTILINESTRING ((70 0, 80 0))
+
+ Exact result
+ GEOMETRYCOLLECTION (
+ MULTILINESTRING ((20 0, 30 0)),
+ MULTILINESTRING ((70 0, 80 0)))
+ Both elements use A's left-to-right coordinate order.
+ Element position—not color—is the SQL result contract.
+
diff --git a/docs/image/ST_SharedPaths/ST_SharedPaths_multipart.svg b/docs/image/ST_SharedPaths/ST_SharedPaths_multipart.svg
new file mode 100644
index 00000000000..be554d459f4
--- /dev/null
+++ b/docs/image/ST_SharedPaths/ST_SharedPaths_multipart.svg
@@ -0,0 +1,86 @@
+
+ ST_SharedPaths with a multipart input and noded output
+ Four-stage diagram. Input A is shown alone as MULTILINESTRING ((26 125, 26 200, 126 200, 126 125, 26 125), (51 150, 101 150, 76 175, 51 150)). Input B is shown alone as LINESTRING (151 100, 126 156.25, 126 125, 90 161, 76 175). The overlay highlights three same-direction shared paths. The output is GEOMETRYCOLLECTION (MULTILINESTRING ((126 156.25, 126 125), (101 150, 90 161), (90 161, 76 175)), MULTILINESTRING EMPTY).
+
+
+
+
+ Multipart input and overlay noding
+ Separate inputs make the components and traversal order visible before they are overlaid
+
+ 1 · Input A · MultiLineString
+ A₁ is the outer loop; A₂ is the inner triangle
+ A₁
+ A₂
+ Both components are shown before B can obscure their edges.
+
+ 2 · Input B · LineString
+ Start at (151 100); finish at (76 175)
+ start (151 100)
+ B vertex (90 161)
+
+ 3 · Overlay and node
+ Green bands identify the shared portions while blue A remains visible between red B's dashes
+ B's (90 161) vertex splits this edge
+ shared loop edge
+ shared triangle edge
+ A
+ B
+ element 0 · same
+
+ 4 · GeometryCollection output
+ Element 0 · same direction
+ (126 156.25, 126 125)
+ (101 150, 90 161)
+ (90 161, 76 175)
+ Element 1 · MULTILINESTRING EMPTY
+
+ Exact inputs and output
+ A = MULTILINESTRING ((26 125, 26 200, 126 200, 126 125, 26 125),
+ (51 150, 101 150, 76 175, 51 150))
+ B = LINESTRING (151 100, 126 156.25, 126 125, 90 161, 76 175)
+ Result = GEOMETRYCOLLECTION (MULTILINESTRING ((126 156.25, 126 125), (101 150, 90 161), (90 161, 76 175)), MULTILINESTRING EMPTY)
+
diff --git a/docs/image/ST_SharedPaths/ST_SharedPaths_point_intersections.svg b/docs/image/ST_SharedPaths/ST_SharedPaths_point_intersections.svg
new file mode 100644
index 00000000000..b539b49db4b
--- /dev/null
+++ b/docs/image/ST_SharedPaths/ST_SharedPaths_point_intersections.svg
@@ -0,0 +1,81 @@
+
+ ST_SharedPaths ignores point-only intersections
+ Two four-stage examples, each showing input A alone, input B alone, their overlay, and the output. In the first, LINESTRING (0 0, 10 0) crosses LINESTRING (5 -5, 5 5) at POINT (5 0). In the second, LINESTRING (0 0, 10 0) touches LINESTRING (10 0, 20 0) at POINT (10 0). Neither case has one-dimensional overlap, so both return GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY).
+
+
+
+
+ Point contact is not a shared path
+ Each input is shown separately before the overlay: only one-dimensional overlap is returned
+
+ Case 1 · Lines cross at an interior point
+ 1 · Input A
+ LINESTRING (0 0, 10 0)
+ left → right
+
+ 2 · Input B
+ LINESTRING (5 -5, 5 5)
+
+ 3 · Overlay
+ Intersection = POINT (5 0)
+
+ 4 · Output
+ A point has no shared length.
+ Element 0 = MULTILINESTRING EMPTY
+ Element 1 = MULTILINESTRING EMPTY
+
+ Case 2 · Lines touch at one endpoint
+ 1 · Input A
+ LINESTRING (0 0, 10 0)
+
+ 2 · Input B
+ LINESTRING (10 0, 20 0)
+
+ 3 · Overlay
+ Intersection = POINT (10 0)
+
+ 4 · Output
+ An endpoint has no shared length.
+ Element 0 = MULTILINESTRING EMPTY
+ Element 1 = MULTILINESTRING EMPTY
+
+ Exact result for both cases
+ GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING EMPTY)
+ ST_Intersection may return the highlighted point, but ST_SharedPaths deliberately does not.
+ Use ST_SharedPaths when shared line length and direction—not simple contact—are the question.
+
diff --git a/flink/src/main/java/org/apache/sedona/flink/Catalog.java b/flink/src/main/java/org/apache/sedona/flink/Catalog.java
index 49391445c23..ac99dc0b16c 100644
--- a/flink/src/main/java/org/apache/sedona/flink/Catalog.java
+++ b/flink/src/main/java/org/apache/sedona/flink/Catalog.java
@@ -223,6 +223,7 @@ public static UserDefinedFunction[] getFuncs() {
new Functions.ST_SimplifyPreserveTopology(),
new Functions.ST_SimplifyVW(),
new Functions.ST_SimplifyPolygonHull(),
+ new Functions.ST_SharedPaths(),
new Functions.ST_Split(),
new Functions.ST_Subdivide(),
new Functions.ST_Segmentize(),
diff --git a/flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java b/flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java
index d45354dd125..58fde5c4eea 100644
--- a/flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java
+++ b/flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java
@@ -3151,6 +3151,28 @@ public Geometry eval(
}
}
+ public static class ST_SharedPaths extends ScalarFunction {
+ @DataTypeHint(
+ value = "RAW",
+ rawSerializer = GeometryTypeSerializer.class,
+ bridgedTo = Geometry.class)
+ public Geometry eval(
+ @DataTypeHint(
+ value = "RAW",
+ rawSerializer = GeometryTypeSerializer.class,
+ bridgedTo = Geometry.class)
+ Object o1,
+ @DataTypeHint(
+ value = "RAW",
+ rawSerializer = GeometryTypeSerializer.class,
+ bridgedTo = Geometry.class)
+ Object o2) {
+ Geometry geom1 = (Geometry) o1;
+ Geometry geom2 = (Geometry) o2;
+ return org.apache.sedona.common.Functions.sharedPaths(geom1, geom2);
+ }
+ }
+
public static class ST_GeometricMedian extends ScalarFunction {
@DataTypeHint(
value = "RAW",
diff --git a/flink/src/test/java/org/apache/sedona/flink/FunctionTest.java b/flink/src/test/java/org/apache/sedona/flink/FunctionTest.java
index 8074120788d..2ffa9c5bd1e 100644
--- a/flink/src/test/java/org/apache/sedona/flink/FunctionTest.java
+++ b/flink/src/test/java/org/apache/sedona/flink/FunctionTest.java
@@ -834,6 +834,18 @@ public void testIntersection() {
assertEquals("POINT (0 0)", result.toString());
}
+ @Test
+ public void testSharedPaths() {
+ Table table =
+ tableEnv.sqlQuery(
+ "SELECT ST_SharedPaths(ST_GeomFromWKT('LINESTRING (0 0, 10 0)'), "
+ + "ST_GeomFromWKT('LINESTRING (15 0, 5 0)'))");
+ Geometry result = (Geometry) first(table).getField(0);
+ assertEquals(
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))",
+ result.toString());
+ }
+
@Test
public void testLength() {
Table polygonTable = createLineStringTable(1);
diff --git a/python/sedona/spark/geopandas/_crs.py b/python/sedona/spark/geopandas/_crs.py
index c7f85509c29..27cb57be7a7 100644
--- a/python/sedona/spark/geopandas/_crs.py
+++ b/python/sedona/spark/geopandas/_crs.py
@@ -19,6 +19,8 @@
from __future__ import annotations
+import warnings
+from types import SimpleNamespace
from typing import Any
from pyproj import CRS
@@ -28,6 +30,36 @@
NO_CRS_OVERRIDE = object()
+def warn_crs_mismatch(left_crs, right_crs, *, stacklevel: int = 3) -> bool:
+ """Emit GeoPandas' standard warning when two geometry CRSs differ.
+
+ ``geopandas.array._crs_mismatch_warn`` adds internal call frames, so its
+ stack level must be higher than the local fallback to attribute both
+ warnings to the same caller.
+ """
+ if left_crs == right_crs:
+ return False
+
+ try:
+ from geopandas.array import _crs_mismatch_warn
+ except ImportError:
+ warnings.warn(
+ f"CRS mismatch between the CRS of left geometries ({left_crs}) "
+ f"and right geometries ({right_crs}).",
+ UserWarning,
+ stacklevel=stacklevel,
+ )
+ else:
+ # GeoPandas' private helper expects objects exposing ``.crs``. Use
+ # immutable values here so distributed objects are never re-inspected.
+ _crs_mismatch_warn(
+ SimpleNamespace(crs=left_crs),
+ SimpleNamespace(crs=right_crs),
+ stacklevel=stacklevel + 1,
+ )
+ return True
+
+
def read_crs_metadata(field: InternalField) -> tuple[bool, CRS | None]:
"""Return whether CRS metadata is present and its decoded value."""
metadata = field.metadata or {}
diff --git a/python/sedona/spark/geopandas/base.py b/python/sedona/spark/geopandas/base.py
index db20cee5a0e..6d718923129 100644
--- a/python/sedona/spark/geopandas/base.py
+++ b/python/sedona/spark/geopandas/base.py
@@ -3622,6 +3622,65 @@ def intersection(self, other, align=None):
"""
return _delegate_to_geometry_column("intersection", self, other, align)
+ def shared_paths(self, other, align=None):
+ """Return the paths shared by each geometry and `other`.
+
+ Input geometries must be ``LineString``, ``LinearRing``, or
+ ``MultiLineString`` geometries. Each result is a
+ ``GeometryCollection`` containing two ``MultiLineString`` geometries.
+ The first contains paths traversed in the same direction by both
+ inputs; the second contains paths traversed in opposite directions.
+
+ The operation works in a 1-to-1 row-wise manner.
+
+ Parameters
+ ----------
+ other : GeoSeries or geometric object
+ The GeoSeries (elementwise) or geometric object to find shared
+ paths with. It must contain only lineal geometries.
+ align : bool | None (default None)
+ If True, automatically align GeoSeries based on their indices.
+ If False, pair values in their existing order. None defaults to
+ True.
+
+ Returns
+ -------
+ GeoSeries
+
+ Examples
+ --------
+ >>> from sedona.spark.geopandas import GeoSeries
+ >>> from shapely.geometry import LineString
+ >>> s = GeoSeries(
+ ... [
+ ... LineString([(0, 0), (1, 0), (2, 0)]),
+ ... LineString([(0, 1), (1, 1), (2, 1)]),
+ ... ]
+ ... )
+ >>> other = LineString([(2, 0), (1, 0), (0, 0)])
+ >>> s.shared_paths(other)
+ 0 GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MUL...
+ 1 GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MUL...
+ dtype: geometry
+
+ Use :meth:`get_geometry` to select the same- or opposite-direction
+ component from each result.
+
+ Notes
+ -----
+ This method follows PostGIS ``ST_SharedPaths`` dimensional semantics:
+ non-empty shared paths retain Z when either input has Z, empty results
+ are two-dimensional, and M coordinates are not retained. The native
+ SQL function rejects mixed SRIDs. At this GeoPandas-compatible layer,
+ differing CRS metadata emits a warning, coordinates are evaluated
+ as-is, and the left CRS is retained.
+
+ See Also
+ --------
+ GeoSeries.get_geometry
+ """
+ return _delegate_to_geometry_column("shared_paths", self, other, align)
+
def shortest_line(self, other, align=None):
"""Returns the shortest line between each geometry in the ``GeoSeries``
and `other`.
diff --git a/python/sedona/spark/geopandas/geoseries.py b/python/sedona/spark/geopandas/geoseries.py
index f0b9d2637d7..52a412f5b7a 100644
--- a/python/sedona/spark/geopandas/geoseries.py
+++ b/python/sedona/spark/geopandas/geoseries.py
@@ -59,6 +59,7 @@
from sedona.spark.geopandas._crs import (
NO_CRS_OVERRIDE,
read_crs_metadata,
+ warn_crs_mismatch,
with_crs_metadata,
)
from sedona.spark.geopandas._explode import expand_geometry_column
@@ -2938,6 +2939,58 @@ def shortest_line(self, other, align=None) -> "GeoSeries":
returns_geom=True,
)
+ def shared_paths(self, other, align=None) -> "GeoSeries":
+ if isinstance(other, BaseGeometry):
+ other_geometry = stc.ST_GeomFromWKB(F.lit(other.wkb))
+ # A standalone Shapely geometry carries no SRID. When CRS metadata
+ # is available, use its scalar SRID so Catalyst can fold this whole
+ # literal expression instead of copying the geometry for every row.
+ has_crs_metadata, metadata_crs = read_crs_metadata(
+ self._internal.data_fields[0]
+ )
+ if has_crs_metadata:
+ scalar_srid = (
+ (metadata_crs.to_epsg() or 0) if metadata_crs is not None else 0
+ )
+ other_geometry = stf.ST_SetSRID(other_geometry, scalar_srid)
+ else:
+ # Spark-created geometry columns may carry an SRID without CRS
+ # metadata, so retain per-row SRID matching for that case.
+ other_geometry = stf.ST_SetSRID(
+ other_geometry, stf.ST_SRID(self.spark.column)
+ )
+ spark_expr = stf.ST_SharedPaths(self.spark.column, other_geometry)
+ return self._query_geometry_column(spark_expr, returns_geom=True)
+
+ if not isinstance(other, (GeoSeries, GeoDataFrame, PandasOnSparkSeries)):
+ raise TypeError(
+ "'other' must be a GeoSeries, GeoDataFrame, "
+ "pandas-on-Spark Series, or geometry"
+ )
+
+ other_series, extended = self._make_series_of_val(other)
+ align = False if extended else align
+ normalize_right_srid = True
+ if isinstance(other_series, GeoSeries):
+ left_crs = self.crs
+ right_crs = other_series.crs
+ normalize_right_srid = warn_crs_mismatch(left_crs, right_crs)
+
+ # GeoPandas warns for a CRS mismatch but still evaluates coordinates
+ # as-is. Assigning the left SRID here preserves that layer contract;
+ # direct ST_SharedPaths calls remain PostGIS-strict about mixed SRIDs.
+ right_geometry = F.col("R")
+ if normalize_right_srid:
+ right_geometry = stf.ST_SetSRID(right_geometry, stf.ST_SRID(F.col("L")))
+ spark_expr = stf.ST_SharedPaths(F.col("L"), right_geometry)
+ return self._row_wise_operation(
+ spark_expr,
+ other_series,
+ align,
+ returns_geom=True,
+ validate_alignment=True,
+ )
+
def snap(self, other, tolerance, align=None) -> "GeoSeries":
if not isinstance(tolerance, (float, int)):
raise NotImplementedError(
@@ -2960,21 +3013,32 @@ def snap(self, other, tolerance, align=None) -> "GeoSeries":
)
return result
- def _boolean_result_preserving_index(
+ def _result_preserving_index(
self,
spark_col: PySparkColumn,
df: pyspark.sql.DataFrame,
index_spark_columns: List[PySparkColumn],
index_fields: List,
index_names: List,
- ) -> pspd.Series:
- """Build a boolean Series while retaining every source index level."""
+ returns_geom: bool = False,
+ keep_name: bool = False,
+ ) -> Union["GeoSeries", pspd.Series]:
+ """Build a Series while retaining every result index level."""
result_index_names = [
f"__index_level_{level}__" for level in range(len(index_spark_columns))
]
+ # Series names are arbitrary hashable Python objects, while Spark
+ # aliases must be strings. Keep the physical name private and apply
+ # the logical Series name only after rebuilding the result.
result_name = SPARK_DEFAULT_SERIES_NAME
+ result_field = None
+ if returns_geom:
+ result_field = self._internal.data_fields[0].copy(name=result_name)
+ result_column = spark_col.alias(result_name, metadata=result_field.metadata)
+ else:
+ result_column = spark_col.alias(result_name)
sdf = df.select(
- spark_col.alias(result_name),
+ result_column,
*[
index_col.alias(result_index_name)
for index_col, result_index_name in zip(
@@ -2994,7 +3058,15 @@ def _boolean_result_preserving_index(
)
for source_field, result_index_name in zip(index_fields, result_index_names)
]
- result_data_field = InternalField.from_struct_field(schema_fields[result_name])
+ if result_field is None:
+ result_field = InternalField.from_struct_field(schema_fields[result_name])
+ else:
+ result_schema_field = schema_fields[result_name]
+ result_field = result_field.copy(
+ spark_type=result_schema_field.dataType,
+ nullable=result_schema_field.nullable,
+ metadata=result_schema_field.metadata,
+ )
internal = InternalFrame(
spark_frame=sdf,
index_spark_columns=[
@@ -3005,15 +3077,37 @@ def _boolean_result_preserving_index(
index_fields=result_index_fields,
column_labels=[(result_name,)],
data_spark_columns=[scol_for(sdf, result_name)],
- data_fields=[result_data_field],
+ data_fields=[result_field],
column_label_names=[None],
)
- return first_series(PandasOnSparkDataFrame(internal)).rename(None)
+ result = first_series(PandasOnSparkDataFrame(internal)).rename(
+ self.name if keep_name else None
+ )
+ return GeoSeries(result) if returns_geom else result
+
+ def _boolean_result_preserving_index(
+ self,
+ spark_col: PySparkColumn,
+ df: pyspark.sql.DataFrame,
+ index_spark_columns: List[PySparkColumn],
+ index_fields: List,
+ index_names: List,
+ ) -> pspd.Series:
+ """Build a boolean Series while retaining every result index level."""
+ return self._result_preserving_index(
+ spark_col,
+ df,
+ index_spark_columns,
+ index_fields,
+ index_names,
+ )
def _align_binary_geometry_series(
self,
other: pspd.Series,
align: Union[bool, None],
+ validate_alignment: bool = True,
+ warning_stacklevel: int = 3,
):
"""Align two geometry Series and preserve every resulting index level."""
position_col = "__binary_geometry_position__"
@@ -3057,68 +3151,72 @@ def _align_binary_geometry_series(
),
F.lit(True).alias(right_present_col),
)
- left_frame = left_frame.orderBy(left_order_col)
- right_frame = right_frame.orderBy(right_order_col)
left_frame = InternalFrame.attach_distributed_sequence_column(
- left_frame, position_col
+ left_frame.orderBy(left_order_col), position_col
)
right_frame = InternalFrame.attach_distributed_sequence_column(
- right_frame, position_col
+ right_frame.orderBy(right_order_col), position_col
)
positional_join = left_frame.join(right_frame, on=position_col, how="outer")
- missing_side = (
- F.col(left_present_col).isNull() | F.col(right_present_col).isNull()
- )
- same_index_structure = len(left_index_aliases) == len(right_index_aliases)
- if same_index_structure:
- index_mismatch = missing_side
- for left_index, right_index in zip(left_index_aliases, right_index_aliases):
- index_mismatch = index_mismatch | ~F.col(left_index).eqNullSafe(
- F.col(right_index)
- )
- else:
- index_mismatch = F.lit(True)
+ indices_match = False
+ if validate_alignment:
+ missing_side = (
+ F.col(left_present_col).isNull() | F.col(right_present_col).isNull()
+ )
+ same_index_structure = len(left_index_aliases) == len(right_index_aliases)
+ if same_index_structure:
+ index_mismatch = missing_side
+ for left_index, right_index in zip(
+ left_index_aliases, right_index_aliases
+ ):
+ index_mismatch = index_mismatch | ~F.col(left_index).eqNullSafe(
+ F.col(right_index)
+ )
+ else:
+ index_mismatch = F.lit(True)
- status = (
- positional_join.select(
- F.col(left_present_col),
- F.col(right_present_col),
- index_mismatch.alias("__index_mismatch__"),
+ status = (
+ positional_join.select(
+ F.col(left_present_col),
+ F.col(right_present_col),
+ index_mismatch.alias("__index_mismatch__"),
+ )
+ .agg(
+ F.count(F.col(left_present_col)).alias("left_count"),
+ F.count(F.col(right_present_col)).alias("right_count"),
+ F.max(F.col("__index_mismatch__").cast("int")).alias(
+ "index_mismatch"
+ ),
+ )
+ .first()
)
- .agg(
- F.count(F.col(left_present_col)).alias("left_count"),
- F.count(F.col(right_present_col)).alias("right_count"),
- F.max(F.col("__index_mismatch__").cast("int")).alias("index_mismatch"),
+ left_count = status["left_count"]
+ right_count = status["right_count"]
+ lengths_match = left_count == right_count
+ indices_match = (
+ same_index_structure
+ and lengths_match
+ and not bool(status["index_mismatch"] or False)
)
- .first()
- )
- left_count = status["left_count"]
- right_count = status["right_count"]
- lengths_match = left_count == right_count
- indices_match = (
- same_index_structure
- and lengths_match
- and not bool(status["index_mismatch"] or False)
- )
- if align is False and not lengths_match:
- raise ValueError(
- "Lengths of inputs do not match. "
- f"Left: {left_count}, Right: {right_count}"
- )
+ if align is False and not lengths_match:
+ raise ValueError(
+ "Lengths of inputs do not match. "
+ f"Left: {left_count}, Right: {right_count}"
+ )
- if align is None and not indices_match:
- warnings.warn(
- "The indices of the left and right GeoSeries' are not equal, "
- "and therefore they will be aligned (reordering and/or "
- "introducing missing values) before executing the operation. "
- "If this alignment is the desired behaviour, you can silence "
- "this warning by passing 'align=True'. If you don't want "
- "alignment and protect yourself of accidentally aligning, "
- "you can pass 'align=False'.",
- stacklevel=3,
- )
+ if align is None and not indices_match:
+ warnings.warn(
+ "The indices of the left and right GeoSeries' are not equal, "
+ "and therefore they will be aligned (reordering and/or "
+ "introducing missing values) before executing the operation. "
+ "If this alignment is the desired behaviour, you can silence "
+ "this warning by passing 'align=True'. If you don't want "
+ "alignment and protect yourself of accidentally aligning, "
+ "you can pass 'align=False'.",
+ stacklevel=warning_stacklevel,
+ )
if align is False or indices_match:
result_index_columns = [
@@ -3355,7 +3453,11 @@ def _geom_equals_exact_series(
result_index_columns,
result_index_fields,
result_index_names,
- ) = self._align_binary_geometry_series(other, align)
+ ) = self._align_binary_geometry_series(
+ other,
+ align,
+ warning_stacklevel=4,
+ )
spark_expr = stp.ST_EqualsExact(F.col("L"), F.col("R"), tolerance)
result = self._boolean_result_preserving_index(
@@ -3367,6 +3469,70 @@ def _geom_equals_exact_series(
)
return _to_bool(result)
+ def _align_single_index_series_lazily(
+ self,
+ other: pspd.Series,
+ align: Union[bool, None],
+ ):
+ """Retain the established one-join plan for single-level indexes."""
+ result_index_column = "__index_level_0__"
+ left_source = self._internal.spark_frame
+ right_source = other._internal.spark_frame
+
+ if align is False:
+ # Natural-order IDs are partition-dependent, so replace them with
+ # distributed sequence IDs before positional pairing.
+ left_frame = InternalFrame.attach_distributed_sequence_column(
+ left_source.drop(NATURAL_ORDER_COLUMN_NAME),
+ NATURAL_ORDER_COLUMN_NAME,
+ ).select(
+ self.spark.column.alias("L"),
+ self._internal.index_spark_columns[0].alias(result_index_column),
+ F.col(NATURAL_ORDER_COLUMN_NAME),
+ )
+ right_frame = InternalFrame.attach_distributed_sequence_column(
+ right_source.drop(NATURAL_ORDER_COLUMN_NAME),
+ NATURAL_ORDER_COLUMN_NAME,
+ ).select(
+ other.spark.column.alias("R"),
+ F.col(NATURAL_ORDER_COLUMN_NAME),
+ )
+ aligned_frame = left_frame.join(
+ right_frame,
+ on=NATURAL_ORDER_COLUMN_NAME,
+ how="outer",
+ )
+ result_index_names = self._internal.index_names
+ else:
+ left_frame = left_source.select(
+ self.spark.column.alias("L"),
+ self._internal.index_spark_columns[0].alias(result_index_column),
+ scol_for(left_source, NATURAL_ORDER_COLUMN_NAME),
+ )
+ right_frame = right_source.select(
+ other.spark.column.alias("R"),
+ other._internal.index_spark_columns[0].alias(result_index_column),
+ )
+ aligned_frame = left_frame.join(
+ right_frame,
+ on=result_index_column,
+ how="outer",
+ )
+ result_index_names = [
+ (
+ self._internal.index_names[0]
+ if self._internal.index_names[0] == other._internal.index_names[0]
+ else None
+ )
+ ]
+
+ return (
+ aligned_frame,
+ [result_index_column],
+ self._internal.index_fields,
+ result_index_names,
+ )
+
def _row_wise_operation(
self,
spark_col: PySparkColumn,
@@ -3375,6 +3541,7 @@ def _row_wise_operation(
returns_geom: bool = False,
default_val: Any = None,
keep_name: bool = False,
+ validate_alignment: bool = False,
):
"""
Helper function to perform a row-wise operation on two GeoSeries.
@@ -3391,46 +3558,41 @@ def _row_wise_operation(
default_val : Any (default None)
The value to use if either L or R is null. If None, nulls are not handled.
- """
- from pyspark.sql.functions import col
- # Note: this is specifically False. None is valid since it defaults to True similar to GeoPandas.
- index_col = (
- NATURAL_ORDER_COLUMN_NAME if align is False else SPARK_DEFAULT_INDEX_NAME
- )
-
- # This code assumes there is only one index (SPARK_DEFAULT_INDEX_NAME)
- # and would need to be updated if Sedona later supports multi-index.
-
- sdf = self._internal.spark_frame
- other_sdf = other._internal.spark_frame
-
- if index_col == NATURAL_ORDER_COLUMN_NAME:
- # NATURAL_ORDER_COLUMN_NAME is not deterministic or sequential, so we instead replace it with a
- # new column of row_numbers to perform the alignment. This mimics the following code.
- # sdf.withColumn(index_col, F.row_number().over(Window.orderBy(NATURAL_ORDER_COLUMN_NAME)))
- sdf = sdf.drop(index_col)
- other_sdf = other_sdf.drop(index_col)
- sdf = self._internal.attach_distributed_sequence_column(sdf, index_col)
- other_sdf = other._internal.attach_distributed_sequence_column(
- other_sdf, index_col
+ validate_alignment : bool, default False
+ If True, validate lengths and index equality eagerly and reproduce
+ GeoPandas' default-alignment warning. Established single-index
+ binary APIs retain their lazy execution path.
+ """
+ # Equal duplicate MultiIndexes pair positionally, while unequal ones
+ # form a Cartesian product per key. Distinguishing those cases needs
+ # the bounded distributed index-status aggregation.
+ validate_alignment = validate_alignment or (
+ align is not False
+ and (
+ len(self._internal.index_spark_columns) > 1
+ or len(other._internal.index_spark_columns) > 1
)
-
- sdf = sdf.select(
- self.spark.column.alias("L"),
- # For the left side:
- # - We always select NATURAL_ORDER_COLUMN_NAME, to avoid having to regenerate it in the result
- # - We always select SPARK_DEFAULT_INDEX_NAME, to retain series index info
- col(NATURAL_ORDER_COLUMN_NAME),
- col(SPARK_DEFAULT_INDEX_NAME),
)
- other_sdf = other_sdf.select(
- other.spark.column.alias("R"),
- # for the right side, we only need the column that we are joining on
- col(index_col),
+ single_index_layout = (
+ len(self._internal.index_spark_columns) == 1
+ and len(other._internal.index_spark_columns) == 1
)
-
- joined_df = sdf.join(other_sdf, on=index_col, how="outer")
+ if not validate_alignment and single_index_layout:
+ alignment_result = self._align_single_index_series_lazily(other, align)
+ else:
+ alignment_result = self._align_binary_geometry_series(
+ other,
+ align,
+ validate_alignment=validate_alignment,
+ warning_stacklevel=4,
+ )
+ (
+ aligned_frame,
+ result_index_columns,
+ result_index_fields,
+ result_index_names,
+ ) = alignment_result
if default_val is not None:
# ps.Series.fillna() doesn't always work for the output for some reason
@@ -3440,9 +3602,12 @@ def _row_wise_operation(
default_val,
).otherwise(spark_col)
- return self._query_geometry_column(
+ return self._result_preserving_index(
spark_col,
- joined_df,
+ aligned_frame,
+ [scol_for(aligned_frame, name) for name in result_index_columns],
+ result_index_fields,
+ result_index_names,
returns_geom=returns_geom,
keep_name=keep_name,
)
diff --git a/python/sedona/spark/geopandas/tools/clip.py b/python/sedona/spark/geopandas/tools/clip.py
index bdbfc292934..be46c6c6806 100644
--- a/python/sedona/spark/geopandas/tools/clip.py
+++ b/python/sedona/spark/geopandas/tools/clip.py
@@ -38,6 +38,7 @@
from sedona.spark.sql import st_constructors as stc
from sedona.spark.sql import st_functions as stf
from sedona.spark.sql import st_predicates as stp
+from sedona.spark.geopandas._crs import warn_crs_mismatch
_POINT_TYPES = ("ST_Point", "ST_MultiPoint")
_LINE_TYPES = ("ST_LineString", "ST_MultiLineString")
@@ -253,23 +254,6 @@ def _as_distributed_mask(mask):
return None
-def _warn_crs_mismatch(obj, mask):
- if obj.crs == mask.crs:
- return
-
- try:
- from geopandas.array import _crs_mismatch_warn
- except ImportError:
- warnings.warn(
- f"CRS mismatch between the CRS of left geometries ({obj.crs}) "
- f"and right geometries ({mask.crs}).",
- UserWarning,
- stacklevel=3,
- )
- else:
- _crs_mismatch_warn(obj, mask, stacklevel=3)
-
-
def _mask_expression(source_sdf, mask, reserved):
"""Return a source frame and a mask column expression."""
rectangle = _mask_is_list_like_rectangle(mask)
@@ -466,7 +450,7 @@ def clip(gdf, mask, keep_geom_type: bool = False, sort: bool = False):
distributed_mask = _as_distributed_mask(mask)
if distributed_mask is not None:
- _warn_crs_mismatch(gdf, distributed_mask)
+ warn_crs_mismatch(gdf.crs, distributed_mask.crs)
mask = distributed_mask
internal, source_sdf, geometry_name = _geometry_context(gdf)
diff --git a/python/sedona/spark/geopandas/tools/overlay.py b/python/sedona/spark/geopandas/tools/overlay.py
index b60b9ce84cd..1cb35f7353f 100644
--- a/python/sedona/spark/geopandas/tools/overlay.py
+++ b/python/sedona/spark/geopandas/tools/overlay.py
@@ -33,7 +33,7 @@
from pyspark.sql import functions as F
from pyspark.sql.types import DoubleType
-from sedona.spark.geopandas._crs import copy_crs_metadata
+from sedona.spark.geopandas._crs import copy_crs_metadata, warn_crs_mismatch
from sedona.spark.sql import st_aggregates as sta
from sedona.spark.sql import st_functions as stf
from sedona.spark.sql import st_predicates as stp
@@ -223,24 +223,6 @@ def _repair_frame(frame: _OverlayFrame, invalid_count: int) -> _OverlayFrame:
return frame._replace(sdf=repaired_sdf)
-def _warn_crs_mismatch(left, right):
- if left.crs == right.crs:
- return
- try:
- from geopandas.array import _crs_mismatch_warn
- except ImportError:
- import warnings
-
- warnings.warn(
- f"CRS mismatch between the CRS of left geometries ({left.crs}) "
- f"and right geometries ({right.crs}).",
- UserWarning,
- stacklevel=3,
- )
- else:
- _crs_mismatch_warn(left, right, stacklevel=3)
-
-
def _candidate_pairs(left: _OverlayFrame, right: _OverlayFrame):
left_geometry = scol_for(left.sdf, left.geometry_name)
right_geometry = scol_for(right.sdf, right.geometry_name)
@@ -718,7 +700,7 @@ def overlay(
raise TypeError("'make_valid' must be a boolean")
make_valid = bool(make_valid)
- _warn_crs_mismatch(df1, df2)
+ warn_crs_mismatch(df1.crs, df2.crs)
left = _normalize_frame(df1, "left")
right = _normalize_frame(df2, "right")
summary = _input_summary(left, right)
diff --git a/python/sedona/spark/sql/st_functions.py b/python/sedona/spark/sql/st_functions.py
index 4a7fbbdffa0..65ab7717fc0 100644
--- a/python/sedona/spark/sql/st_functions.py
+++ b/python/sedona/spark/sql/st_functions.py
@@ -2223,6 +2223,26 @@ def ST_SimplifyPolygonHull(
return _call_st_function("ST_SimplifyPolygonHull", args)
+@validate_argument_types
+def ST_SharedPaths(a: ColumnOrName, b: ColumnOrName) -> Column:
+ """Return the paths shared by two lineal Geometry columns.
+
+ The result is a GeometryCollection containing two MultiLineStrings. The first contains paths
+ traversed in the same direction by both inputs; the second contains paths traversed in opposite
+ directions. Path coordinates follow the direction of ``a``. Inputs must have matching SRIDs;
+ the result retains that SRID. As in PostGIS, a non-empty result retains Z when either input has
+ Z, while M values are dropped.
+
+ :param a: LineString or MultiLineString Geometry column.
+ :type a: ColumnOrName
+ :param b: Other LineString or MultiLineString Geometry column.
+ :type b: ColumnOrName
+ :return: Shared paths grouped by direction as a GeometryCollection column.
+ :rtype: Column
+ """
+ return _call_st_function("ST_SharedPaths", (a, b))
+
+
@validate_argument_types
def ST_Split(input: ColumnOrName, blade: ColumnOrName) -> Column:
"""Split input geometry by the blade geometry.
diff --git a/python/tests/geopandas/test_geoseries.py b/python/tests/geopandas/test_geoseries.py
index cdfd2629ff2..74794f3c829 100644
--- a/python/tests/geopandas/test_geoseries.py
+++ b/python/tests/geopandas/test_geoseries.py
@@ -46,6 +46,11 @@
import pytest
from packaging.version import parse as parse_version
+requires_geopandas_shared_paths = pytest.mark.skipif(
+ parse_version(gpd.__version__) < parse_version("1.0.0"),
+ reason=f"Tests require geopandas>=1.0.0, but found v{gpd.__version__}",
+)
+
@pytest.mark.skipif(
parse_version(shapely.__version__) < parse_version("2.0.0"),
@@ -388,6 +393,11 @@ def test_fillna(self):
]
)
self.check_sgpd_equals_gpd(result, expected)
+
+ numeric_name = sgpd.GeoSeries(gpd.GeoSeries([Point(1, 1), None], name=0))
+ numeric_name_result = numeric_name.fillna(Point(0, 0))
+ numeric_name_expected = gpd.GeoSeries([Point(1, 1), Point(0, 0)], name=0)
+ self.check_sgpd_equals_gpd(numeric_name_result, numeric_name_expected)
result = s.fillna(Polygon([(0, 1), (2, 1), (1, 2)]))
expected = gpd.GeoSeries(
[
@@ -4869,6 +4879,233 @@ def test_shortest_line(self):
)
self.check_sgpd_equals_gpd(df_result, expected)
+ @requires_geopandas_shared_paths
+ def test_shared_paths(self):
+ reference = LineString([(0, 0), (2, 0), (2, 1)])
+ geometries = [
+ LineString([(0, 0), (2, 0), (2, 2)]),
+ LineString([(2, 2), (2, 0), (0, 0)]),
+ MultiLineString(
+ [
+ [(0, 0), (2, 0)],
+ [(2, 2), (2, 1)],
+ ]
+ ),
+ LineString(),
+ None,
+ ]
+ index = pd.MultiIndex.from_tuples(
+ [
+ ("same", 1),
+ ("opposite", 2),
+ ("multi", 3),
+ ("empty", 4),
+ ("null", 5),
+ ],
+ names=["kind", "row"],
+ )
+ source = GeoSeries(
+ geometries,
+ index=index,
+ crs="EPSG:3857",
+ name="roads",
+ )
+ expected = gpd.GeoSeries(
+ geometries,
+ index=index,
+ crs="EPSG:3857",
+ name="roads",
+ ).shared_paths(reference)
+
+ result = source.shared_paths(reference)
+
+ self.check_sgpd_equals_gpd(result, expected)
+ assert result.name is None
+ assert result.crs == source.crs == expected.crs
+ actual = result.to_geopandas()
+ for label in index[:-1]:
+ actual_collection = actual.loc[label]
+ expected_collection = expected.loc[label]
+ assert actual_collection.geom_type == "GeometryCollection"
+ assert len(actual_collection.geoms) == 2
+ for component in range(2):
+ assert actual_collection.geoms[component].equals(
+ expected_collection.geoms[component]
+ )
+ assert actual.loc[("null", 5)] is None
+
+ srids = result._internal.spark_frame.select(
+ stf.ST_SRID(result.spark.column).alias("srid")
+ ).collect()
+ assert {row.srid for row in srids if row.srid is not None} == {3857}
+
+ frame_source = GeoSeries(
+ geometries,
+ index=pd.Index(["same", "opposite", "multi", "empty", "null"]),
+ crs="EPSG:3857",
+ )
+ frame_expected = gpd.GeoSeries(
+ geometries,
+ index=pd.Index(["same", "opposite", "multi", "empty", "null"]),
+ crs="EPSG:3857",
+ ).shared_paths(reference)
+ frame_result = frame_source.to_geoframe().shared_paths(reference)
+ assert isinstance(frame_result, GeoSeries)
+ self.check_sgpd_equals_gpd(frame_result, frame_expected)
+ assert frame_result.crs == frame_source.crs
+
+ with pytest.raises(TypeError, match="'other' must be"):
+ source.shared_paths(None)
+
+ def test_shared_paths_avoids_unnecessary_srid_copies(self, monkeypatch):
+ left = GeoSeries(
+ [LineString([(0, 0), (2, 0)])],
+ crs="EPSG:3857",
+ )
+ right = GeoSeries(
+ [LineString([(0, 0), (1, 0)])],
+ crs="EPSG:3857",
+ )
+ reference = LineString([(0, 0), (1, 0)])
+ original_set_srid = stf.ST_SetSRID
+ set_srid_arguments = []
+
+ def recording_set_srid(geometry, srid):
+ set_srid_arguments.append(srid)
+ return original_set_srid(geometry, srid)
+
+ monkeypatch.setattr(stf, "ST_SetSRID", recording_set_srid)
+
+ matching_result = left.shared_paths(right, align=False)
+ assert set_srid_arguments == []
+ assert matching_result.to_geopandas().iloc[0] is not None
+
+ scalar_result = left.shared_paths(reference)
+ assert set_srid_arguments == [3857]
+ assert scalar_result.to_geopandas().iloc[0] is not None
+
+ @requires_geopandas_shared_paths
+ def test_shared_paths_duplicate_multiindex_alignment(self):
+ left_index = pd.MultiIndex.from_tuples(
+ [("a", 1), ("a", 1), ("b", 2)], names=["group", "row"]
+ )
+ right_index = pd.MultiIndex.from_tuples(
+ [("a", 1), ("a", 1), ("c", 3)], names=["group", "row"]
+ )
+ left_geometries = [
+ LineString([(0, 0), (2, 0)]),
+ LineString([(2, 0), (0, 0)]),
+ LineString([(0, 1), (2, 1)]),
+ ]
+ right_geometries = [
+ LineString([(0, 0), (1, 0)]),
+ LineString([(2, 0), (1, 0)]),
+ LineString([(0, 3), (2, 3)]),
+ ]
+ left = GeoSeries(left_geometries, index=left_index, crs="EPSG:4326")
+ right = GeoSeries(right_geometries, index=right_index, crs="EPSG:4326")
+ expected_left = gpd.GeoSeries(
+ left_geometries, index=left_index, crs="EPSG:4326"
+ )
+ expected_right = gpd.GeoSeries(
+ right_geometries, index=right_index, crs="EPSG:4326"
+ )
+
+ result = left.shared_paths(right, align=True)
+ expected = expected_left.shared_paths(expected_right, align=True)
+
+ actual = result.to_geopandas()
+ assert len(actual) == 6
+ pd.testing.assert_index_equal(actual.index, expected.index)
+ assert result.crs == left.crs
+ for actual_collection, expected_collection in zip(actual, expected):
+ if actual_collection is None or expected_collection is None:
+ assert actual_collection is None and expected_collection is None
+ continue
+ for component in range(2):
+ assert actual_collection.geoms[component].equals(
+ expected_collection.geoms[component]
+ )
+
+ positional_result = left.shared_paths(right, align=False)
+ positional_expected = expected_left.shared_paths(expected_right, align=False)
+ positional_actual = positional_result.to_geopandas()
+ pd.testing.assert_index_equal(
+ positional_actual.index, positional_expected.index
+ )
+ for actual_collection, expected_collection in zip(
+ positional_actual, positional_expected
+ ):
+ for component in range(2):
+ assert actual_collection.geoms[component].equals(
+ expected_collection.geoms[component]
+ )
+
+ @requires_geopandas_shared_paths
+ def test_shared_paths_default_alignment_warning_and_length_validation(self):
+ left = GeoSeries(
+ [LineString([(0, 0), (2, 0)])],
+ index=pd.Index(["left"], name="feature"),
+ )
+ right = GeoSeries(
+ [LineString([(0, 0), (1, 0)])],
+ index=pd.Index(["right"], name="feature"),
+ )
+
+ with pytest.warns(
+ UserWarning,
+ match="The indices of the left and right GeoSeries' are not equal",
+ ):
+ result = left.shared_paths(right)
+ expected = gpd.GeoSeries(
+ [LineString([(0, 0), (2, 0)])],
+ index=pd.Index(["left"], name="feature"),
+ ).shared_paths(
+ gpd.GeoSeries(
+ [LineString([(0, 0), (1, 0)])],
+ index=pd.Index(["right"], name="feature"),
+ ),
+ align=True,
+ )
+ self.check_sgpd_equals_gpd(result, expected)
+
+ with pytest.warns(UserWarning, match="CRS mismatch") as warning_info:
+ crs_result = GeoSeries(
+ [LineString([(0, 0), (2, 0)])], crs="EPSG:4326"
+ ).shared_paths(
+ GeoSeries([LineString([(0, 0), (1, 0)])], crs="EPSG:3857"),
+ align=False,
+ )
+ crs_warning = next(
+ warning
+ for warning in warning_info
+ if "CRS mismatch" in str(warning.message)
+ )
+ assert crs_warning.filename == __file__
+ assert crs_result.crs.to_epsg() == 4326
+ with pytest.warns(UserWarning, match="CRS mismatch"):
+ crs_expected = gpd.GeoSeries(
+ [LineString([(0, 0), (2, 0)])], crs="EPSG:4326"
+ ).shared_paths(
+ gpd.GeoSeries([LineString([(0, 0), (1, 0)])], crs="EPSG:3857"),
+ align=False,
+ )
+ self.check_sgpd_equals_gpd(crs_result, crs_expected)
+
+ with pytest.raises(
+ ValueError,
+ match=r"Lengths of inputs do not match\. Left: 1, Right: 2",
+ ):
+ left.shared_paths(
+ GeoSeries(
+ [
+ LineString([(0, 0), (1, 0)]),
+ LineString([(0, 1), (1, 1)]),
+ ]
+ ),
+ align=False,
+ )
+
def test_intersection_all(self):
s = GeoSeries([box(0, 0, 2, 2), box(1, 1, 3, 3)])
result = s.intersection_all()
@@ -5167,9 +5404,92 @@ def test_binary_operation_with_projected_multiindex(self):
result = GeoSeries([Point(0, 0), Point(1, 1)], index=index).geom_equals(
Point(0, 0)
)
- expected = pd.Series([True, False], index=pd.Index(["a", "b"], name="group"))
+ expected = pd.Series([True, False], index=index)
self.check_pd_series_equal(result, expected)
+ def test_binary_operation_single_index_alignment_stays_lazy(self, monkeypatch):
+ left_geometries = [box(0, 0, 2, 2), box(3, 3, 5, 5)]
+ right_geometries = [box(1, 1, 4, 4), box(0, 0, 1, 1)]
+ left = GeoSeries(left_geometries, index=["a", "b"])
+ right = GeoSeries(right_geometries, index=["b", "a"])
+ sequence_attachments = []
+ original_attach = InternalFrame.attach_distributed_sequence_column
+
+ def recording_attach(sdf, column_name):
+ sequence_attachments.append(column_name)
+ return original_attach(sdf, column_name)
+
+ monkeypatch.setattr(
+ InternalFrame,
+ "attach_distributed_sequence_column",
+ recording_attach,
+ )
+
+ result = left.intersection(right, align=True)
+ assert sequence_attachments == []
+ expected = gpd.GeoSeries(left_geometries, index=["a", "b"]).intersection(
+ gpd.GeoSeries(right_geometries, index=["b", "a"]), align=True
+ )
+ self.check_sgpd_equals_gpd(result, expected)
+
+ @pytest.mark.parametrize(
+ ("method", "kwargs"),
+ [
+ ("intersection", {}),
+ ("shortest_line", {}),
+ ("snap", {"tolerance": 0.25}),
+ ],
+ )
+ def test_binary_geometry_operations_preserve_duplicate_multiindex(
+ self, method, kwargs
+ ):
+ left_index = pd.MultiIndex.from_tuples(
+ [("a", 1), ("a", 1), ("b", 2)], names=["group", "row"]
+ )
+ right_index = pd.MultiIndex.from_tuples(
+ [("a", 1), ("a", 1), ("c", 3)], names=["group", "row"]
+ )
+ left = GeoSeries(
+ [
+ LineString([(0, 0), (2, 0)]),
+ LineString([(0, 0), (2, 0)]),
+ LineString([(0, 1), (2, 1)]),
+ ],
+ index=left_index,
+ )
+ right = GeoSeries(
+ [
+ LineString([(0, 0), (1, 0)]),
+ LineString([(1, 0), (2, 0)]),
+ LineString([(0, 3), (2, 3)]),
+ ],
+ index=right_index,
+ )
+
+ result = getattr(left, method)(right, align=True, **kwargs)
+ expected = getattr(gpd.GeoSeries(left.to_geopandas()), method)(
+ gpd.GeoSeries(right.to_geopandas()),
+ align=True,
+ **kwargs,
+ )
+
+ self.check_sgpd_equals_gpd(result, expected)
+ assert len(result) == 6
+
+ identical_right = GeoSeries(
+ right.to_geopandas().array,
+ index=left_index,
+ )
+ identical_result = getattr(left, method)(identical_right, align=True, **kwargs)
+ identical_expected = getattr(gpd.GeoSeries(left.to_geopandas()), method)(
+ gpd.GeoSeries(identical_right.to_geopandas()),
+ align=True,
+ **kwargs,
+ )
+
+ self.check_sgpd_equals_gpd(identical_result, identical_expected)
+ assert len(identical_result) == 3
+
def test_geom_equals_exact(self):
s = GeoSeries([Point(0, 1.1), Point(0, 1.0), Point(0, 1.2)])
diff --git a/python/tests/geopandas/test_match_geopandas_series.py b/python/tests/geopandas/test_match_geopandas_series.py
index d3b95193893..cd96d0170e0 100644
--- a/python/tests/geopandas/test_match_geopandas_series.py
+++ b/python/tests/geopandas/test_match_geopandas_series.py
@@ -1969,6 +1969,46 @@ def test_shortest_line(self):
)
self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+ def test_shared_paths(self):
+ if parse_version(gpd.__version__) < parse_version("1.0.0"):
+ pytest.skip("geopandas < 1.0.0 does not support shared_paths")
+
+ lines = [
+ LineString([(0, 0), (0.5, 0.5), (1, 0)]),
+ LineString([(0, 1), (0.5, 0.5), (1, 1)]),
+ MultiLineString(
+ [
+ [(0, 0), (0.5, 0.5)],
+ [(1, 1), (0.5, 0.5)],
+ ]
+ ),
+ LineString(),
+ None,
+ ]
+ reference = LineString([(0, 0), (0.5, 0.5), (0, 1)])
+
+ sgpd_result = GeoSeries(lines, crs="EPSG:3857").shared_paths(reference)
+ gpd_result = gpd.GeoSeries(lines, crs="EPSG:3857").shared_paths(reference)
+ self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+
+ other_lines = [
+ LineString([(0, 0), (0.5, 0.5), (1, 0)]),
+ LineString([(1, 1), (0.5, 0.5), (0, 1)]),
+ MultiLineString(
+ [
+ [(0.5, 0.5), (0, 0)],
+ [(0.5, 0.5), (1, 1)],
+ ]
+ ),
+ LineString(),
+ None,
+ ]
+ sgpd_result = GeoSeries(lines).shared_paths(GeoSeries(other_lines), align=False)
+ gpd_result = gpd.GeoSeries(lines).shared_paths(
+ gpd.GeoSeries(other_lines), align=False
+ )
+ self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+
def test_intersection_all(self):
pass
diff --git a/python/tests/sql/test_dataframe_api.py b/python/tests/sql/test_dataframe_api.py
index 522534eb834..60a6d6bb5aa 100644
--- a/python/tests/sql/test_dataframe_api.py
+++ b/python/tests/sql/test_dataframe_api.py
@@ -713,6 +713,16 @@
"",
"POLYGON ((2 0, 1 0, 1 1, 2 1, 2 0))",
),
+ (
+ stf.ST_SharedPaths,
+ (
+ lambda: f.expr("ST_GeomFromWKT('LINESTRING (0 0, 10 0)')"),
+ lambda: f.expr("ST_GeomFromWKT('LINESTRING (15 0, 5 0)')"),
+ ),
+ "null",
+ "",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))",
+ ),
(stf.ST_InterpolatePoint, ("linem", "point"), "linestringm_and_point", "", 1.0),
(stf.ST_IsCollection, ("geom",), "geom_collection", "", True),
(stf.ST_IsClosed, ("geom",), "closed_linestring_geom", "", True),
@@ -1524,6 +1534,8 @@
(stf.ST_InteriorRingN, ("", 0.0)),
(stf.ST_Intersection, (None, "")),
(stf.ST_Intersection, ("", None)),
+ (stf.ST_SharedPaths, (None, "")),
+ (stf.ST_SharedPaths, ("", None)),
(stf.ST_IsClosed, (None,)),
(stf.ST_IsEmpty, (None,)),
(stf.ST_IsLineStringCCW, (None,)),
diff --git a/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctions.java b/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctions.java
index 44087387d5f..751a73a6ac7 100644
--- a/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctions.java
+++ b/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctions.java
@@ -529,6 +529,14 @@ public void test_ST_Intersection() {
"POINT (1 1)");
}
+ @Test
+ public void test_ST_SharedPaths() {
+ registerUDF("ST_SharedPaths", byte[].class, byte[].class);
+ verifySqlSingleRes(
+ "select sedona.ST_AsText(sedona.ST_SharedPaths(sedona.ST_GeomFromText('LINESTRING(0 0, 10 0)'), sedona.ST_GeomFromText('LINESTRING(15 0, 5 0)')))",
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))");
+ }
+
@Test
public void test_ST_IsClosed() {
registerUDF("ST_IsClosed", byte[].class);
diff --git a/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctionsV2.java b/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctionsV2.java
index 73c2da9d15f..6757377904e 100644
--- a/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctionsV2.java
+++ b/snowflake-tester/src/test/java/org/apache/sedona/snowflake/snowsql/TestFunctionsV2.java
@@ -509,6 +509,14 @@ public void test_ST_Intersection() {
"POINT(1 1)");
}
+ @Test
+ public void test_ST_SharedPaths() {
+ registerUDFV2("ST_SharedPaths", String.class, String.class);
+ verifySqlSingleRes(
+ "select ST_AsText(sedona.ST_SharedPaths(ST_GeometryFromWKT('LINESTRING(0 0, 10 0)'), ST_GeometryFromWKT('LINESTRING(15 0, 5 0)')))",
+ "GEOMETRYCOLLECTION(MULTILINESTRING EMPTY,MULTILINESTRING((5 0,10 0)))");
+ }
+
@Test
public void test_ST_IsClosed() {
registerUDFV2("ST_IsClosed", String.class);
diff --git a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFs.java b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFs.java
index 46390a8bfa5..f6eae3f91f5 100644
--- a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFs.java
+++ b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFs.java
@@ -558,6 +558,13 @@ public static byte[] ST_Intersection(byte[] leftGeometry, byte[] rightGeometry)
GeometrySerde.deserialize(leftGeometry), GeometrySerde.deserialize(rightGeometry)));
}
+ @UDFAnnotations.ParamMeta(argNames = {"leftGeometry", "rightGeometry"})
+ public static byte[] ST_SharedPaths(byte[] leftGeometry, byte[] rightGeometry) {
+ return GeometrySerde.serialize(
+ Functions.sharedPaths(
+ GeometrySerde.deserialize(leftGeometry), GeometrySerde.deserialize(rightGeometry)));
+ }
+
@UDFAnnotations.ParamMeta(argNames = {"leftGeometry", "rightGeometry"})
public static boolean ST_Intersects(byte[] leftGeometry, byte[] rightGeometry) {
return Predicates.intersects(
diff --git a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFsV2.java b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFsV2.java
index f7977a3b418..5cdcb68dfd5 100644
--- a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFsV2.java
+++ b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/UDFsV2.java
@@ -703,6 +703,17 @@ public static String ST_Intersection(String leftGeometry, String rightGeometry)
GeometrySerde.deserGeoJson(leftGeometry), GeometrySerde.deserGeoJson(rightGeometry)));
}
+ @UDFAnnotations.GeometryOnly
+ @UDFAnnotations.ParamMeta(
+ argNames = {"leftGeometry", "rightGeometry"},
+ argTypes = {"Geometry", "Geometry"},
+ returnTypes = "Geometry")
+ public static String ST_SharedPaths(String leftGeometry, String rightGeometry) {
+ return GeometrySerde.serGeoJson(
+ Functions.sharedPaths(
+ GeometrySerde.deserGeoJson(leftGeometry), GeometrySerde.deserGeoJson(rightGeometry)));
+ }
+
@UDFAnnotations.ParamMeta(
argNames = {"leftGeometry", "rightGeometry"},
argTypes = {"Geometry", "Geometry"})
diff --git a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/annotations/UDFAnnotations.java b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/annotations/UDFAnnotations.java
index b480ad74712..79f694da0e1 100644
--- a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/annotations/UDFAnnotations.java
+++ b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/annotations/UDFAnnotations.java
@@ -36,6 +36,11 @@ public class UDFAnnotations {
@Target(ElementType.METHOD)
public static @interface Volatile {}
+ /** Marks a UDFV2 method that must not be generated for Snowflake GEOGRAPHY arguments. */
+ @Retention(RetentionPolicy.RUNTIME)
+ @Target(ElementType.METHOD)
+ public static @interface GeometryOnly {}
+
@Retention(RetentionPolicy.RUNTIME)
@Target(ElementType.METHOD)
public static @interface ParamMeta {
diff --git a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGenerator.java b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGenerator.java
index ff6d39a243a..b82dffa48db 100644
--- a/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGenerator.java
+++ b/snowflake/src/main/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGenerator.java
@@ -51,6 +51,11 @@ public static String buildUDFDDL(
if (!method.isAnnotationPresent(UDFAnnotations.ParamMeta.class)) {
throw new RuntimeException("Missing ParamMeta annotation for method: " + method.getName());
}
+ if (method.isAnnotationPresent(UDFAnnotations.GeometryOnly.class)
+ && "GEOGRAPHY".equals(Constants.snowflakeTypeMap.get("Geometry"))) {
+ throw new IllegalArgumentException(
+ "Cannot generate GEOGRAPHY DDL for @GeometryOnly method: " + method.getName());
+ }
String[] argNames = method.getAnnotation(UDFAnnotations.ParamMeta.class).argNames();
Parameter[] argTypesRaw = method.getParameters();
String argTypesCustom[] = method.getAnnotation(UDFAnnotations.ParamMeta.class).argTypes();
@@ -101,11 +106,22 @@ public static List buildAll(
ddlList.add(buildUDFDDL(method, configs, stageName, isNativeApp, appRoleName));
}
}
- // Replace Geometry with GEOGRAPHY and generate DDL for UDFsV2 again
- Constants.snowflakeTypeMap.replace("Geometry", "GEOGRAPHY");
- for (Method method : udfV2Methods()) {
- if (method.getModifiers() == (Modifier.PUBLIC | Modifier.STATIC)) {
- ddlList.add(buildUDFDDL(method, configs, stageName, isNativeApp, appRoleName));
+ // Replace Geometry with GEOGRAPHY and generate DDL for UDFsV2 again. Restore the shared type
+ // map so repeated generator calls and direct buildUDFDDL calls still target GEOMETRY by
+ // default.
+ String originalGeometryType = Constants.snowflakeTypeMap.put("Geometry", "GEOGRAPHY");
+ try {
+ for (Method method : udfV2Methods()) {
+ if (method.getModifiers() == (Modifier.PUBLIC | Modifier.STATIC)
+ && !method.isAnnotationPresent(UDFAnnotations.GeometryOnly.class)) {
+ ddlList.add(buildUDFDDL(method, configs, stageName, isNativeApp, appRoleName));
+ }
+ }
+ } finally {
+ if (originalGeometryType == null) {
+ Constants.snowflakeTypeMap.remove("Geometry");
+ } else {
+ Constants.snowflakeTypeMap.put("Geometry", originalGeometryType);
}
}
return ddlList;
diff --git a/snowflake/src/test/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGeneratorTest.java b/snowflake/src/test/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGeneratorTest.java
new file mode 100644
index 00000000000..1980f546cdf
--- /dev/null
+++ b/snowflake/src/test/java/org/apache/sedona/snowflake/snowsql/ddl/UDFDDLGeneratorTest.java
@@ -0,0 +1,108 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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.apache.sedona.snowflake.snowsql.ddl;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertFalse;
+import static org.junit.Assert.assertThrows;
+import static org.junit.Assert.assertTrue;
+
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.stream.Collectors;
+import org.apache.sedona.snowflake.snowsql.UDFsV2;
+import org.junit.Test;
+
+public class UDFDDLGeneratorTest {
+
+ @Test
+ public void geometryOnlyFunctionsAreNotGeneratedForGeography() {
+ String originalGeometryType = Constants.snowflakeTypeMap.put("Geometry", "GEOMETRY");
+
+ try {
+ assertSharedPathsTargets(UDFDDLGenerator.buildAll(configs(), "@ApacheSedona", false, ""));
+ } finally {
+ restoreGeometryType(originalGeometryType);
+ }
+ }
+
+ @Test
+ public void geometryOnlyFunctionsRejectDirectGeographyGeneration() throws NoSuchMethodException {
+ String originalGeometryType = Constants.snowflakeTypeMap.put("Geometry", "GEOGRAPHY");
+
+ try {
+ IllegalArgumentException error =
+ assertThrows(
+ IllegalArgumentException.class,
+ () ->
+ UDFDDLGenerator.buildUDFDDL(
+ UDFsV2.class.getMethod("ST_SharedPaths", String.class, String.class),
+ configs(),
+ "@ApacheSedona",
+ false,
+ ""));
+
+ assertEquals(
+ "Cannot generate GEOGRAPHY DDL for @GeometryOnly method: ST_SharedPaths",
+ error.getMessage());
+ } finally {
+ restoreGeometryType(originalGeometryType);
+ }
+ }
+
+ @Test
+ public void buildAllRestoresGeometryTypeAndCanBeRepeated() {
+ String originalGeometryType = Constants.snowflakeTypeMap.put("Geometry", "GEOMETRY");
+
+ try {
+ assertSharedPathsTargets(UDFDDLGenerator.buildAll(configs(), "@ApacheSedona", false, ""));
+ assertEquals("GEOMETRY", Constants.snowflakeTypeMap.get("Geometry"));
+
+ assertSharedPathsTargets(UDFDDLGenerator.buildAll(configs(), "@ApacheSedona", false, ""));
+ assertEquals("GEOMETRY", Constants.snowflakeTypeMap.get("Geometry"));
+ } finally {
+ restoreGeometryType(originalGeometryType);
+ }
+ }
+
+ private static Map configs() {
+ Map configs = new HashMap<>();
+ configs.put(Constants.SEDONA_VERSION, "test");
+ configs.put(Constants.GEOTOOLS_VERSION, "test");
+ return configs;
+ }
+
+ private static void assertSharedPathsTargets(List ddls) {
+ List sharedPathsDdls =
+ ddls.stream().filter(ddl -> ddl.contains(".ST_SharedPaths ")).collect(Collectors.toList());
+
+ assertEquals(2, sharedPathsDdls.size());
+ assertTrue(sharedPathsDdls.stream().anyMatch(ddl -> ddl.contains(" GEOMETRY")));
+ assertFalse(sharedPathsDdls.stream().anyMatch(ddl -> ddl.contains(" GEOGRAPHY")));
+ }
+
+ private static void restoreGeometryType(String geometryType) {
+ if (geometryType == null) {
+ Constants.snowflakeTypeMap.remove("Geometry");
+ } else {
+ Constants.snowflakeTypeMap.put("Geometry", geometryType);
+ }
+ }
+}
diff --git a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
index 895b03ff22a..47e2e88f991 100644
--- a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
+++ b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
@@ -240,6 +240,7 @@ object Catalog extends AbstractCatalog with Logging {
val overlayExprs: Seq[FunctionDescription] = Seq(
function[ST_Difference](),
function[ST_Intersection](),
+ function[ST_SharedPaths](),
function[ST_Split](),
function[ST_SubDivide](),
function[ST_SubDivideExplode](),
diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
index 305506bdc8e..7c4b325076b 100644
--- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
+++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
@@ -1405,6 +1405,20 @@ private[apache] case class ST_Difference(inputExpressions: Seq[Expression])
}
}
+/**
+ * Return the paths shared by two lineal geometries, separated into same-direction and
+ * opposite-direction paths.
+ *
+ * @param inputExpressions
+ */
+private[apache] case class ST_SharedPaths(inputExpressions: Seq[Expression])
+ extends InferredExpression(Functions.sharedPaths _) {
+
+ protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+ copy(inputExpressions = newChildren)
+ }
+}
+
/**
* Return the symmetrical difference between geometry A and B
*
diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
index b34a6a608ea..cf229fa46f7 100644
--- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
+++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
@@ -754,6 +754,9 @@ object st_functions {
def ST_SimplifyPreserveTopology(geometry: String, distanceTolerance: Double): Column =
wrapExpression[ST_SimplifyPreserveTopology](geometry, distanceTolerance)
+ def ST_SharedPaths(a: Column, b: Column): Column = wrapExpression[ST_SharedPaths](a, b)
+ def ST_SharedPaths(a: String, b: String): Column = wrapExpression[ST_SharedPaths](a, b)
+
def ST_Split(input: Column, blade: Column): Column = wrapExpression[ST_Split](input, blade)
def ST_Split(input: String, blade: String): Column = wrapExpression[ST_Split](input, blade)
diff --git a/spark/common/src/test/scala/org/apache/sedona/sql/dataFrameAPITestScala.scala b/spark/common/src/test/scala/org/apache/sedona/sql/dataFrameAPITestScala.scala
index 40f90338e72..d7fd7df1e9b 100644
--- a/spark/common/src/test/scala/org/apache/sedona/sql/dataFrameAPITestScala.scala
+++ b/spark/common/src/test/scala/org/apache/sedona/sql/dataFrameAPITestScala.scala
@@ -1161,6 +1161,19 @@ class dataFrameAPITestScala extends TestBaseScala {
assert(actualResult == expectedResult)
}
+ it("Passed ST_SharedPaths") {
+ val lineDf = sparkSession.sql(
+ "SELECT ST_GeomFromWKT('LINESTRING (0 0, 10 0)') AS a, ST_GeomFromWKT('LINESTRING (15 0, 5 0)') AS b")
+ val expectedResult =
+ "GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((5 0, 10 0)))"
+
+ val stringResult = lineDf.select(ST_SharedPaths("a", "b")).first().getAs[Geometry](0)
+ assertEquals(expectedResult, stringResult.toText)
+
+ val columnResult = lineDf.select(ST_SharedPaths($"a", $"b")).first().getAs[Geometry](0)
+ assertEquals(expectedResult, columnResult.toText)
+ }
+
it("Passed ST_Union") {
val polygonDf = sparkSession.sql(
"SELECT ST_GeomFromWKT('POLYGON ((-3 -3, 3 -3, 3 3, -3 3, -3 -3))') AS a, ST_GeomFromWKT('POLYGON ((-2 1, 2 1, 2 4, -2 4, -2 1))') AS b")
diff --git a/spark/common/src/test/scala/org/apache/sedona/sql/functionTestScala.scala b/spark/common/src/test/scala/org/apache/sedona/sql/functionTestScala.scala
index 75fb6f1664b..4870eb5e4b0 100644
--- a/spark/common/src/test/scala/org/apache/sedona/sql/functionTestScala.scala
+++ b/spark/common/src/test/scala/org/apache/sedona/sql/functionTestScala.scala
@@ -1603,6 +1603,41 @@ class functionTestScala
expected = "POLYGON ((0.5 10.7, 2.6 20, 12.6 20, 12.6 12.5, 10.1 10, 5.4 8.4, 0.5 10.7))"
assert(expected.equals(actual))
}
+
+ it("Should pass ST_SharedPaths") {
+ val result = sparkSession
+ .sql("""
+ |SELECT ST_AsText(ST_SharedPaths(
+ | ST_GeomFromWKT('LINESTRING (0 0, 10 0)'),
+ | ST_GeomFromWKT('LINESTRING (5 0, 15 0)')))
+ |""".stripMargin)
+ .first()
+ .getString(0)
+
+ assertEquals(
+ "GEOMETRYCOLLECTION (MULTILINESTRING ((5 0, 10 0)), MULTILINESTRING EMPTY)",
+ result)
+
+ val dimensionalResult = sparkSession
+ .sql("""
+ |WITH shared AS (
+ | SELECT ST_SharedPaths(
+ | ST_SetSRID(ST_GeomFromWKT('LINESTRING ZM (0 1 5 4, 0 0 6 5, 1 0 7 6)'), 4326),
+ | ST_SetSRID(ST_GeomFromWKT('LINESTRING M (0 -1 8, 0 0 9, 1 0 10)'), 4326)
+ | ) AS geom
+ |)
+ |SELECT
+ | ST_SRID(geom),
+ | ST_HasZ(ST_GeometryN(ST_GeometryN(geom, 0), 0)),
+ | ST_HasM(ST_GeometryN(ST_GeometryN(geom, 0), 0))
+ |FROM shared
+ |""".stripMargin)
+ .first()
+
+ assertEquals(4326, dimensionalResult.getInt(0))
+ assertTrue(dimensionalResult.getBoolean(1))
+ assertFalse(dimensionalResult.getBoolean(2))
+ }
}
it("Should pass ST_ForcePolygonCW") {