Merge pull request #321 from ie3-institute/ms/#320-adding-new-feature-to-GeoUtils

Adding new feature to geo utils

Author: t-ober

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `OsmEntity` and `OsmContainer` to provide a simple, lightweight representation of openstreetmap data
 - QuantityUtils previously implemented in SIMONA [#288](https://github.com/ie3-institute/PowerSystemUtils/issues/288)
 - Enhanced `RichQuantityDouble` with new units and generic method [#312](https://github.com/ie3-institute/PowerSystemUtils/issues/312), [#322](https://github.com/ie3-institute/PowerSystemUtils/issues/322)
+- Calculation of bounding box to `GeoUtils` [#320](https://github.com/ie3-institute/PowerSystemUtils/issues/320)
 
 ### Changed
 - Refactored `GeoUtils`, moved them to the scala package and tailored them toward the `loactiontec.jts` Geometries used in the `OsmContainer` [#163](https://github.com/ie3-institute/PowerSystemUtils/issues/163)

src/main/java/edu/ie3/util/geo/GeoUtils.java

@@ -5,7 +5,8 @@
 */
 package edu.ie3.util.geo;
 
-import static edu.ie3.util.quantities.PowerSystemUnits.*;
+import static edu.ie3.util.quantities.PowerSystemUnits.METRE;
+import static edu.ie3.util.quantities.PowerSystemUnits.RADIAN;
 import static java.lang.Math.*;
 
 import edu.ie3.util.exceptions.GeoException;
@@ -16,15 +17,7 @@
 import javax.measure.quantity.Angle;
 import javax.measure.quantity.Length;
 import org.locationtech.jts.algorithm.ConvexHull;
-import org.locationtech.jts.geom.Coordinate;
-import org.locationtech.jts.geom.Geometry;
-import org.locationtech.jts.geom.GeometryCollection;
-import org.locationtech.jts.geom.GeometryFactory;
-import org.locationtech.jts.geom.LineString;
-import org.locationtech.jts.geom.LinearRing;
-import org.locationtech.jts.geom.Point;
-import org.locationtech.jts.geom.Polygon;
-import org.locationtech.jts.geom.PrecisionModel;
+import org.locationtech.jts.geom.*;
 import org.locationtech.jts.geom.impl.CoordinateArraySequence;
 import org.locationtech.jts.math.Vector2D;
 import tech.units.indriya.ComparableQuantity;
@@ -34,7 +27,7 @@ public class GeoUtils {
   public static final GeometryFactory DEFAULT_GEOMETRY_FACTORY =
       new GeometryFactory(new PrecisionModel(), 4326);
 
-  private static final ComparableQuantity<Length> EARTH_RADIUS =
+  public static final ComparableQuantity<Length> EARTH_RADIUS =
       Quantities.getQuantity(6378137.0, METRE);
 
   protected GeoUtils() {
@@ -211,6 +204,45 @@ public static ComparableQuantity<Length> calcHaversine(LineString lineString) {
     return y;
   }
 
+  /**
+   * Method for calculating a bounding box around a point. This method uses an inverse haversine to
+   * turn a distance into x- and y-deltas. These deltas are used to create an envelope which
+   * represents the bounding box. The methode can be found here: <a
+   * href="https://math.stackexchange.com/questions/474602/reverse-use-of-haversine-formula">source</a>
+   *
+   * @param coordinate the starting point for the calculation
+   * @param distance maximal distance in x- and y-direction
+   * @return envelope
+   */
+  public static Envelope calculateBoundingBox(
+      Point coordinate, ComparableQuantity<Length> distance) {
+    // y-degrees are evenly spaced, so we can just divide a distance
+    // by the earth's radius to get a y-delta in radians
+    double deltaY = distance.divide(EARTH_RADIUS).getValue().doubleValue();
+
+    // because the spacing between x-degrees change between the equator
+    // and the poles, we need to calculate the x-delta using the inverse
+    // haversine formula
+    double sinus = sin(deltaY / 2);
+    double squaredSinus = pow(sinus, 2);
+    double cosine = cos(toRadians(coordinate.getY()));
+    double squaredCosine = pow(cosine, 2);
+
+    double deltaX = 2 * asin(sqrt(squaredSinus / squaredCosine));
+
+    // turning the deltas to degree
+    double deltaXDegree = toDegrees(deltaX);
+    double deltaYDegree = toDegrees(deltaY);
+
+    // calculating minimums and maximums for longitude and latitude and returning them as an
+    // envelope
+    return new Envelope(
+        coordinate.getX() - deltaXDegree,
+        coordinate.getX() + deltaXDegree,
+        coordinate.getY() - deltaYDegree,
+        coordinate.getY() + deltaYDegree);
+  }
+
   /**
    * Builds a convex hull from a set of latitude/longitude coordinates.
    *

src/test/groovy/edu/ie3/util/geo/GeoUtilsTest.groovy

@@ -5,14 +5,8 @@
  */
 package edu.ie3.util.geo
 
-import static edu.ie3.util.quantities.PowerSystemUnits.METRE
-
-import edu.ie3.util.quantities.PowerSystemUnits
 import edu.ie3.util.quantities.QuantityUtil
-import org.locationtech.jts.geom.Coordinate
-import org.locationtech.jts.geom.GeometryFactory
-import org.locationtech.jts.geom.LineString
-import org.locationtech.jts.geom.PrecisionModel
+import org.locationtech.jts.geom.*
 import org.locationtech.jts.io.geojson.GeoJsonReader
 import spock.lang.Shared
 import spock.lang.Specification
@@ -21,6 +15,8 @@ import tech.units.indriya.quantity.Quantities
 
 import javax.measure.quantity.Length
 
+import static edu.ie3.util.quantities.PowerSystemUnits.METRE
+
 class GeoUtilsTest extends Specification {
 
     @Shared
@@ -294,4 +290,42 @@ class GeoUtilsTest extends Specification {
             Math.abs(actual.x - 7.5) < 1e-9
             Math.abs(actual.y - 50d) < 1e-9
     }
+
+    def "GeoUtils calculates x-delta (longitude) correctly"() {
+        given:
+        Point coordinate = GeoUtils.buildPoint(50, 7)
+        ComparableQuantity<Length> distanceX = GeoUtils.calcHaversine(50, 7, 50, 6)
+
+        when:
+        Envelope envelope = GeoUtils.calculateBoundingBox(coordinate, distanceX)
+
+        then:
+        envelope.minX == 6
+        envelope.maxX == 8
+
+        double yMin = 49.35721717797476
+        double yMax = 50.64278282202524
+
+        envelope.minY == yMin
+        envelope.maxY == yMax
+    }
+
+    def "GeoUtils calculates y-delta (latitude) correctly"() {
+        given:
+        Point coordinate = GeoUtils.buildPoint(50, 7)
+        ComparableQuantity<Length> distanceX = GeoUtils.calcHaversine(51, 7, 52, 7)
+
+        when:
+        Envelope envelope = GeoUtils.calculateBoundingBox(coordinate, distanceX)
+
+        then:
+        double xMin = 5.444248126340358
+        double xMax = 8.555751873659641
+
+        envelope.minX == xMin
+        envelope.maxX == xMax
+
+        envelope.minY == 49
+        envelope.maxY == 51
+    }
 }