From 431d2a2ae0b67aeb98f8d550165ae04028a2b9ec Mon Sep 17 00:00:00 2001 From: Michael Barry Date: Sat, 26 Oct 2024 08:15:17 -0400 Subject: [PATCH] Add utilities for computing feature size in meters (#1078) --- .../onthegomap/planetiler/geo/GeoUtils.java | 90 ++++++++++++++++++- .../planetiler/reader/SourceFeature.java | 27 ++++-- .../planetiler/FeatureCollectorTest.java | 33 +++++++ .../planetiler/geo/GeoUtilsTest.java | 34 +++++++ 4 files changed, 177 insertions(+), 7 deletions(-) diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java index f0e0c131..7b714692 100644 --- a/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java @@ -52,11 +52,13 @@ public class GeoUtils { private static final LineString[] EMPTY_LINE_STRING_ARRAY = new LineString[0]; private static final Polygon[] EMPTY_POLYGON_ARRAY = new Polygon[0]; private static final Point[] EMPTY_POINT_ARRAY = new Point[0]; - private static final double WORLD_RADIUS_METERS = 6_378_137; - public static final double WORLD_CIRCUMFERENCE_METERS = Math.PI * 2 * WORLD_RADIUS_METERS; + private static final double WORLD_RADIUS_METERS_AT_EQUATOR = 6_378_137; + private static final double AVERAGE_WORLD_RADIUS_METERS = 6_371_008.8; + public static final double WORLD_CIRCUMFERENCE_METERS = Math.PI * 2 * WORLD_RADIUS_METERS_AT_EQUATOR; private static final double RADIANS_PER_DEGREE = Math.PI / 180; private static final double DEGREES_PER_RADIAN = 180 / Math.PI; private static final double LOG2 = Math.log(2); + private static final double AREA_FACTOR = AVERAGE_WORLD_RADIUS_METERS * AVERAGE_WORLD_RADIUS_METERS / 2; /** * Transform web mercator coordinates where top-left corner of the planet is (0,0) and bottom-right is (1,1) to * latitude/longitude coordinates. @@ -575,6 +577,90 @@ public class GeoUtils { return new WKTReader(JTS_FACTORY); } + /** Returns the distance in meters between 2 lat/lon coordinates using the haversine formula. */ + public static double metersBetween(double fromLon, double fromLat, double toLon, double toLat) { + double sinDeltaLat = Math.sin((toLat - fromLat) * RADIANS_PER_DEGREE / 2); + double sinDeltaLon = Math.sin((toLon - fromLon) * RADIANS_PER_DEGREE / 2); + double a = sinDeltaLat * sinDeltaLat + + sinDeltaLon * sinDeltaLon * Math.cos(fromLat * RADIANS_PER_DEGREE) * Math.cos(toLat * RADIANS_PER_DEGREE); + return AVERAGE_WORLD_RADIUS_METERS * 2 * Math.asin(Math.sqrt(a)); + } + + /** Returns the sum of the length of all edges using {@link #metersBetween(double, double, double, double)}. */ + public static double lineLengthMeters(CoordinateSequence sequence) { + double total = 0; + int numEdges = sequence.size() - 1; + for (int i = 0; i < numEdges; i++) { + double fromLon = sequence.getX(i); + double toLon = sequence.getX(i + 1); + double fromLat = sequence.getY(i); + double toLat = sequence.getY(i + 1); + total += metersBetween(fromLon, fromLat, toLon, toLat); + } + + return total; + } + + /** + * Returns the approximate area in meters of a polygon in lat/lon degree coordinates. + * + * @see "Some Algorithms for Polygons on a Sphere", JPL + * Publication 07-03, Jet Propulsion * Laboratory, Pasadena, CA, June 2007. + */ + public static double ringAreaMeters(CoordinateSequence ring) { + double total = 0; + var numEdges = ring.size() - 1; + for (int i = 0; i < numEdges; i++) { + double lowerX = ring.getX(i) * RADIANS_PER_DEGREE; + double midY = ring.getY(i + 1 == numEdges ? 0 : i + 1) * RADIANS_PER_DEGREE; + double upperX = ring.getX(i + 2 >= numEdges ? (i + 2) % numEdges : i + 2) * RADIANS_PER_DEGREE; + total += (upperX - lowerX) * Math.sin(midY); + } + return Math.abs(total) * AREA_FACTOR; + } + + /** + * Returns the approximate area in meters of a polygon, or all polygons contained within a multigeometry using + * {@link #ringAreaMeters(CoordinateSequence)}. + */ + public static double areaInMeters(Geometry latLonGeom) { + return switch (latLonGeom) { + case Polygon poly -> { + double result = ringAreaMeters(poly.getExteriorRing().getCoordinateSequence()); + for (int i = 0; i < poly.getNumInteriorRing(); i++) { + result -= ringAreaMeters(poly.getInteriorRingN(i).getCoordinateSequence()); + } + yield result; + } + case GeometryCollection collection -> { + double result = 0; + for (int i = 0; i < collection.getNumGeometries(); i++) { + result += areaInMeters(collection.getGeometryN(i)); + } + yield result; + } + case null, default -> 0; + }; + } + + /** + * Returns the approximate length in meters of a line, or all lines contained within a multigeometry using + * {@link #lineLengthMeters(CoordinateSequence)}. + */ + public static double lengthInMeters(Geometry latLonGeom) { + return switch (latLonGeom) { + case LineString line -> lineLengthMeters(line.getCoordinateSequence()); + case GeometryCollection collection -> { + double result = 0; + for (int i = 0; i < collection.getNumGeometries(); i++) { + result += lengthInMeters(collection.getGeometryN(i)); + } + yield result; + } + case null, default -> 0; + }; + } + /** Helper class to sort polygons by area of their outer shell. */ private record PolyAndArea(Polygon poly, double area) implements Comparable { diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/reader/SourceFeature.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/reader/SourceFeature.java index 10a04a6a..f2d51fb1 100644 --- a/planetiler-core/src/main/java/com/onthegomap/planetiler/reader/SourceFeature.java +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/reader/SourceFeature.java @@ -43,7 +43,8 @@ public abstract class SourceFeature implements WithTags, WithGeometryType, WithS private Geometry validPolygon = null; private double area = Double.NaN; private double length = Double.NaN; - private double size = Double.NaN; + private double areaMeters = Double.NaN; + private double lengthMeters = Double.NaN; private LineSplitter lineSplitter; /** @@ -283,7 +284,7 @@ public abstract class SourceFeature implements WithTags, WithGeometryType, WithS * {@code 1} means the area of the entire planet. */ public double area() throws GeometryException { - return Double.isNaN(area) ? (area = canBePolygon() ? polygon().getArea() : 0) : area; + return Double.isNaN(area) ? (area = canBePolygon() ? Math.abs(polygon().getArea()) : 0) : area; } /** @@ -297,11 +298,27 @@ public abstract class SourceFeature implements WithTags, WithGeometryType, WithS } /** - * Returns and caches sqrt of {@link #area()} if polygon or {@link #length()} if a line string. + * Returns the sqrt of {@link #area()} if polygon or {@link #length()} if a line string. */ public double size() throws GeometryException { - return Double.isNaN(size) ? (size = canBePolygon() ? Math.sqrt(Math.abs(area())) : canBeLine() ? length() : 0) : - size; + return canBePolygon() ? Math.sqrt(Math.abs(area())) : canBeLine() ? length() : 0; + } + + /** Returns and caches the approximate area of the geometry in square meters. */ + public double areaMeters() throws GeometryException { + return Double.isNaN(areaMeters) ? (areaMeters = + (isPoint() || canBePolygon() || canBeLine()) ? GeoUtils.areaInMeters(latLonGeometry()) : 0) : areaMeters; + } + + /** Returns and caches the approximate length of the geometry in meters. */ + public double lengthMeters() throws GeometryException { + return Double.isNaN(lengthMeters) ? (lengthMeters = + (isPoint() || canBePolygon() || canBeLine()) ? GeoUtils.lengthInMeters(latLonGeometry()) : 0) : lengthMeters; + } + + /** Returns the sqrt of {@link #areaMeters()} if polygon or {@link #lengthMeters()} if a line string. */ + public double sizeMeters() throws GeometryException { + return canBePolygon() ? Math.sqrt(Math.abs(areaMeters())) : canBeLine() ? lengthMeters() : 0; } /** Returns the ID of the source that this feature came from. */ diff --git a/planetiler-core/src/test/java/com/onthegomap/planetiler/FeatureCollectorTest.java b/planetiler-core/src/test/java/com/onthegomap/planetiler/FeatureCollectorTest.java index 799e2ae8..8fa6e394 100644 --- a/planetiler-core/src/test/java/com/onthegomap/planetiler/FeatureCollectorTest.java +++ b/planetiler-core/src/test/java/com/onthegomap/planetiler/FeatureCollectorTest.java @@ -890,4 +890,37 @@ class FeatureCollectorTest { assertFalse(iter.hasNext()); } + + + @Test + void testSizeInMetersOfLine() throws GeometryException { + var sourceLine = newReaderFeature(newLineString( + 0, 0, + 1, 0 + ), Map.of()); + + assertEquals(111_195, sourceLine.lengthMeters(), 1d); + assertEquals(111_195, sourceLine.sizeMeters(), 1d); + assertEquals(0, sourceLine.areaMeters()); + } + + + @Test + void testSizeInMetersOfPolygon() throws GeometryException { + var sourceLine = newReaderFeature(rectangle(0, 1), Map.of()); + + assertEquals(0, sourceLine.lengthMeters()); + assertEquals(111192, sourceLine.sizeMeters(), 1d); + assertEquals(Math.pow(111192.25757749, 2), sourceLine.areaMeters(), 1d); + } + + + @Test + void testSizeInMetersOfPoint() throws GeometryException { + var sourceLine = newReaderFeature(newPoint(0, 1), Map.of()); + + assertEquals(0, sourceLine.lengthMeters()); + assertEquals(0, sourceLine.sizeMeters()); + assertEquals(0, sourceLine.areaMeters()); + } } diff --git a/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java b/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java index 64809ec0..3f799c25 100644 --- a/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java +++ b/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java @@ -457,4 +457,38 @@ class GeoUtilsTest { assertEquals(line1, GeoUtils.getLongestLine(newMultiLineString(line1))); assertEquals(line2, GeoUtils.getLongestLine(newMultiLineString(line1, line2))); } + + + @Test + void areaInMeters() { + assertEquals(0, GeoUtils.areaInMeters(newLineString(0, 0, 1, 1))); + assertEquals(0, GeoUtils.areaInMeters(newPoint(0, 1))); + assertEquals(111_194, Math.sqrt(GeoUtils.areaInMeters(rectangle(-0.5, 0.5))), 1d); + assertEquals(GeoUtils.areaInMeters(rectangle(0, 3)) - GeoUtils.areaInMeters(rectangle(1, 2)), + GeoUtils.areaInMeters(newPolygon(rectangleCoordList(0, 3), List.of(rectangleCoordList(1, 2)))), 1d); + assertEquals(96_963, Math.sqrt(GeoUtils.areaInMeters(rectangle(40, 41))), 1d); + assertEquals(10_387, Math.sqrt(GeoUtils.areaInMeters(rectangle(89, 90))), 1d); + assertEquals(10_387, Math.sqrt(GeoUtils.areaInMeters(rectangle(-89, -90))), 1d); + assertEquals(5.1e14, + GeoUtils.areaInMeters(newPolygon(-180, -90, 180, -90, 180, 90, -180, 90, -180, -90)), 1e11); + } + + @Test + void lengthInMeters() { + assertEquals(0, GeoUtils.lengthInMeters(rectangle(0, 1))); + assertEquals(0, GeoUtils.lengthInMeters(newPoint(0, 0))); + assertEquals(97_129, GeoUtils.lengthInMeters(newLineString( + -75.343, 39.984, + -75.534, 39.123 + )), 1); + assertEquals(97_129 * 3, GeoUtils.lengthInMeters(newLineString( + -75.343, 39.984, + -75.534, 39.123, + -75.343, 39.984, + -75.534, 39.123 + )), 1); + assertEquals(15051, GeoUtils.lengthInMeters(newLineString( + 47.234, 24.235, + 47.234 + 0.1, 24.235 - 0.1)), 1d); + } }