diff --git a/src/main/java/gov/usgs/earthquake/nshmp/geo/Locations.java b/src/main/java/gov/usgs/earthquake/nshmp/geo/Locations.java index 2ec9ce204ca1b742c9c9e4a58fb3f229172a83ea..d9b0a3cc8788d6930bb43053723903a6ca0f80b3 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/geo/Locations.java +++ b/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"> diff --git a/src/test/java/gov/usgs/earthquake/nshmp/geo/LocationsTests.java b/src/test/java/gov/usgs/earthquake/nshmp/geo/LocationsTests.java index 5e49dd20375d358177f9351a717377878e3a5579..53f7e084a0130ccb92d710ade1d8d469a83f2ae0 100644 --- a/src/test/java/gov/usgs/earthquake/nshmp/geo/LocationsTests.java +++ b/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 *