Add utilities for computing feature size in meters (#1078)

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 <a href="https://trs.jpl.nasa.gov/handle/2014/40409">"Some Algorithms for Polygons on a Sphere", JPL
+   *      Publication 07-03, Jet Propulsion * Laboratory, Pasadena, CA, June 2007</a>.
+   */
+  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<PolyAndArea> {
 

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. */

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());
+  }
 }

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);
+  }
 }