From 681112ac25bca8d0afebb8f2a8926a1a1c8eb46e Mon Sep 17 00:00:00 2001 From: Mike Barry Date: Fri, 25 Oct 2024 07:06:56 -0400 Subject: [PATCH 1/3] size in meters --- .../onthegomap/planetiler/geo/GeoUtils.java | 91 ++++++++++++++++++- .../planetiler/geo/GeoUtilsTest.java | 34 +++++++ 2 files changed, 123 insertions(+), 2 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 f0e0c13189..7b07ff7178 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,91 @@ public static WKTReader wktReader() { 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(LinearRing ring) { + double total = 0; + var cs = ring.getCoordinateSequence(); + var numEdges = cs.size() - 1; + for (int i = 0; i < numEdges; i++) { + double lowerX = cs.getX(i) * RADIANS_PER_DEGREE; + double midY = cs.getY(i + 1 == numEdges ? 0 : i + 1) * RADIANS_PER_DEGREE; + double upperX = cs.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(LinearRing)}. + */ + public static double areaInMeters(Geometry latLonGeom) { + return switch (latLonGeom) { + case Polygon poly -> { + double result = ringAreaMeters(poly.getExteriorRing()); + for (int i = 0; i < poly.getNumInteriorRing(); i++) { + result -= ringAreaMeters(poly.getInteriorRingN(i)); + } + 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/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java b/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java index 64809ec0e7..3f799c255a 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 @@ void getLongestLine() { 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); + } } From 120c92eeb36db76d06da6ef765e27565aa91cce9 Mon Sep 17 00:00:00 2001 From: Mike Barry Date: Sat, 26 Oct 2024 07:08:08 -0400 Subject: [PATCH 2/3] accessor on feature --- .../planetiler/reader/SourceFeature.java | 27 ++++++++++++--- .../planetiler/FeatureCollectorTest.java | 33 +++++++++++++++++++ 2 files changed, 55 insertions(+), 5 deletions(-) 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 10a04a6a7f..f2d51fb197 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 final Geometry validatedPolygon() throws GeometryException { * {@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 double length() throws GeometryException { } /** - * 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 799e2ae81f..8fa6e39449 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 @@ void testPointAlongLine() { 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()); + } } From 2341f8c1fb4abeabcca80788c59ddf8258ea2b3c Mon Sep 17 00:00:00 2001 From: Mike Barry Date: Sat, 26 Oct 2024 07:26:20 -0400 Subject: [PATCH 3/3] change to coordinate sequence --- .../com/onthegomap/planetiler/geo/GeoUtils.java | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 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 7b07ff7178..7b71469277 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 @@ -607,14 +607,13 @@ public static double lineLengthMeters(CoordinateSequence sequence) { * @see "Some Algorithms for Polygons on a Sphere", JPL * Publication 07-03, Jet Propulsion * Laboratory, Pasadena, CA, June 2007. */ - public static double ringAreaMeters(LinearRing ring) { + public static double ringAreaMeters(CoordinateSequence ring) { double total = 0; - var cs = ring.getCoordinateSequence(); - var numEdges = cs.size() - 1; + var numEdges = ring.size() - 1; for (int i = 0; i < numEdges; i++) { - double lowerX = cs.getX(i) * RADIANS_PER_DEGREE; - double midY = cs.getY(i + 1 == numEdges ? 0 : i + 1) * RADIANS_PER_DEGREE; - double upperX = cs.getX(i + 2 >= numEdges ? (i + 2) % numEdges : i + 2) * RADIANS_PER_DEGREE; + 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; @@ -622,14 +621,14 @@ public static double ringAreaMeters(LinearRing ring) { /** * Returns the approximate area in meters of a polygon, or all polygons contained within a multigeometry using - * {@link #ringAreaMeters(LinearRing)}. + * {@link #ringAreaMeters(CoordinateSequence)}. */ public static double areaInMeters(Geometry latLonGeom) { return switch (latLonGeom) { case Polygon poly -> { - double result = ringAreaMeters(poly.getExteriorRing()); + double result = ringAreaMeters(poly.getExteriorRing().getCoordinateSequence()); for (int i = 0; i < poly.getNumInteriorRing(); i++) { - result -= ringAreaMeters(poly.getInteriorRingN(i)); + result -= ringAreaMeters(poly.getInteriorRingN(i).getCoordinateSequence()); } yield result; }