diff --git a/CHANGELOG.md b/CHANGELOG.md index 946feb16..e9731ac4 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/main/java/edu/ie3/util/geo/GeoUtils.java b/src/main/java/edu/ie3/util/geo/GeoUtils.java index d276aada..2735cc40 100644 --- a/src/main/java/edu/ie3/util/geo/GeoUtils.java +++ b/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 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 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 calcOrderedCoordinateDistances( */ public static ComparableQuantity calcHaversine( double latA, double lngA, double latB, double lngB) { - - ComparableQuantity dLat = Quantities.getQuantity(toRadians(latB - latA), RADIAN); - ComparableQuantity 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); } /** diff --git a/src/main/scala/edu/ie3/util/geo/RichGeometries.scala b/src/main/scala/edu/ie3/util/geo/RichGeometries.scala index ddf1c678..86dd8094 100644 --- a/src/main/scala/edu/ie3/util/geo/RichGeometries.scala +++ b/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