This is an automated email from the ASF dual-hosted git repository. afs pushed a commit to branch main in repository https://gitbox.apache.org/repos/asf/jena.git
commit f784ca447a6feef699e790205bf16424a0055ae5 Author: Edmond Chuc <[email protected]> AuthorDate: Tue Sep 29 14:10:11 2026 +1000 GH-4265: Support geographic CRS area calculations --- .../geosparql/implementation/GeographicArea.java | 267 +++++++++++++++++++++ .../geosparql/implementation/GeometryArea.java | 58 ++++- .../geosparql/implementation/GeometryWrapper.java | 20 +- .../filter_functions/AreaFFTest.java | 58 ++++- .../filter_functions/MetricAreaFFTest.java | 39 ++- .../geosparql/implementation/GeometryAreaTest.java | 203 ++++++++++++++-- 6 files changed, 588 insertions(+), 57 deletions(-) diff --git a/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeographicArea.java b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeographicArea.java new file mode 100644 index 0000000000..315f1b604f --- /dev/null +++ b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeographicArea.java @@ -0,0 +1,267 @@ +/* + * 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 + * + * https://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. + * + * SPDX-License-Identifier: Apache-2.0 + */ +package org.apache.jena.geosparql.implementation; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.Map; + +import org.apache.sis.geometry.DirectPosition2D; +import org.apache.sis.measure.Units; +import org.apache.sis.referencing.CRS; +import org.apache.sis.referencing.CommonCRS; +import org.apache.sis.referencing.GeodeticCalculator; +import org.apache.sis.referencing.GeodeticException; +import org.apache.sis.referencing.crs.AbstractCRS; +import org.apache.sis.referencing.crs.DefaultProjectedCRS; +import org.apache.sis.referencing.cs.AxesConvention; +import org.apache.sis.referencing.operation.DefaultConversion; +import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.LinearRing; +import org.locationtech.jts.geom.Polygon; +import org.opengis.geometry.DirectPosition; +import org.opengis.parameter.ParameterValueGroup; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.crs.GeographicCRS; +import org.opengis.referencing.cs.CartesianCS; +import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.OperationMethod; +import org.opengis.referencing.operation.TransformException; +import org.opengis.util.FactoryException; + +/** Approximates ellipsoidal polygon area in square metres without changing the input geometry. */ +final class GeographicArea { + // These limits bound the projection and boundary approximation. They are not + // presented as an absolute area-error guarantee, which also depends on geometry. + private static final double MAX_CENTER_ANGLE_DEGREES = 89.0; + private static final double MAX_CHORD_ERROR_METRES = 0.25; + private static final double MAX_SEGMENT_METRES = 50_000.0; + private static final int MAX_RECURSION_DEPTH = 24; + private static final int MAX_COORDINATES = 50_000; + + private GeographicArea() { + } + + static double calculate(GeometryWrapper geometry) { + Polygon polygon = (Polygon) geometry.getParsingGeometry(); + CoordinateReferenceSystem horizontal = CRS.getHorizontalComponent(geometry.getSrsInfo().getCrs()); + if (!(horizontal instanceof GeographicCRS sourceCrs) + || sourceCrs.getCoordinateSystem().getDimension() != 2) { + throw new UnitsConversionException("Geographic area requires a two-dimensional geographic horizontal CRS."); + } + + try { + GeographicCRS normalizedCrs = (GeographicCRS) AbstractCRS.castOrCopy(sourceCrs) + .forConvention(AxesConvention.NORMALIZED); + MathTransform toNormalized = CRS.findOperation(sourceCrs, normalizedCrs, null).getMathTransform(); + double[] center = center(polygon, toNormalized); + + DefaultMathTransformFactory factory = DefaultMathTransformFactory.provider(); + OperationMethod method = factory.getOperationMethod("Lambert Azimuthal Equal Area"); + ParameterValueGroup parameters = method.getParameters().createValue(); + parameters.parameter("latitude_of_center").setValue(center[1]); + parameters.parameter("longitude_of_center").setValue(center[0]); + DefaultConversion conversion = new DefaultConversion(Map.of("name", "Local equal-area conversion"), + method, null, parameters); + CartesianCS cartesian = (CartesianCS) CommonCRS.WGS84.universal(0, 0).getCoordinateSystem(); + DefaultProjectedCRS projectedCrs = new DefaultProjectedCRS(Map.of("name", "Local equal-area CRS"), + sourceCrs, conversion, cartesian); + MathTransform projection = CRS.findOperation(sourceCrs, projectedCrs, null).getMathTransform(); + + CoordinateBudget budget = new CoordinateBudget(); + GeodeticCalculator calculator = GeodeticCalculator.create(sourceCrs); + return projectPolygon(polygon, sourceCrs, projection, calculator, budget).getArea(); + } catch (FactoryException | TransformException | GeodeticException e) { + throw new UnitsConversionException("Geographic area could not be calculated reliably.", e); + } + } + + private static double[] center(Polygon polygon, MathTransform toNormalized) throws TransformException { + List<double[]> vertices = new ArrayList<>(); + List<Double> longitudes = new ArrayList<>(); + double[] sum = new double[3]; + double minLatitude = Double.POSITIVE_INFINITY; + double maxLatitude = Double.NEGATIVE_INFINITY; + Coordinate[] shell = polygon.getExteriorRing().getCoordinates(); + for (int i = 0; i < shell.length - 1; i++) { + double[] lonLat = transform(shell[i], toNormalized); + double lon = Math.toRadians(lonLat[0]); + double lat = Math.toRadians(lonLat[1]); + if (!Double.isFinite(lon) || !Double.isFinite(lat) || Math.abs(lat) > Math.PI / 2) { + throw new UnitsConversionException("Geographic area requires finite coordinates within latitude limits."); + } + vertices.add(new double[] {lon, lat}); + if (vertices.size() > MAX_COORDINATES) { + throw new UnitsConversionException("Geographic area exceeds the supported edge approximation limit."); + } + longitudes.add(Math.atan2(Math.sin(lon), Math.cos(lon))); + minLatitude = Math.min(minLatitude, lat); + maxLatitude = Math.max(maxLatitude, lat); + sum[0] += Math.cos(lat) * Math.cos(lon); + sum[1] += Math.cos(lat) * Math.sin(lon); + sum[2] += Math.sin(lat); + } + Collections.sort(longitudes); + double largestGap = -1; + double arcStart = 0; + for (int i = 0; i < longitudes.size(); i++) { + double current = longitudes.get(i); + double next = i + 1 < longitudes.size() ? longitudes.get(i + 1) + : longitudes.get(0) + 2 * Math.PI; + if (next - current > largestGap) { + largestGap = next - current; + arcStart = next; + } + } + double centerLon = arcStart + (2 * Math.PI - largestGap) / 2; + double centerLat = (minLatitude + maxLatitude) / 2; + double bestMinimumCosine = minimumCosine(vertices, centerLon, centerLat); + double norm = Math.hypot(Math.hypot(sum[0], sum[1]), sum[2]); + if (norm >= 1e-9) { + double meanLon = Math.atan2(sum[1], sum[0]); + double meanLat = Math.atan2(sum[2], Math.hypot(sum[0], sum[1])); + double meanMinimumCosine = minimumCosine(vertices, meanLon, meanLat); + if (meanMinimumCosine > bestMinimumCosine) { + centerLon = meanLon; + centerLat = meanLat; + bestMinimumCosine = meanMinimumCosine; + } + } + double maxAngleCosine = Math.cos(Math.toRadians(MAX_CENTER_ANGLE_DEGREES)); + if (bestMinimumCosine < maxAngleCosine) { + throw new UnitsConversionException("Geographic area extends beyond one supported local hemisphere."); + } + return new double[] {Math.toDegrees(Math.atan2(Math.sin(centerLon), Math.cos(centerLon))), + Math.toDegrees(centerLat)}; + } + + private static double minimumCosine(List<double[]> vertices, double centerLon, double centerLat) { + double minimum = 1; + for (double[] vertex : vertices) { + double lon = vertex[0]; + double lat = vertex[1]; + double cosine = Math.sin(centerLat) * Math.sin(lat) + + Math.cos(centerLat) * Math.cos(lat) * Math.cos(lon - centerLon); + minimum = Math.min(minimum, cosine); + } + return minimum; + } + + private static Polygon projectPolygon(Polygon polygon, GeographicCRS sourceCrs, + MathTransform projection, GeodeticCalculator calculator, CoordinateBudget budget) throws TransformException { + GeometryFactory factory = polygon.getFactory(); + LinearRing shell = projectRing(polygon.getExteriorRing(), sourceCrs, projection, calculator, budget); + LinearRing[] holes = new LinearRing[polygon.getNumInteriorRing()]; + for (int i = 0; i < holes.length; i++) { + holes[i] = projectRing(polygon.getInteriorRingN(i), sourceCrs, projection, calculator, budget); + } + return factory.createPolygon(shell, holes); + } + + private static LinearRing projectRing(LinearRing ring, GeographicCRS sourceCrs, + MathTransform projection, GeodeticCalculator calculator, CoordinateBudget budget) throws TransformException { + Coordinate[] source = ring.getCoordinates(); + List<Coordinate> target = new ArrayList<>(); + target.add(project(source[0], projection)); + for (int i = 1; i < source.length; i++) { + Coordinate end = project(source[i], projection); + appendGeodesic(source[i - 1], source[i], target.get(target.size() - 1), end, + sourceCrs, projection, calculator, target, budget, 0); + } + target.set(target.size() - 1, new Coordinate(target.get(0))); + return ring.getFactory().createLinearRing(target.toArray(Coordinate[]::new)); + } + + private static void appendGeodesic(Coordinate start, Coordinate end, Coordinate projectedStart, + Coordinate projectedEnd, GeographicCRS sourceCrs, MathTransform projection, + GeodeticCalculator calculator, List<Coordinate> target, CoordinateBudget budget, int depth) throws TransformException { + calculator.setStartPoint(new DirectPosition2D(sourceCrs, start.x, start.y)); + calculator.setEndPoint(new DirectPosition2D(sourceCrs, end.x, end.y)); + double length = calculator.getGeodesicDistance(); + double metres = calculator.getDistanceUnit().getConverterTo(Units.METRE).convert(length); + if (!Double.isFinite(metres)) { + throw new UnitsConversionException("Geographic area has a non-finite polygon edge."); + } + if (metres == 0) { + add(target, projectedEnd, budget); + return; + } + double azimuth = calculator.getStartingAzimuth(); + calculator.setStartingAzimuth(azimuth); + calculator.setGeodesicDistance(length / 2); + DirectPosition middle = calculator.getEndPoint(); + Coordinate midpoint = new Coordinate(middle.getOrdinate(0), middle.getOrdinate(1)); + Coordinate projectedMidpoint = project(midpoint, projection); + double chordError = distanceToSegment(projectedMidpoint, projectedStart, projectedEnd); + if (metres > MAX_SEGMENT_METRES || chordError > MAX_CHORD_ERROR_METRES) { + if (depth >= MAX_RECURSION_DEPTH) { + throw new UnitsConversionException("Geographic area exceeds the supported edge approximation limit."); + } + appendGeodesic(start, midpoint, projectedStart, projectedMidpoint, + sourceCrs, projection, calculator, target, budget, depth + 1); + appendGeodesic(midpoint, end, projectedMidpoint, projectedEnd, + sourceCrs, projection, calculator, target, budget, depth + 1); + } else { + add(target, projectedEnd, budget); + } + } + + private static void add(List<Coordinate> target, Coordinate point, CoordinateBudget budget) { + if (++budget.count > MAX_COORDINATES) { + throw new UnitsConversionException("Geographic area exceeds the supported edge approximation limit."); + } + target.add(point); + } + + private static double distanceToSegment(Coordinate point, Coordinate start, Coordinate end) { + double dx = end.x - start.x; + double dy = end.y - start.y; + double lengthSquared = dx * dx + dy * dy; + if (lengthSquared == 0) { + return point.distance(start); + } + double fraction = Math.max(0, Math.min(1, + ((point.x - start.x) * dx + (point.y - start.y) * dy) / lengthSquared)); + return Math.hypot(point.x - start.x - fraction * dx, + point.y - start.y - fraction * dy); + } + + private static Coordinate project(Coordinate point, MathTransform transform) throws TransformException { + double[] projected = transform(point, transform); + if (!Double.isFinite(projected[0]) || !Double.isFinite(projected[1])) { + throw new UnitsConversionException("Geographic area has a coordinate outside the local projection."); + } + return new Coordinate(projected[0], projected[1]); + } + + private static double[] transform(Coordinate point, MathTransform transform) throws TransformException { + double[] result = new double[2]; + transform.transform(new double[] {point.x, point.y}, 0, result, 0, 1); + return result; + } + + private static final class CoordinateBudget { + int count; + } +} diff --git a/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryArea.java b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryArea.java index 501e9f40d8..3e0cf2e1a1 100644 --- a/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryArea.java +++ b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryArea.java @@ -24,20 +24,26 @@ import javax.measure.IncommensurableException; import javax.measure.Unit; import javax.measure.quantity.Area; +import org.apache.jena.geosparql.configuration.GeoSPARQLConfig; import org.apache.sis.measure.Quantities; +import org.apache.sis.measure.Units; import org.apache.sis.referencing.CRS; +import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Geometry; -import org.locationtech.jts.geom.MultiPolygon; import org.locationtech.jts.geom.Polygon; import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.crs.GeographicCRS; +import org.opengis.referencing.cs.CartesianCS; import org.opengis.referencing.cs.CoordinateSystem; /** - * Calculates planar area for polygonal geometries and converts to target area units. + * Calculates the area of a Polygon and converts to target area units. * - * <p>Area is calculated for {@link Polygon} and {@link MultiPolygon} geometries; - * empty geometries and non-polygonal types return zero. Non-empty polygonal geometries - * require a non-geographic horizontal CRS with equivalent linear units on both axes. + * <p>Empty geometries and types other than {@link Polygon} return zero, including + * collections containing polygons. Non-empty geometries require finite X/Y + * coordinates. For geographic CRSs, area is approximated on the source ellipsoid + * using a local equal-area projection. Other non-empty Polygons require a Cartesian + * horizontal CRS with equivalent linear units on both axes. */ final class GeometryArea { private GeometryArea() { @@ -46,19 +52,44 @@ final class GeometryArea { static double calculate(GeometryWrapper geometry, String targetUnitUri) { Unit<Area> targetUnit = AreaUnitsOfMeasure.getUnit(targetUnitUri); Geometry xyGeometry = geometry.getXYGeometry(); - if (!(xyGeometry instanceof Polygon || xyGeometry instanceof MultiPolygon) - || xyGeometry.isEmpty()) { + for (Coordinate coordinate : xyGeometry.getCoordinates()) { + if (!coordinate.isValid()) { + throw new UnitsConversionException("Area requires finite X/Y coordinates."); + } + } + // GeoSPARQL 1.1 requires zero for every geometry type other than Polygon. + if (!(xyGeometry instanceof Polygon polygon) || polygon.isEmpty()) { return 0.0; } - if (geometry.getSrsInfo().isGeographic()) { - throw new UnitsConversionException("Area is not supported for geographic coordinate reference systems."); + double sourceArea; + Unit<Area> sourceUnit; + if (CRS.getHorizontalComponent(geometry.getSrsInfo().getCrs()) instanceof GeographicCRS) { + if (!geometry.getSrsInfo().isSRSRecognised()) { + throw new UnitsConversionException("Area requires a recognised geographic coordinate reference system."); + } + if (!GeoSPARQLConfig.ALLOW_GEOMETRY_SRS_TRANSFORMATION) { + throw new UnitsConversionException("Geographic area requires geometry SRS transformation to be enabled."); + } + sourceArea = GeographicArea.calculate(geometry); + sourceUnit = Units.SQUARE_METRE; + } else { + sourceUnit = equivalentHorizontalAxisUnits(geometry).getUnit() + .pow(2).asType(Area.class); + sourceArea = polygon.getArea(); } - Unit<Area> sourceUnit = equivalentHorizontalAxisUnits(geometry).getUnit() - .pow(2).asType(Area.class); - return Quantities.create(xyGeometry.getArea(), sourceUnit) + double area = Quantities.create(sourceArea, sourceUnit) .to(targetUnit).getValue().doubleValue(); + if (!Double.isFinite(area)) { + throw new UnitsConversionException("Area result is not finite."); + } + return area; } + /** + * Returns the common linear unit of a two-dimensional Cartesian horizontal CRS. + * JTS calculates area from coordinate values as though the axes were orthogonal. + * This method rejects non-Cartesian axes because it does not correct for skew. + */ private static UnitsOfMeasure equivalentHorizontalAxisUnits(GeometryWrapper geometry) { CoordinateReferenceSystem horizontalCrs = CRS.getHorizontalComponent( geometry.getSrsInfo().getCrs()); @@ -68,6 +99,9 @@ final class GeometryArea { "Area requires a two-dimensional horizontal source coordinate system."); } CoordinateSystem coordinateSystem = horizontalCrs.getCoordinateSystem(); + if (!(coordinateSystem instanceof CartesianCS)) { + throw new UnitsConversionException("Area requires Cartesian horizontal source axes."); + } Unit<?> firstAxisUnit = coordinateSystem.getAxis(0).getUnit(); Unit<?> secondAxisUnit = coordinateSystem.getAxis(1).getUnit(); try { diff --git a/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryWrapper.java b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryWrapper.java index 411d5319e9..3e7c8b70c5 100644 --- a/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryWrapper.java +++ b/jena-geosparql/src/main/java/org/apache/jena/geosparql/implementation/GeometryWrapper.java @@ -810,20 +810,30 @@ public class GeometryWrapper implements Serializable { return new GeometryWrapper(parsingGeo, xyGeo, srsInfo.getSrsURI(), geometryDatatypeURI, dimensionInfo); } - /** Polygon or MultiPolygon area defaulting to square metres. */ + /** Area of a {@link org.locationtech.jts.geom.Polygon}, defaulting to square metres. */ public double area() { return area(Unit_URI.SQUARE_METRE_QUDT); } /** - * Returns Polygon or MultiPolygon area in the requested area unit. + * Returns the area of a {@link org.locationtech.jts.geom.Polygon} in the requested area unit. * - * Empty and non-polygonal geometries return zero. Nonempty polygons require - * a non-geographic horizontal CRS with equivalent linear units on both axes. + * <p>Polygon holes are subtracted. Lines have zero area even when their + * coordinates close to form a ring. + * + * <p>Empty geometries and types other than Polygon return zero, including + * collections containing polygons. Geographic polygon edges follow shortest + * geodesics on the source ellipsoid; their area is approximated with a local + * equal-area projection. The geographic CRS must be recognised and the + * geometry SRS transformation setting must be enabled. Polygons too large + * for a reliable local projection produce an error. Other nonempty polygons + * require a Cartesian horizontal CRS with equivalent linear units on both axes. + * Nonempty geometries with non-finite X/Y coordinates produce an error. * * @param unitsURI URI of an explicit area unit. * @throws org.apache.jena.geosparql.implementation.registry.UnitsURIException if the unit URI is unknown. - * @throws UnitsConversionException if the source CRS or target units are unsupported. + * @throws UnitsConversionException if the source CRS or target units are unsupported, + * or the geometry's X/Y coordinates or calculated area are non-finite. */ public double area(String unitsURI) { return GeometryArea.calculate(this, unitsURI); diff --git a/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/AreaFFTest.java b/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/AreaFFTest.java index 267a4386e3..924e7ee60c 100644 --- a/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/AreaFFTest.java +++ b/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/AreaFFTest.java @@ -23,8 +23,10 @@ package org.apache.jena.geosparql.geof.nontopological.filter_functions; import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertNull; import static org.junit.Assert.assertThrows; +import static org.junit.Assert.assertTrue; import org.apache.jena.geosparql.configuration.GeoSPARQLConfig; +import org.apache.jena.geosparql.implementation.datatype.GMLDatatype; import org.apache.jena.geosparql.implementation.datatype.WKTDatatype; import org.apache.jena.geosparql.implementation.vocabulary.Unit_URI; import org.apache.jena.graph.Node; @@ -53,19 +55,20 @@ public class AreaFFTest { } @Test - public void nonPolygonalAndEmptyGeometriesReturnZero() { - for (String wkt : new String[] { "POINT (1 2)", PROJECTED + "POLYGON EMPTY" }) { + public void nonPolygonAndEmptyGeometriesReturnZero() { + for (String wkt : new String[] { "POINT (1 2)", PROJECTED + "POLYGON EMPTY", + PROJECTED + "MULTIPOLYGON(((0 0,2 0,2 2,0 2,0 0)))", + "GEOMETRYCOLLECTION(POLYGON((0 0,2 0,2 2,0 2,0 0)))" }) { assertEquals(NodeValue.makeDouble(0).asNode(), evaluate("'" + wkt + "'^^geo:wktLiteral", SQUARE_METRE)); } } @Test - public void geographicPolygonRaisesExpressionError() { - String wkt = "POLYGON ((0 0, 1 0, 1 1, 0 0))"; - assertNull(evaluate("'" + wkt + "'^^geo:wktLiteral", SQUARE_METRE)); - assertThrows(ExprEvalException.class, - () -> function.exec(NodeValue.makeNode(wkt, WKTDatatype.INSTANCE), - NodeValue.makeNode(NodeFactory.createURI(Unit_URI.SQUARE_METRE_QUDT)))); + public void geographicPolygonReturnsArea() { + String wkt = "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"; + Node result = evaluate("'" + wkt + "'^^geo:wktLiteral", SQUARE_METRE); + double area = ((Number) result.getLiteralValue()).doubleValue(); + assertTrue(area > 12_000_000_000.0 && area < 12_500_000_000.0); } @Test @@ -79,6 +82,22 @@ public class AreaFFTest { assertEquals(NodeValue.makeDouble(12).asNode(), evaluate(gml, SQUARE_METRE)); } + @Test + public void geographicGmlUsesDeclaredAxisOrder() { + String wkt = "'<http://www.opengis.net/def/crs/EPSG/0/4326> " + + "POLYGON((0 0,0 2,1 2,1 0,0 0))'^^geo:wktLiteral"; + String gml = """ + '<gml:Polygon xmlns:gml="http://www.opengis.net/gml/3.2"; + srsName="http://www.opengis.net/def/crs/EPSG/0/4326";> + <gml:exterior><gml:LinearRing><gml:posList>0 0 0 2 1 2 1 0 0 0</gml:posList></gml:LinearRing></gml:exterior> + </gml:Polygon>'^^geo:gmlLiteral + """.replace("\n", " "); + + double expected = ((Number) evaluate(wkt, SQUARE_METRE).getLiteralValue()).doubleValue(); + double actual = ((Number) evaluate(gml, SQUARE_METRE).getLiteralValue()).doubleValue(); + assertEquals(expected, actual, expected * 1e-8); + } + @Test public void malformedOrUnboundGeometryLeavesBindUnbound() { for (String value : new String[] { "42", "'POINT (1 2)'", "<urn:geometry>", "'invalid'^^geo:wktLiteral", "?missing" }) { @@ -95,6 +114,29 @@ public class AreaFFTest { } } + @Test + public void nonFiniteCoordinatesAndAreaRaiseExpressionErrors() { + NodeValue squareMetre = NodeValue.makeNode(NodeFactory.createURI(Unit_URI.SQUARE_METRE_QUDT)); + for (String wkt : new String[] { + PROJECTED + "LINESTRING(0 0,1e309 1)", + PROJECTED + "POLYGON((0 0,1 0,1e309 1,0 1,0 0))", + PROJECTED + "POLYGON((0 0,1e200 0,1e200 1e200,0 1e200,0 0))" }) { + assertThrows(wkt, ExprEvalException.class, + () -> function.exec(NodeValue.makeNode(wkt, WKTDatatype.INSTANCE), squareMetre)); + assertNull(wkt, evaluate("'" + wkt + "'^^geo:wktLiteral", SQUARE_METRE)); + } + + String gml = """ + <gml:Polygon xmlns:gml="http://www.opengis.net/gml/3.2"; + srsName="http://www.opengis.net/def/crs/EPSG/0/27700";> + <gml:exterior><gml:LinearRing><gml:posList>0 0 1 0 NaN 1 0 1 0 0</gml:posList></gml:LinearRing></gml:exterior> + </gml:Polygon> + """.replace("\n", " "); + assertThrows(ExprEvalException.class, + () -> function.exec(NodeValue.makeNode(gml, GMLDatatype.INSTANCE), squareMetre)); + assertNull(evaluate("'" + gml + "'^^geo:gmlLiteral", SQUARE_METRE)); + } + @Test public void wrongArityIsRejectedAtQueryBuild() { String geometry = "'POINT EMPTY'^^geo:wktLiteral"; diff --git a/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/MetricAreaFFTest.java b/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/MetricAreaFFTest.java index af9019e246..d7b73a60eb 100644 --- a/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/MetricAreaFFTest.java +++ b/jena-geosparql/src/test/java/org/apache/jena/geosparql/geof/nontopological/filter_functions/MetricAreaFFTest.java @@ -23,8 +23,10 @@ package org.apache.jena.geosparql.geof.nontopological.filter_functions; import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertNull; import static org.junit.Assert.assertThrows; +import static org.junit.Assert.assertTrue; import org.apache.jena.geosparql.configuration.GeoSPARQLConfig; +import org.apache.jena.geosparql.implementation.datatype.GMLDatatype; import org.apache.jena.geosparql.implementation.datatype.WKTDatatype; import org.apache.jena.graph.Node; import org.apache.jena.graph.NodeFactory; @@ -50,17 +52,20 @@ public class MetricAreaFFTest { } @Test - public void nonPolygonalAndEmptyGeometriesReturnZero() { - for (String wkt : new String[] { "POINT (1 2)", PROJECTED + "POLYGON EMPTY" }) { + public void nonPolygonAndEmptyGeometriesReturnZero() { + for (String wkt : new String[] { "POINT (1 2)", PROJECTED + "POLYGON EMPTY", + PROJECTED + "MULTIPOLYGON(((0 0,2 0,2 2,0 2,0 0)))", + "GEOMETRYCOLLECTION(POLYGON((0 0,2 0,2 2,0 2,0 0)))" }) { assertEquals(NodeValue.makeDouble(0).asNode(), evaluate("'" + wkt + "'^^geo:wktLiteral")); } } @Test - public void geographicPolygonRaisesExpressionError() { - String wkt = "POLYGON ((0 0, 1 0, 1 1, 0 0))"; - assertNull(evaluate("'" + wkt + "'^^geo:wktLiteral")); - assertThrows(ExprEvalException.class, () -> exec(NodeValue.makeNode(wkt, WKTDatatype.INSTANCE))); + public void geographicPolygonReturnsSquareMetres() { + String wkt = "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"; + Node result = evaluate("'" + wkt + "'^^geo:wktLiteral"); + double area = ((Number) result.getLiteralValue()).doubleValue(); + assertTrue(area > 12_000_000_000.0 && area < 12_500_000_000.0); } @Test @@ -89,6 +94,28 @@ public class MetricAreaFFTest { } } + @Test + public void nonFiniteCoordinatesAndAreaRaiseExpressionErrors() { + for (String wkt : new String[] { + PROJECTED + "LINESTRING(0 0,1e309 1)", + PROJECTED + "POLYGON((0 0,1 0,1e309 1,0 1,0 0))", + PROJECTED + "POLYGON((0 0,1e200 0,1e200 1e200,0 1e200,0 0))" }) { + assertThrows(wkt, ExprEvalException.class, + () -> function.exec(NodeValue.makeNode(wkt, WKTDatatype.INSTANCE))); + assertNull(wkt, evaluate("'" + wkt + "'^^geo:wktLiteral")); + } + + String gml = """ + <gml:Polygon xmlns:gml="http://www.opengis.net/gml/3.2"; + srsName="http://www.opengis.net/def/crs/EPSG/0/27700";> + <gml:exterior><gml:LinearRing><gml:posList>0 0 1 0 NaN 1 0 1 0 0</gml:posList></gml:LinearRing></gml:exterior> + </gml:Polygon> + """.replace("\n", " "); + assertThrows(ExprEvalException.class, + () -> function.exec(NodeValue.makeNode(gml, GMLDatatype.INSTANCE))); + assertNull(evaluate("'" + gml + "'^^geo:gmlLiteral")); + } + @Test public void wrongArityIsRejectedAtQueryBuild() { String geometry = "'POINT EMPTY'^^geo:wktLiteral"; diff --git a/jena-geosparql/src/test/java/org/apache/jena/geosparql/implementation/GeometryAreaTest.java b/jena-geosparql/src/test/java/org/apache/jena/geosparql/implementation/GeometryAreaTest.java index c39434160d..c4c8b7f3f0 100644 --- a/jena-geosparql/src/test/java/org/apache/jena/geosparql/implementation/GeometryAreaTest.java +++ b/jena-geosparql/src/test/java/org/apache/jena/geosparql/implementation/GeometryAreaTest.java @@ -21,6 +21,7 @@ package org.apache.jena.geosparql.implementation; import org.apache.jena.geosparql.configuration.GeoSPARQLConfig; +import org.apache.jena.geosparql.implementation.datatype.GMLDatatype; import org.apache.jena.geosparql.implementation.datatype.WKTDatatype; import org.apache.jena.geosparql.implementation.registry.UnitsURIException; import org.apache.jena.geosparql.implementation.vocabulary.Unit_URI; @@ -33,10 +34,16 @@ import java.util.List; import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertThrows; +import static org.junit.Assert.assertTrue; public class GeometryAreaTest { private static final String EPSG_32634 = "http://www.opengis.net/def/crs/EPSG/0/32634";; private static final String EPSG_2227 = "http://www.opengis.net/def/crs/EPSG/0/2227";; + private static final String EPSG_4326 = "http://www.opengis.net/def/crs/EPSG/0/4326";; + private static final String EPSG_4267 = "http://www.opengis.net/def/crs/EPSG/0/4267";; + private static final String EPSG_4275 = "http://www.opengis.net/def/crs/EPSG/0/4275";; + private static final String EPSG_4807 = "http://www.opengis.net/def/crs/EPSG/0/4807";; + private static final String EPSG_4979 = "http://www.opengis.net/def/crs/EPSG/0/4979";; private static final String UNKNOWN_UNIT = "http://example.com/unit/unknown";; @BeforeClass @@ -54,16 +61,12 @@ public class GeometryAreaTest { } @Test - public void areaAccumulatesMultiPolygonShellsAndHoles() { - GeometryWrapper multiPolygon = geometry("<" + EPSG_32634 + "> MULTIPOLYGON(" - + "((500000 4600000,500010 4600000,500010 4600010," - + "500000 4600010,500000 4600000)," - + "(500004 4600004,500006 4600004,500006 4600006," - + "500004 4600006,500004 4600004))," - + "((500020 4600000,500023 4600000,500023 4600004," - + "500020 4600004,500020 4600000)))"); + public void polygonHolesAreExcluded() { + GeometryWrapper polygon = geometry("<" + EPSG_32634 + "> POLYGON(" + + "(0 0,10 0,10 10,0 10,0 0)," + + "(4 4,6 4,6 6,4 6,4 4))"); - assertEquals(108.0, multiPolygon.area(), 0.0); + assertEquals(96.0, polygon.area(), 0.0); } @Test @@ -86,17 +89,57 @@ public class GeometryAreaTest { } @Test - public void areaIneligibleGeometryTypesReturnZeroWithoutSourceCrsCalculation() { + public void nonPolygonTypesReturnZeroWithoutSourceCrsCalculation() { for (String wkt : List.of( "POINT(1 1)", "LINESTRING(0 0,1 1)", "MULTIPOINT((0 0),(1 1))", "MULTILINESTRING((0 0,1 1),(2 2,3 3))", - "GEOMETRYCOLLECTION(POLYGON((0 0,0 1,1 1,1 0,0 0)))")) { + "MULTIPOLYGON(((0 0,2 0,2 2,0 2,0 0)))", + "GEOMETRYCOLLECTION(POINT(1 1),LINESTRING(0 0,1 1))", + "GEOMETRYCOLLECTION(POLYGON EMPTY,POINT(1 1))", + "GEOMETRYCOLLECTION(POLYGON((0 0,2 0,2 2,0 2,0 0)))")) { assertEquals(0.0, geometry(wkt).area(), 0.0); + assertEquals(0.0, geometry("<" + EPSG_32634 + "> " + wkt).area(), 0.0); } } + @Test + public void nonFiniteXYCoordinatesAreRejected() { + for (String wkt : List.of( + "LINESTRING(0 0,1e309 1)", + "<" + EPSG_32634 + "> POLYGON((0 0,1 0,1e309 1,0 1,0 0))")) { + assertThrows(UnitsConversionException.class, () -> geometry(wkt).area()); + } + + String gml = """ + <gml:Polygon xmlns:gml="http://www.opengis.net/gml/3.2"; + srsName="http://www.opengis.net/def/crs/EPSG/0/32634";> + <gml:exterior><gml:LinearRing><gml:posList>0 0 1 0 NaN 1 0 1 0 0</gml:posList></gml:LinearRing></gml:exterior> + </gml:Polygon> + """; + GeometryWrapper polygon = GeometryWrapper.extract(gml, GMLDatatype.URI); + assertThrows(UnitsConversionException.class, polygon::area); + } + + @Test + public void nonFinitePlanarAreaResultIsRejected() { + GeometryWrapper polygon = geometry("<" + EPSG_32634 + "> POLYGON((" + + "0 0,1e200 0,1e200 1e200,0 1e200,0 0))"); + + assertThrows(UnitsConversionException.class, polygon::area); + } + + @Test + public void areaOverflowDuringUnitConversionIsRejected() { + GeometryWrapper polygon = geometry("<" + EPSG_32634 + "> POLYGON((" + + "0 0,2e151 0,2e151 2e151,0 2e151,0 0))"); + + assertTrue(Double.isFinite(polygon.area())); + assertThrows(UnitsConversionException.class, + () -> polygon.area(Unit_URI.SQUARE_MILLIMETRE_QUDT)); + } + @Test public void eligibleAreaRejectsDifferentHorizontalAxisScales() throws Exception { GeometryWrapper polygon = geometryWithCrs( @@ -110,6 +153,31 @@ public class GeometryAreaTest { polygon::area); } + @Test + public void eligibleAreaRejectsNonCartesianHorizontalAxes() throws Exception { + GeometryWrapper polygon = geometryWithCrs( + "<" + EPSG_32634 + "> POLYGON((0 0,1 0,1 1,0 1,0 0))", + CRS.fromWKT("ENGCRS[\"Skewed axes\"," + + "EDATUM[\"Engineering datum\"],CS[affine,2]," + + "AXIS[\"x\",east,ORDER[1],LENGTHUNIT[\"metre\",1]]," + + "AXIS[\"y\",northEast,ORDER[2],LENGTHUNIT[\"metre\",1]]]")); + + UnitsConversionException error = assertThrows(UnitsConversionException.class, polygon::area); + assertTrue(error.getMessage().contains("Cartesian")); + } + + @Test + public void cartesianEngineeringCrsStillSupportsArea() throws Exception { + GeometryWrapper polygon = geometryWithCrs( + "<" + EPSG_32634 + "> POLYGON((0 0,1 0,1 1,0 1,0 0))", + CRS.fromWKT("ENGCRS[\"Cartesian engineering grid\"," + + "EDATUM[\"Engineering datum\"],CS[Cartesian,2]," + + "AXIS[\"x\",east,ORDER[1],LENGTHUNIT[\"metre\",1]]," + + "AXIS[\"y\",north,ORDER[2],LENGTHUNIT[\"metre\",1]]]")); + + assertEquals(1.0, polygon.area(), 0.0); + } + @Test public void projectedNonMetreAreaConvertsTheSourceAreaQuantity() { GeometryWrapper polygon = geometry("<" + EPSG_2227 + "> POLYGON((" @@ -120,16 +188,6 @@ public class GeometryAreaTest { polygon.area(), 1e-12); } - @Test - public void multiPolygonAreaIsAccumulatedBeforeUnitConversion() { - GeometryWrapper multiPolygon = geometry("<" + EPSG_32634 + "> MULTIPOLYGON(" - + "((0 0,0.02 0,0.02 0.02,0 0.02,0 0))," - + "((1 0,1.02 0,1.02 0.02,1 0.02,1 0)))"); - - assertEquals(0.0000000008, - multiPolygon.area(Unit_URI.SQUARE_KILOMETRE_QUDT), 1e-20); - } - @Test public void compoundCrsUsesItsTwoDimensionalHorizontalAxisUnits() throws Exception { CoordinateReferenceSystem compoundCrs = CRS.compound( @@ -142,11 +200,104 @@ public class GeometryAreaTest { } @Test - public void eligibleGeographicAreaIsRejected() { - for (String wkt : List.of( - "POLYGON((0 0,0 1,1 1,1 0,0 0))", - "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)))")) { - assertThrows(UnitsConversionException.class, geometry(wkt)::area); + public void geographicAreaUsesTheDeclaredAxisOrderAndConvertsUnits() { + GeometryWrapper crs84 = geometry("POLYGON((0 0,2 0,2 1,0 1,0 0))"); + GeometryWrapper epsg4326 = geometry("<" + EPSG_4326 + "> POLYGON((" + + "0 0,0 2,1 2,1 0,0 0))"); + + double squareMetres = crs84.area(); + assertTrue(squareMetres > 24_000_000_000.0 && squareMetres < 25_000_000_000.0); + assertEquals(squareMetres, epsg4326.area(), squareMetres * 1e-8); + assertEquals(squareMetres / 1_000_000, crs84.area(Unit_URI.SQUARE_KILOMETRE_QUDT), 1e-6); + } + + @Test + public void geographicPolygonExcludesHoles() { + String square = "POLYGON((0 0,1 0,1 1,0 1,0 0))"; + GeometryWrapper polygon = geometry(square); + GeometryWrapper withHole = geometry("POLYGON((0 0,1 0,1 1,0 1,0 0)," + + "(0.25 0.25,0.25 0.75,0.75 0.75,0.75 0.25,0.25 0.25))"); + + // WGS84 geodesic reference areas from pyproj.Geod.polygon_area_perimeter. + assertEquals(12_308_778_361.469, polygon.area(), 100_000.0); + assertEquals(9_231_614_224.815, withHole.area(), 100_000.0); + } + + @Test + public void geographicAreaSupportsDatelineAndPolarPolygons() { + double dateline = geometry("POLYGON((179 10,-179 10,-179 11,179 11,179 10))").area(); + double polar = geometry("POLYGON((-90 80,0 80,90 80,180 80,-90 80))").area(); + + // WGS84 geodesic references from pyproj.Geod.polygon_area_perimeter; + // tolerances allow for projected-edge approximation. + assertEquals(24_218_606_783.655, dateline, 250_000.0); + assertEquals(2_507_270_031_169.875, polar, 2_500_000.0); + } + + @Test + public void geographicAreaUsesTheSourceDatum() { + String square = "POLYGON((0 0,1 0,1 1,0 1,0 0))"; + double wgs84 = geometry(square).area(); + double nad27 = geometry("<" + EPSG_4267 + "> " + square).area(); + + assertTrue(nad27 > 0 && Double.isFinite(nad27)); + assertTrue(Math.abs(nad27 - wgs84) > 10_000.0); + } + + @Test + public void parisPrimeMeridianAndGradUnitsMatchGreenwichDegrees() { + // NTF (Paris) uses grads and the Paris meridian; NTF uses degrees and Greenwich. + GeometryWrapper paris = geometry("<" + EPSG_4807 + "> POLYGON((" + + "48 0,48 1,49 1,49 0,48 0))"); + GeometryWrapper greenwich = geometry("<" + EPSG_4275 + "> POLYGON((" + + "43.2 2.33722917,43.2 3.23722917,44.1 3.23722917," + + "44.1 2.33722917,43.2 2.33722917))"); + + double area = greenwich.area(); + assertTrue(area > 0); + assertEquals(area, paris.area(), area * 1e-5); + } + + @Test + public void geographicThreeDimensionalCrsIgnoresHeight() { + GeometryWrapper twoDimensional = geometry("<" + EPSG_4326 + "> POLYGON((" + + "0 0,0 2,1 2,1 0,0 0))"); + GeometryWrapper threeDimensional = geometry("<" + EPSG_4979 + "> POLYGON Z ((" + + "0 0 10,0 2 20,1 2 30,1 0 40,0 0 10))"); + + assertEquals(twoDimensional.area(), threeDimensional.area(), twoDimensional.area() * 1e-8); + } + + @Test + public void compoundGeographicCrsUsesItsHorizontalComponent() throws Exception { + CoordinateReferenceSystem compoundCrs = CRS.compound( + CRS.forCode(EPSG_4326), CRS.forCode("EPSG:5703")); + GeometryWrapper polygon = geometryWithCrs("<" + EPSG_32634 + "> POLYGON Z ((" + + "0 0 5,0 1 5,1 1 5,1 0 5,0 0 5))", compoundCrs); + + assertEquals(geometry("POLYGON((0 0,1 0,1 1,0 1,0 0))").area(), + polygon.area(), 100_000.0); + } + + @Test + public void geographicAreaRejectsPolygonBeyondOneLocalHemisphere() { + GeometryWrapper wide = geometry("POLYGON((0 0,179.9 0,179.9 1,0 1,0 0))"); + + assertThrows(UnitsConversionException.class, wide::area); + } + + @Test + public void geographicAreaRejectsUnknownCrsAndDisabledTransformation() { + GeometryWrapper unknown = geometry("<http://example.com/crs/unknown> " + + "POLYGON((0 0,1 0,1 1,0 1,0 0))"); + assertThrows(UnitsConversionException.class, unknown::area); + + GeoSPARQLConfig.allowGeometrySRSTransformation(false); + try { + assertThrows(UnitsConversionException.class, + () -> geometry("POLYGON((0 0,1 0,1 1,0 1,0 0))").area()); + } finally { + GeoSPARQLConfig.allowGeometrySRSTransformation(true); } }
