Implement DistanceCalculator with multiple methods and refactor coordinate handling

- Introduced `DistanceCalculator` for 2D and 3D distance computation between coordinates and map grids.
- Added support for various calculation methods: Haversine, spherical law of cosines, and extended UTM Euclidean.
- Implemented the `Distance`, `DistanceCalculationMethod`, and `DistanceDimension` classes for structured result handling and flexible computation options.
- Extended `CoordinateConverter` for UTM projection enhancements, including forced zone and latitude band support.
- Updated `BalloonLiveParser` and `GridCoordinate` to accommodate new distance models and UTM latitude band handling.
This commit is contained in:
Jan Meinl
2026-05-16 04:33:47 +02:00
parent 384eb5729e
commit fa888405bd
7 changed files with 422 additions and 13 deletions
@@ -0,0 +1,47 @@
package dev.coph.flightscore.backend.coordinate.distance;
import lombok.Getter;
import lombok.experimental.Accessors;
/**
* Result of a {@link DistanceCalculator} call. Components that were not
* requested via {@link DistanceDimension} (or that could not be computed,
* e.g. {@link #vertical()} when an altitude was missing) are {@code null}.
*/
@Getter
@Accessors(fluent = true)
public class Distance {
/** 2D horizontal distance in metres, or {@code null} if not requested. */
private final Double horizontal;
/** Vertical separation in metres, or {@code null} if not computed. */
private final Double vertical;
/** 3D slant distance in metres, or {@code null} if not requested / computable. */
private final Double slant3d;
public Distance(Double horizontal, Double vertical, Double slant3d) {
this.horizontal = horizontal;
this.vertical = vertical;
this.slant3d = slant3d;
}
@Override
public String toString() {
StringBuilder sb = new StringBuilder("Distance{");
boolean first = true;
if (horizontal != null) {
sb.append("horizontal=").append(String.format("%.2fm", horizontal));
first = false;
}
if (vertical != null) {
if (!first) sb.append(", ");
sb.append("vertical=").append(String.format("%.2fm", vertical));
first = false;
}
if (slant3d != null) {
if (!first) sb.append(", ");
sb.append("slant3d=").append(String.format("%.2fm", slant3d));
}
return sb.append('}').toString();
}
}
@@ -0,0 +1,27 @@
package dev.coph.flightscore.backend.coordinate.distance;
/**
* Algorithms supported by {@link DistanceCalculator}.
*/
public enum DistanceCalculationMethod {
/** Great-circle distance via the haversine formula (spherical earth). */
HAVERSINE("Haversine formula (spherical great circle)"),
/** Great-circle distance via the spherical law of cosines. */
SPHERICAL_LAW_OF_COSINES("Spherical law of cosines (great circle)"),
/**
* Planar (Euclidean) distance on the UTM grid. Both points are projected
* with the same central meridian (the source point's UTM zone) so distances
* remain continuous across zone boundaries.
*/
UTM_EUCLIDEAN("Extended UTM Euclidean (single-zone planar)");
private final String description;
DistanceCalculationMethod(String description) {
this.description = description;
}
public String description() {
return description;
}
}
@@ -0,0 +1,239 @@
package dev.coph.flightscore.backend.coordinate.distance;
import dev.coph.flightscore.backend.coordinate.Altitude;
import dev.coph.flightscore.backend.coordinate.Coordinate;
import dev.coph.flightscore.backend.map.CoordinateConverter;
import dev.coph.flightscore.backend.map.GridCoordinate;
import dev.coph.flightscore.backend.map.MapDatum;
import dev.coph.flightscore.backend.map.MapGrid;
/**
* Computes 2D and / or 3D distances between any pair of {@link Coordinate} and
* / or {@link GridCoordinate}.
* <p>
* Three formulas are available via {@link DistanceCalculationMethod}:
* <ul>
* <li>{@link DistanceCalculationMethod#HAVERSINE} the standard great-circle
* formula on a spherical earth.</li>
* <li>{@link DistanceCalculationMethod#SPHERICAL_LAW_OF_COSINES} the
* trigonometric great-circle formula.</li>
* <li>{@link DistanceCalculationMethod#UTM_EUCLIDEAN} planar distance on an
* <i>extended</i> UTM grid: the second point is projected with the first
* point's central meridian, which eliminates the discontinuity at zone
* boundaries (the wedge-shaped overlap).</li>
* </ul>
* Conversions between coordinate systems and datums are delegated to
* {@link CoordinateConverter}, so the calculator works for every datum and
* grid known to {@link MapDatum} / {@link MapGrid}.
* <h2>Choosing what to compute</h2>
* Pass a {@link DistanceDimension} to skip work you don't need:
* <ul>
* <li>{@link DistanceDimension#TWO_D} only the horizontal distance is
* computed (no altitude lookup, no slant).</li>
* <li>{@link DistanceDimension#THREE_D} only the 3D slant distance is
* exposed; horizontal is computed internally but not returned.</li>
* <li>{@link DistanceDimension#BOTH} both are computed and returned. This
* is the default for overloads that don't take a dimension.</li>
* </ul>
* Vertical separation is only computable when both inputs carry an altitude
* (GPS preferred, falling back to barometric); otherwise the slant component
* is {@code null} even if it was requested.
*/
public final class DistanceCalculator {
/** Mean radius of the WGS84 ellipsoid (IUGG mean radius R1). */
private static final double EARTH_RADIUS_METERS = 6_371_008.8;
private static final double DEG_TO_RAD = Math.PI / 180.0;
private DistanceCalculator() {
}
// ------------------------------------------------------------------
// Coordinate × Coordinate
// ------------------------------------------------------------------
public static Distance distance(Coordinate from, Coordinate to,
DistanceCalculationMethod method) {
return distance(from, to, method, DistanceDimension.BOTH, MapDatum.WGS84);
}
public static Distance distance(Coordinate from, Coordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension) {
return distance(from, to, method, dimension, MapDatum.WGS84);
}
/**
* @param datum datum the input coordinates are referenced to. Used by
* {@link DistanceCalculationMethod#UTM_EUCLIDEAN} for the
* projection; ignored by the spherical methods, which always
* assume WGS84 lat/lon.
*/
public static Distance distance(Coordinate from, Coordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension,
MapDatum datum) {
Double horizontal = dimension.needsHorizontal() || dimension.needsSlant()
? horizontalDistance(from, to, method, datum)
: null;
return assemble(from, to, horizontal, dimension);
}
// ------------------------------------------------------------------
// GridCoordinate × GridCoordinate
// ------------------------------------------------------------------
public static Distance distance(GridCoordinate from, GridCoordinate to,
DistanceCalculationMethod method) {
return distance(from, to, method, DistanceDimension.BOTH, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(GridCoordinate from, GridCoordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension) {
return distance(from, to, method, dimension, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(GridCoordinate from, GridCoordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension,
MapDatum datum, MapGrid grid) {
return distance(from, datum, grid, to, datum, grid, method, dimension);
}
/**
* Distance between two grid coordinates that may live in different
* {@code (datum, grid)} systems. Both are converted back to geographic
* coordinates with {@link CoordinateConverter#toGeographic} before the
* chosen method is applied.
*/
public static Distance distance(GridCoordinate from, MapDatum fromDatum, MapGrid fromGrid,
GridCoordinate to, MapDatum toDatum, MapGrid toGrid,
DistanceCalculationMethod method,
DistanceDimension dimension) {
Coordinate fromCoord = CoordinateConverter.toGeographic(from, fromDatum, fromGrid);
Coordinate toCoord = CoordinateConverter.toGeographic(to, toDatum, toGrid);
return distance(fromCoord, toCoord, method, dimension, fromDatum);
}
// ------------------------------------------------------------------
// Mixed Coordinate × GridCoordinate
// ------------------------------------------------------------------
public static Distance distance(Coordinate from, GridCoordinate to,
DistanceCalculationMethod method) {
return distance(from, to, method, DistanceDimension.BOTH, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(Coordinate from, GridCoordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension) {
return distance(from, to, method, dimension, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(Coordinate from, GridCoordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension,
MapDatum datum, MapGrid grid) {
return distance(from, CoordinateConverter.toGeographic(to, datum, grid), method, dimension, datum);
}
public static Distance distance(GridCoordinate from, Coordinate to,
DistanceCalculationMethod method) {
return distance(from, to, method, DistanceDimension.BOTH, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(GridCoordinate from, Coordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension) {
return distance(from, to, method, dimension, MapDatum.WGS84, MapGrid.LAT_LONG);
}
public static Distance distance(GridCoordinate from, Coordinate to,
DistanceCalculationMethod method,
DistanceDimension dimension,
MapDatum datum, MapGrid grid) {
return distance(CoordinateConverter.toGeographic(from, datum, grid), to, method, dimension, datum);
}
// ------------------------------------------------------------------
// Horizontal computation (single method dispatch)
// ------------------------------------------------------------------
private static double horizontalDistance(Coordinate from, Coordinate to,
DistanceCalculationMethod method, MapDatum datum) {
return switch (method) {
case HAVERSINE -> haversine(from, to);
case SPHERICAL_LAW_OF_COSINES -> sphericalLawOfCosines(from, to);
case UTM_EUCLIDEAN -> utmEuclidean(from, to, datum);
};
}
private static double haversine(Coordinate from, Coordinate to) {
double lat1 = from.latitude() * DEG_TO_RAD;
double lat2 = to.latitude() * DEG_TO_RAD;
double dLat = lat2 - lat1;
double dLon = (to.longitude() - from.longitude()) * DEG_TO_RAD;
double sinHalfLat = Math.sin(dLat / 2.0);
double sinHalfLon = Math.sin(dLon / 2.0);
double a = sinHalfLat * sinHalfLat
+ Math.cos(lat1) * Math.cos(lat2) * sinHalfLon * sinHalfLon;
return EARTH_RADIUS_METERS * 2.0 * Math.asin(Math.min(1.0, Math.sqrt(a)));
}
private static double sphericalLawOfCosines(Coordinate from, Coordinate to) {
double lat1 = from.latitude() * DEG_TO_RAD;
double lat2 = to.latitude() * DEG_TO_RAD;
double dLon = (to.longitude() - from.longitude()) * DEG_TO_RAD;
double cosCentral = Math.sin(lat1) * Math.sin(lat2)
+ Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon);
cosCentral = Math.max(-1.0, Math.min(1.0, cosCentral));
return EARTH_RADIUS_METERS * Math.acos(cosCentral);
}
private static double utmEuclidean(Coordinate from, Coordinate to, MapDatum datum) {
// Project the source naturally; this fixes the zone (and hemisphere) of
// the planar frame in which the distance is measured.
GridCoordinate fromUtm = CoordinateConverter.convert(from, datum, MapDatum.WGS84, MapGrid.UTM);
// Project the target into the same zone / hemisphere as the source, so
// both points share one central meridian and the wedge-shaped overlap
// between zones causes no jump in the planar coordinates.
Coordinate toWgs84 = CoordinateConverter.transformDatum(to, datum, MapDatum.WGS84);
GridCoordinate toUtm = fromUtm.utmZone() != null
? CoordinateConverter.toUtmGrid(toWgs84, MapDatum.WGS84,
fromUtm.utmZone(), fromUtm.southernHemisphere())
: CoordinateConverter.toGrid(toWgs84, MapDatum.WGS84, MapGrid.UTM);
return Math.hypot(toUtm.easting() - fromUtm.easting(),
toUtm.northing() - fromUtm.northing());
}
// ------------------------------------------------------------------
// Result assembly
// ------------------------------------------------------------------
private static Distance assemble(Coordinate from, Coordinate to,
Double horizontal, DistanceDimension dimension) {
Double vertical = null;
Double slant = null;
if (dimension.needsSlant()) {
Altitude fromAltitude = selectAltitude(from);
Altitude toAltitude = selectAltitude(to);
if (fromAltitude != null && toAltitude != null && horizontal != null) {
vertical = Math.abs(toAltitude.meters() - fromAltitude.meters());
slant = Math.hypot(horizontal, vertical);
}
}
Double exposedHorizontal = dimension.needsHorizontal() ? horizontal : null;
return new Distance(exposedHorizontal, vertical, slant);
}
private static Altitude selectAltitude(Coordinate coord) {
return coord.gpsAltitude() != null ? coord.gpsAltitude() : coord.barometricAltitude();
}
}
@@ -0,0 +1,27 @@
package dev.coph.flightscore.backend.coordinate.distance;
/**
* Selects which distance components {@link DistanceCalculator} actually
* computes. Components that are not requested are skipped entirely; the
* corresponding accessors on {@link Distance} return {@code null}.
*/
public enum DistanceDimension {
/** Only the 2D horizontal (great-circle / planar) distance. */
TWO_D,
/**
* Only the 3D slant distance through space. Internally still requires the
* horizontal component to combine with the vertical separation, but the
* horizontal value itself is not exposed on the returned {@link Distance}.
*/
THREE_D,
/** Both the horizontal and the 3D slant distance. */
BOTH;
boolean needsHorizontal() {
return this == TWO_D || this == BOTH;
}
boolean needsSlant() {
return this == THREE_D || this == BOTH;
}
}
@@ -80,9 +80,9 @@ public final class CoordinateConverter {
switch (grid.projectionMethod()) {
case NONE:
return new GridCoordinate(coordinate.longitude(), coordinate.latitude(),
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, false);
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, null, false);
case TRANSVERSE_MERCATOR:
return transverseMercatorForward(coordinate, datum, grid);
return transverseMercatorForward(coordinate, datum, grid, null, false);
case SWISS_OBLIQUE_MERCATOR:
return swissObliqueMercatorForward(coordinate, datum, grid);
default:
@@ -90,6 +90,25 @@ public final class CoordinateConverter {
}
}
/**
* Projects a {@link Coordinate} (referenced to {@code datum}) onto UTM,
* forcing the projection into the given {@code zone} and hemisphere instead
* of using the natural zone for the coordinate's longitude.
* <p>
* This is the "extended UTM" projection used to compute Euclidean distances
* between points that straddle a zone boundary: both points are projected
* with the same central meridian, which keeps the planar geometry consistent
* across the wedge-shaped overlap region. The returned coordinate carries
* the forced {@code zone} and hemisphere; the latitude band letter is still
* derived from the coordinate's actual latitude.
*/
public static GridCoordinate toUtmGrid(Coordinate coordinate, MapDatum datum, int zone, boolean southern) {
if (zone < 1 || zone > 60) {
throw new IllegalArgumentException("UTM zone must be in [1, 60], was " + zone);
}
return transverseMercatorForward(coordinate, datum, MapGrid.UTM, zone, southern);
}
/**
* Inverse projection of a grid coordinate (projected with {@code grid},
* referenced to {@code datum}) onto a geographic {@link Coordinate}. Both
@@ -268,7 +287,8 @@ public final class CoordinateConverter {
// Transverse Mercator (Krueger series, Karney's algorithm)
// ------------------------------------------------------------------
private static GridCoordinate transverseMercatorForward(Coordinate coordinate, MapDatum datum, MapGrid grid) {
private static GridCoordinate transverseMercatorForward(Coordinate coordinate, MapDatum datum, MapGrid grid,
Integer forcedZone, boolean forcedSouthern) {
double a = datum.semiMajorAxis();
double f = datum.flattening();
double e = datum.firstEccentricity();
@@ -281,8 +301,13 @@ public final class CoordinateConverter {
double centralMeridian;
double falseNorthing = grid.falseNorthing();
if (grid.isZoneDependent()) {
zone = utmZone(coordinate.longitude());
southern = coordinate.latitude() < 0.0;
if (forcedZone != null) {
zone = forcedZone;
southern = forcedSouthern;
} else {
zone = utmZone(coordinate.longitude());
southern = coordinate.latitude() < 0.0;
}
centralMeridian = Math.toRadians(zone * 6.0 - 183.0);
falseNorthing = southern ? 10_000_000.0 : 0.0;
} else {
@@ -299,8 +324,9 @@ public final class CoordinateConverter {
double[] series = tmForwardSeries(phi, deltaLambda, e, alpha);
double easting = falseEasting + k0 * rectifyingRadius * series[1];
double northing = falseNorthing + k0 * rectifyingRadius * (series[0] - xiOrigin);
Character letter = zone != null ? utmLatitudeBandLetter(coordinate.latitude()) : null;
return new GridCoordinate(easting, northing,
coordinate.barometricAltitude(), coordinate.gpsAltitude(), zone, southern);
coordinate.barometricAltitude(), coordinate.gpsAltitude(), zone, letter, southern);
}
private static Coordinate transverseMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
@@ -424,6 +450,32 @@ public final class CoordinateConverter {
return Math.max(1, Math.min(60, zone));
}
private static Character utmLatitudeBandLetter(double latitude) {
if (latitude < -80.0 || latitude > 84.0) {
return null;
}
if (latitude < -72.0) return 'C';
if (latitude < -64.0) return 'D';
if (latitude < -56.0) return 'E';
if (latitude < -48.0) return 'F';
if (latitude < -40.0) return 'G';
if (latitude < -32.0) return 'H';
if (latitude < -24.0) return 'J';
if (latitude < -16.0) return 'K';
if (latitude < -8.0) return 'L';
if (latitude < 0.0) return 'M';
if (latitude < 8.0) return 'N';
if (latitude < 16.0) return 'P';
if (latitude < 24.0) return 'Q';
if (latitude < 32.0) return 'R';
if (latitude < 40.0) return 'S';
if (latitude < 48.0) return 'T';
if (latitude < 56.0) return 'U';
if (latitude < 64.0) return 'V';
if (latitude < 72.0) return 'W';
return 'X';
}
// ------------------------------------------------------------------
// Swiss oblique Mercator (rigorous swisstopo algorithm)
// ------------------------------------------------------------------
@@ -449,7 +501,7 @@ public final class CoordinateConverter {
double easting = c.r * pseudoLongitude + grid.falseEasting();
double northing = c.r * atanh(sinPseudoLatitude) + grid.falseNorthing();
return new GridCoordinate(easting, northing,
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, false);
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, null, false);
}
private static Coordinate swissObliqueMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
@@ -30,31 +30,48 @@ public class GridCoordinate {
private final Altitude gpsAltitude;
/** UTM zone number (1..60) for zone dependent grids, otherwise {@code null}. */
private final Integer utmZone;
/** UTM latitude band letter (C-X, excluding I and O) for zone dependent grids, otherwise {@code null}. */
private final Character utmLatitudeBandLetter;
/** Whether the coordinate lies on the southern hemisphere (relevant for UTM false northing). */
private final boolean southernHemisphere;
public GridCoordinate(double easting, double northing,
Altitude barometricAltitude, Altitude gpsAltitude,
Integer utmZone, boolean southernHemisphere) {
Integer utmZone, Character utmLatitudeBandLetter, boolean southernHemisphere) {
this.easting = easting;
this.northing = northing;
this.barometricAltitude = barometricAltitude;
this.gpsAltitude = gpsAltitude;
this.utmZone = utmZone;
this.utmLatitudeBandLetter = utmLatitudeBandLetter;
this.southernHemisphere = southernHemisphere;
}
public GridCoordinate(double easting, double northing,
Altitude barometricAltitude, Altitude gpsAltitude) {
this(easting, northing, barometricAltitude, gpsAltitude, null, false);
this(easting, northing, barometricAltitude, gpsAltitude, null, null, false);
}
public GridCoordinate(double easting, double northing, Altitude altitude) {
this(easting, northing, altitude, altitude, null, false);
this(easting, northing, altitude, altitude, null, null, false);
}
public GridCoordinate(double easting, double northing) {
this(easting, northing, null, null, null, false);
this(easting, northing, null, null, null, null, false);
}
/**
* Combined UTM zone designator ({@code "32U"}, {@code "17S"}, ...) for
* zone dependent grids, or {@code null} otherwise. If the zone is known
* but the latitude band letter is not, only the number is returned.
*/
public String utmZoneDesignator() {
if (utmZone == null) {
return null;
}
return utmLatitudeBandLetter != null
? utmZone + Character.toString(utmLatitudeBandLetter)
: Integer.toString(utmZone);
}
@Override
@@ -62,7 +79,7 @@ public class GridCoordinate {
return "GridCoordinate{easting=" + easting + ", northing=" + northing
+ ", barometricAltitude=" + barometricAltitude
+ ", gpsAltitude=" + gpsAltitude
+ (utmZone != null ? ", utmZone=" + utmZone + ", southernHemisphere=" + southernHemisphere : "")
+ (utmZone != null ? ", utmZone=" + utmZoneDesignator() + ", southernHemisphere=" + southernHemisphere : "")
+ '}';
}
}
@@ -502,7 +502,7 @@ public class BalloonLiveParser implements TrackParser {
Altitude baro = declaredAltitude != null ? declaredAltitude : reference.barometricAltitude();
Altitude gps = declaredAltitude != null ? declaredAltitude : reference.gpsAltitude();
GridCoordinate decl = new GridCoordinate(eVal, nVal, baro, gps,
refGrid.utmZone(), refGrid.southernHemisphere());
refGrid.utmZone(), refGrid.utmLatitudeBandLetter(), refGrid.southernHemisphere());
return CoordinateConverter.toGeographic(decl, MapDatum.WGS84, grid);
} catch (Exception e) {
return null;