Fixing runtime issues of GeoUtils.calcHaversine

CHANGELOG.md

@@ -6,6 +6,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased/Snapshot]
 
+### Fixed
+- Fix catastrophic runtime of `GeoUtils.calcHaversine` due to Quantities [#508](https://github.com/ie3-institute/PowerSystemUtils/issues/508)
+
 ## [2.2.1]
 
 ### Changed

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

@@ -5,15 +5,13 @@
 */
 package edu.ie3.util.geo;
 
-import static edu.ie3.util.quantities.PowerSystemUnits.METRE;
-import static edu.ie3.util.quantities.PowerSystemUnits.RADIAN;
+import static edu.ie3.util.quantities.PowerSystemUnits.*;
 import static java.lang.Math.*;
 
 import edu.ie3.util.exceptions.GeoException;
 import java.util.*;
 import java.util.stream.IntStream;
 import javax.measure.Quantity;
-import javax.measure.quantity.Angle;
 import javax.measure.quantity.Length;
 import org.locationtech.jts.algorithm.ConvexHull;
 import org.locationtech.jts.geom.*;
@@ -26,8 +24,10 @@ public class GeoUtils {
   public static final GeometryFactory DEFAULT_GEOMETRY_FACTORY =
       new GeometryFactory(new PrecisionModel(), 4326);
 
+  public static final double EARTH_RADIUS_M = 6378137.0;
+
   public static final ComparableQuantity<Length> EARTH_RADIUS =
-      Quantities.getQuantity(6378137.0, METRE);
+      Quantities.getQuantity(EARTH_RADIUS_M, METRE);
 
   protected GeoUtils() {
     throw new IllegalStateException("Utility classes cannot be instantiated.");
@@ -125,6 +125,39 @@ public static List<CoordinateDistance> calcOrderedCoordinateDistances(
         .toList();
   }
 
+  /**
+   * Calculates between two coordinates on earth's surface (great circle distance). Does not use
+   * Quantities for faster calculation.
+   *
+   * @param latA latitude of coordinate a
+   * @param lngA longitude of coordinate a
+   * @param latB latitude of coordinate b
+   * @param lngB longitude of coordinate b
+   * @return The distance between both coordinates in metres
+   */
+  public static double calcHaversineMetres(double latA, double lngA, double latB, double lngB) {
+    double dLatRad = toRadians(latB - latA);
+    double dLonRad = toRadians(lngB - lngA);
+    double a =
+        sin(dLatRad / 2) * sin(dLatRad / 2)
+            + cos(toRadians(latA)) * cos(toRadians(latB)) * sin(dLonRad / 2) * sin(dLonRad / 2);
+    double c = 2 * atan2(sqrt(a), sqrt(1 - a));
+    return EARTH_RADIUS_M * c;
+  }
+
+  /**
+   * Calculates the distance between two coordinates on earth's surface (great circle distance).
+   * Does not use Quantities for faster calculation.
+   *
+   * @param coordinateA coordinate a
+   * @param coordinateB coordinate b
+   * @return the distance between the coordinates in metres
+   */
+  public static double calcHaversineMetres(Coordinate coordinateA, Coordinate coordinateB) {
+    return calcHaversineMetres(
+        coordinateA.getY(), coordinateA.getX(), coordinateB.getY(), coordinateB.getX());
+  }
+
   /**
    * Calculates between two coordinates on earth's surface (great circle distance).
    *
@@ -136,17 +169,7 @@ public static List<CoordinateDistance> calcOrderedCoordinateDistances(
    */
   public static ComparableQuantity<Length> calcHaversine(
       double latA, double lngA, double latB, double lngB) {
-
-    ComparableQuantity<Angle> dLat = Quantities.getQuantity(toRadians(latB - latA), RADIAN);
-    ComparableQuantity<Angle> dLon = Quantities.getQuantity(toRadians(lngB - lngA), RADIAN);
-    double a =
-        sin(dLat.getValue().doubleValue() / 2) * sin(dLat.getValue().doubleValue() / 2)
-            + cos(toRadians(latA))
-                * cos(toRadians(latB))
-                * sin(dLon.getValue().doubleValue() / 2)
-                * sin(dLon.getValue().doubleValue() / 2);
-    double c = 2 * atan2(sqrt(a), sqrt(1 - a));
-    return EARTH_RADIUS.multiply(c);
+    return Quantities.getQuantity(calcHaversineMetres(latA, lngA, latB, lngB), METRE);
   }
 
   /**

src/main/scala/edu/ie3/util/geo/RichGeometries.scala

@@ -6,7 +6,6 @@
 package edu.ie3.util.geo
 
 import edu.ie3.util.exceptions.GeoException
-import edu.ie3.util.geo.GeoUtils
 import edu.ie3.util.quantities.QuantityUtils.RichQuantityDouble
 import org.locationtech.jts.geom.{
   Coordinate,
@@ -19,7 +18,7 @@ import tech.units.indriya.ComparableQuantity
 
 import javax.measure.quantity.{Area, Length}
 import scala.math.abs
-import scala.util.{Failure, Success, Try}
+import scala.util.{Failure, Try}
 
 object RichGeometries {
 
@@ -37,6 +36,19 @@ object RichGeometries {
     ): ComparableQuantity[Length] =
       GeoUtils.calcHaversine(coordinate, coordinateB)
 
+    /** Calculates the great circle distance between two coordinates. Does not
+      * use Quantities for faster calculation.
+      *
+      * @param coordinateB
+      *   coordinate b
+      * @return
+      *   the distance between the coordinates in metres
+      */
+    def haversineDistanceMetres(
+        coordinateB: Coordinate
+    ): Double =
+      GeoUtils.calcHaversineMetres(coordinate, coordinateB)
+
     /** Checks if the coordinate lies between two coordinates a and b by
       * comparing the distances between a and b with the sum of distances
       * between the coordinate and a and the coordinate and b
@@ -48,20 +60,17 @@ object RichGeometries {
       * @param epsilon
       *   permitted relative deviation
       * @return
-      *   whether or not the coordinate lies between
+      *   whether the coordinate lies between
       */
     def isBetween(
         a: Coordinate,
         b: Coordinate,
         epsilon: Double = 1e-12
     ): Boolean = {
-      val distance = a.haversineDistance(b)
-      val distancePassingMe = a
-        .haversineDistance(coordinate)
-        .add(coordinate.haversineDistance(b))
-        .getValue
-        .doubleValue
-      abs(1 - (distancePassingMe / distance.getValue.doubleValue())) < epsilon
+      val distance = a.haversineDistanceMetres(b)
+      val distancePassingMe = a.haversineDistanceMetres(coordinate)
+        + coordinate.haversineDistanceMetres(b)
+      abs(1 - (distancePassingMe / distance)) < epsilon
     }
 
     /** Creates a [[Point]] from this coordinate
@@ -127,7 +136,7 @@ object RichGeometries {
     }
 
     /** Checks whether the polygon contains the coordinate. Uses "covers()"
-      * insted of "contains()" so borders are included.
+      * instead of "contains()" so borders are included.
       *
       * @param coordinate
       *   the coordinate to check