Merge branch 'polygon-area-145' into 'main'

Compute polygon area

Closes #145

See merge request !255

src/main/java/gov/usgs/earthquake/nshmp/geo/Locations.java

@@ -15,6 +15,7 @@ import static java.lang.Math.max;
 import static java.lang.Math.min;
 import static java.lang.Math.sin;
 import static java.lang.Math.sqrt;
+import static java.lang.Math.tan;
 
 import java.awt.geom.Line2D;
 import java.awt.geom.Path2D;
@@ -56,6 +57,55 @@ public final class Locations {
 
   private Locations() {}
 
+  /**
+   * Compute the signed area of a closed path in km². Method returns zero if
+   * supplied border has less than three coordinates. Supplied polygon does not
+   * need to be closed (the first and last points do not need to be the same).
+   * Calculation uses {@link Coordinates#EARTH_RADIUS_MEAN} to compute area and
+   * does not consider depths of border coordinates at this time.
+   *
+   * @param border
+   * @param radius
+   * @see com.google.maps.android.SphericalUtil
+   */
+  public static double area(LocationList border) {
+    if (border.size() < 3) {
+      return 0;
+    }
+    double total = 0;
+    Location prev = border.last();
+    double prevTanLat = tan((PI / 2.0 - prev.latRad) / 2.0);
+    double prevLon = prev.lonRad;
+    // For each edge, accumulate the signed area of the triangle
+    // formed by the North Pole and that edge ("polar triangle").
+    for (Location loc : border) {
+      double tanLat = tan((PI / 2.0 - loc.latRad) / 2.0);
+      double lon = loc.lonRad;
+      total += polarTriangleArea(tanLat, lon, prevTanLat, prevLon);
+      prevTanLat = tanLat;
+      prevLon = lon;
+    }
+    return total * EARTH_RADIUS_MEAN * EARTH_RADIUS_MEAN;
+  }
+
+  /*
+   * Returns the signed area of a triangle which has North Pole as a vertex.
+   * Formula derived from "Area of a spherical triangle given two edges and the
+   * included angle" as per "Spherical Trigonometry" by Todhunter, page 71,
+   * section 103, point 2.
+   *
+   * See http://books.google.com/books?id=3uBHAAAAIAAJ&pg=PA71
+   *
+   * The arguments named "tan" are tan((pi/2 - latitude)/2).
+   */
+  private static double polarTriangleArea(
+      double tan1, double lon1,
+      double tan2, double lon2) {
+    double dLon = lon1 - lon2;
+    double t = tan1 * tan2;
+    return 2.0 * atan2(t * sin(dLon), 1.0 + t * cos(dLon));
+  }
+
   /**
    * Calculates the angle between two {@code Location}s using the <a
    * href="http://en.wikipedia.org/wiki/Haversine_formula" target="_blank">

src/test/java/gov/usgs/earthquake/nshmp/geo/LocationsTests.java

@@ -14,6 +14,8 @@ import org.junit.jupiter.api.Test;
 
 import com.google.common.collect.Lists;
 
+import gov.usgs.earthquake.nshmp.Maths;
+
 class LocationsTests {
 
   // private static LocationList locs1, locs2;
@@ -526,6 +528,33 @@ class LocationsTests {
     });
   }
 
+  // closed border
+  static final LocationList BORDER_1 = LocationList.builder()
+      .add(-119.41882, 34.99054)
+      .add(-118.96453, 35.20572)
+      .add(-119.41640, 35.52851)
+      .add(-119.90303, 35.12355)
+      .add(-119.41882, 34.99054)
+      .build();
+
+  // open border
+  static final LocationList BORDER_2 = LocationList.copyOf(BORDER_1.subList(0, 4));
+
+  // not a border
+  static final LocationList BORDER_3 = LocationList.copyOf(BORDER_1.subList(0, 1));
+
+  static double expectedArea = 2550.6996;
+
+  @Test
+  void testArea() {
+    double actualArea1 = Maths.round(Locations.area(BORDER_1), 4);
+    assertEquals(expectedArea, actualArea1);
+    double actualArea2 = Maths.round(Locations.area(BORDER_2), 4);
+    assertEquals(expectedArea, actualArea2);
+    double actualArea3 = Locations.area(BORDER_3);
+    assertEquals(0.0, actualArea3);
+  }
+
   /**
    * DEVELOPER NOTE
    *