Add Coordinate Converter with high-precision transformations and integrate Geographic/Projected Coordinate models

- Introduced `CoordinateConverter` for high-precision datum and map grid transformations using Helmert and advanced projection algorithms.
- Added `GeographicCoordinate`, `GridCoordinate`, and `ProjectionMethod` classes to represent and manage coordinate data.
- Enhanced `MapGrid` with specific grid definitions and parameters.
- Updated `BalloonLiveParser` to include reference coordinate handling.
This commit is contained in:
Jan Meinl
2026-05-15 04:12:52 +02:00
parent 02d4935e01
commit 7bcf1bb0d6
8 changed files with 790 additions and 43 deletions
@@ -25,4 +25,11 @@ public class Coordinate {
this.barometricAltitude = altitude;
this.gpsAltitude = altitude;
}
public Coordinate(double latitude, double longitude) {
this.latitude = latitude;
this.longitude = longitude;
this.barometricAltitude = null;
this.gpsAltitude = null;
}
}
@@ -0,0 +1,502 @@
package dev.coph.flightscore.backend.map;
/**
* High precision conversion utility between any combination of {@link MapDatum}
* and {@link MapGrid}.
* <p>
* The conversion of a full {@code (datum, grid)} pair into any other
* {@code (datum, grid)} pair runs through three stages:
* <ol>
* <li>inverse projection of the source grid onto geographic coordinates
* (on the source datum's ellipsoid),</li>
* <li>a geocentric 7-parameter Helmert datum transformation
* (source datum &rarr; WGS84 &rarr; target datum), carried out in 3D so
* the ellipsoidal height is transformed rigorously,</li>
* <li>forward projection onto the target grid (on the target datum's
* ellipsoid).</li>
* </ol>
* <h2>Accuracy of the formulas</h2>
* <ul>
* <li>Transverse Mercator is evaluated with the Kr&uuml;ger series to order
* {@code n^6} (Karney's algorithm), which is accurate to well below a
* millimetre.</li>
* <li>The Swiss grid uses the rigorous swisstopo algorithm (ellipsoid
* &rarr; sphere &rarr; oblique cylinder), accurate to a few millimetres.</li>
* <li>Geographic &harr; geocentric conversions are exact (forward) and
* iterated to machine precision (inverse).</li>
* </ul>
* The overall accuracy is therefore limited only by the published
* {@link MapDatum} Helmert parameters (roughly metre level for the mean
* parameter sets stored in the enum).
*/
public final class CoordinateConverter {
private static final double ARCSEC_TO_RAD = Math.PI / (180.0 * 3600.0);
/** Iteration tolerance for the Newton/fixed-point solvers (radians). */
private static final double TOLERANCE = 1.0e-14;
private static final int MAX_ITERATIONS = 100;
private CoordinateConverter() {
}
// ------------------------------------------------------------------
// Public API
// ------------------------------------------------------------------
/**
* Converts a grid coordinate from one {@code (datum, grid)} combination into
* any other {@code (datum, grid)} combination.
*
* @param source the coordinate in the source grid
* @param sourceDatum datum the source coordinate is referenced to
* @param sourceGrid grid the source coordinate is projected with
* @param targetDatum datum the result shall be referenced to
* @param targetGrid grid the result shall be projected with
* @return the coordinate expressed in the target {@code (datum, grid)} combination
*/
public static GridCoordinate convert(GridCoordinate source,
MapDatum sourceDatum, MapGrid sourceGrid,
MapDatum targetDatum, MapGrid targetGrid) {
GeographicCoordinate sourceGeographic = toGeographic(source, sourceDatum, sourceGrid);
GeographicCoordinate targetGeographic = transformDatum(sourceGeographic, sourceDatum, targetDatum);
return toGrid(targetGeographic, targetDatum, targetGrid);
}
/**
* Converts a geographic coordinate from its datum onto any target
* {@code (datum, grid)} combination. This is the natural entry point when
* the input is plain longitude/latitude (e.g. a GPS track point) and avoids
* having to wrap it into a {@link MapGrid#LAT_LONG} {@link GridCoordinate}.
*/
public static GridCoordinate convert(GeographicCoordinate source, MapDatum sourceDatum,
MapDatum targetDatum, MapGrid targetGrid) {
return toGrid(transformDatum(source, sourceDatum, targetDatum), targetDatum, targetGrid);
}
/**
* Projects a geographic coordinate (referenced to {@code datum}) onto {@code grid}.
*/
public static GridCoordinate toGrid(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
switch (grid.projectionMethod()) {
case NONE:
return new GridCoordinate(geographic.longitude(), geographic.latitude(), geographic.ellipsoidalHeight());
case TRANSVERSE_MERCATOR:
return transverseMercatorForward(geographic, datum, grid);
case SWISS_OBLIQUE_MERCATOR:
return swissObliqueMercatorForward(geographic, datum, grid);
default:
throw new IllegalArgumentException("Unsupported projection method: " + grid.projectionMethod());
}
}
/**
* Inverse projection of a grid coordinate (projected with {@code grid},
* referenced to {@code datum}) onto geographic coordinates.
*/
public static GeographicCoordinate toGeographic(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
switch (mapGrid.projectionMethod()) {
case NONE:
return new GeographicCoordinate(grid.northing(), grid.easting(), grid.ellipsoidalHeight());
case TRANSVERSE_MERCATOR:
return transverseMercatorInverse(grid, datum, mapGrid);
case SWISS_OBLIQUE_MERCATOR:
return swissObliqueMercatorInverse(grid, datum, mapGrid);
default:
throw new IllegalArgumentException("Unsupported projection method: " + mapGrid.projectionMethod());
}
}
/**
* Transforms a geographic coordinate from one datum to another using a
* geocentric 7-parameter Helmert transformation via WGS84.
*/
public static GeographicCoordinate transformDatum(GeographicCoordinate geographic, MapDatum from, MapDatum to) {
if (from == to) {
return geographic;
}
double[] geocentric = geographicToGeocentric(geographic, from);
geocentric = helmertToWgs84(geocentric, from);
geocentric = helmertFromWgs84(geocentric, to);
return geocentricToGeographic(geocentric, to);
}
// ------------------------------------------------------------------
// Datum transformation (geocentric Helmert)
// ------------------------------------------------------------------
private static double[] geographicToGeocentric(GeographicCoordinate geographic, MapDatum datum) {
double a = datum.semiMajorAxis();
double e2 = datum.firstEccentricitySquared();
double phi = Math.toRadians(geographic.latitude());
double lambda = Math.toRadians(geographic.longitude());
double h = geographic.ellipsoidalHeight();
double sinPhi = Math.sin(phi);
double cosPhi = Math.cos(phi);
double primeVerticalRadius = a / Math.sqrt(1.0 - e2 * sinPhi * sinPhi);
double x = (primeVerticalRadius + h) * cosPhi * Math.cos(lambda);
double y = (primeVerticalRadius + h) * cosPhi * Math.sin(lambda);
double z = (primeVerticalRadius * (1.0 - e2) + h) * sinPhi;
return new double[]{x, y, z};
}
private static GeographicCoordinate geocentricToGeographic(double[] xyz, MapDatum datum) {
double a = datum.semiMajorAxis();
double e2 = datum.firstEccentricitySquared();
double x = xyz[0];
double y = xyz[1];
double z = xyz[2];
double lambda = Math.atan2(y, x);
double p = Math.hypot(x, y);
if (p < 1.0e-11) { // on or extremely close to the polar axis
double phiPole = z >= 0.0 ? Math.PI / 2.0 : -Math.PI / 2.0;
double semiMinor = a * Math.sqrt(1.0 - e2);
double height = Math.abs(z) - semiMinor;
return new GeographicCoordinate(Math.toDegrees(phiPole), Math.toDegrees(lambda), height);
}
double phi = Math.atan2(z, p * (1.0 - e2));
double height = 0.0;
for (int i = 0; i < MAX_ITERATIONS; i++) {
double sinPhi = Math.sin(phi);
double primeVerticalRadius = a / Math.sqrt(1.0 - e2 * sinPhi * sinPhi);
height = p / Math.cos(phi) - primeVerticalRadius;
double phiNext = Math.atan2(z, p * (1.0 - e2 * primeVerticalRadius / (primeVerticalRadius + height)));
if (Math.abs(phiNext - phi) < TOLERANCE) {
phi = phiNext;
break;
}
phi = phiNext;
}
return new GeographicCoordinate(Math.toDegrees(phi), Math.toDegrees(lambda), height);
}
private static double[] helmertToWgs84(double[] xyz, MapDatum datum) {
double[][] m = helmertMatrix(datum);
double x = xyz[0];
double y = xyz[1];
double z = xyz[2];
return new double[]{
datum.deltaX() + m[0][0] * x + m[0][1] * y + m[0][2] * z,
datum.deltaY() + m[1][0] * x + m[1][1] * y + m[1][2] * z,
datum.deltaZ() + m[2][0] * x + m[2][1] * y + m[2][2] * z
};
}
private static double[] helmertFromWgs84(double[] xyz, MapDatum datum) {
double[][] inverse = invert3x3(helmertMatrix(datum));
double x = xyz[0] - datum.deltaX();
double y = xyz[1] - datum.deltaY();
double z = xyz[2] - datum.deltaZ();
return new double[]{
inverse[0][0] * x + inverse[0][1] * y + inverse[0][2] * z,
inverse[1][0] * x + inverse[1][1] * y + inverse[1][2] * z,
inverse[2][0] * x + inverse[2][1] * y + inverse[2][2] * z
};
}
/**
* Builds the {@code (1 + s) * R} matrix of the Helmert transformation
* (Position Vector rotation convention).
*/
private static double[][] helmertMatrix(MapDatum datum) {
double scale = 1.0 + datum.scaleDifference() * 1.0e-6;
double rx = datum.rotationX() * ARCSEC_TO_RAD;
double ry = datum.rotationY() * ARCSEC_TO_RAD;
double rz = datum.rotationZ() * ARCSEC_TO_RAD;
return new double[][]{
{scale, -scale * rz, scale * ry},
{scale * rz, scale, -scale * rx},
{-scale * ry, scale * rx, scale}
};
}
private static double[][] invert3x3(double[][] m) {
double determinant =
m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
- m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
double inverseDeterminant = 1.0 / determinant;
double[][] inverse = new double[3][3];
inverse[0][0] = (m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inverseDeterminant;
inverse[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inverseDeterminant;
inverse[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inverseDeterminant;
inverse[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inverseDeterminant;
inverse[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inverseDeterminant;
inverse[1][2] = (m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inverseDeterminant;
inverse[2][0] = (m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inverseDeterminant;
inverse[2][1] = (m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inverseDeterminant;
inverse[2][2] = (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inverseDeterminant;
return inverse;
}
// ------------------------------------------------------------------
// Transverse Mercator (Krueger series, Karney's algorithm)
// ------------------------------------------------------------------
private static GridCoordinate transverseMercatorForward(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
double a = datum.semiMajorAxis();
double f = datum.flattening();
double e = datum.firstEccentricity();
double n = f / (2.0 - f);
double phi = Math.toRadians(geographic.latitude());
Integer zone = null;
boolean southern = false;
double centralMeridian;
double falseNorthing = grid.falseNorthing();
if (grid.isZoneDependent()) {
zone = utmZone(geographic.longitude());
southern = geographic.latitude() < 0.0;
centralMeridian = Math.toRadians(zone * 6.0 - 183.0);
falseNorthing = southern ? 10_000_000.0 : 0.0;
} else {
centralMeridian = Math.toRadians(grid.centralMeridian());
}
double k0 = grid.scaleFactor();
double falseEasting = grid.falseEasting();
double deltaLambda = Math.toRadians(geographic.longitude()) - centralMeridian;
double rectifyingRadius = rectifyingRadius(a, n);
double[] alpha = kruegerAlpha(n);
double xiOrigin = tmForwardSeries(Math.toRadians(grid.latitudeOfOrigin()), 0.0, e, alpha)[0];
double[] series = tmForwardSeries(phi, deltaLambda, e, alpha);
double xi = series[0];
double eta = series[1];
double easting = falseEasting + k0 * rectifyingRadius * eta;
double northing = falseNorthing + k0 * rectifyingRadius * (xi - xiOrigin);
return new GridCoordinate(easting, northing, geographic.ellipsoidalHeight(), zone, southern);
}
private static GeographicCoordinate transverseMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
double a = datum.semiMajorAxis();
double f = datum.flattening();
double e = datum.firstEccentricity();
double e2 = datum.firstEccentricitySquared();
double n = f / (2.0 - f);
double centralMeridian;
double falseNorthing = mapGrid.falseNorthing();
if (mapGrid.isZoneDependent()) {
if (grid.utmZone() == null) {
throw new IllegalArgumentException("A UTM grid coordinate must carry its zone for the inverse projection");
}
centralMeridian = Math.toRadians(grid.utmZone() * 6.0 - 183.0);
falseNorthing = grid.southernHemisphere() ? 10_000_000.0 : 0.0;
} else {
centralMeridian = Math.toRadians(mapGrid.centralMeridian());
}
double k0 = mapGrid.scaleFactor();
double falseEasting = mapGrid.falseEasting();
double rectifyingRadius = rectifyingRadius(a, n);
double[] alpha = kruegerAlpha(n);
double[] beta = kruegerBeta(n);
double xiOrigin = tmForwardSeries(Math.toRadians(mapGrid.latitudeOfOrigin()), 0.0, e, alpha)[0];
double xi = (grid.northing() - falseNorthing) / (k0 * rectifyingRadius) + xiOrigin;
double eta = (grid.easting() - falseEasting) / (k0 * rectifyingRadius);
double xiPrime = xi;
double etaPrime = eta;
for (int j = 1; j <= 6; j++) {
xiPrime -= beta[j] * Math.sin(2 * j * xi) * Math.cosh(2 * j * eta);
etaPrime -= beta[j] * Math.cos(2 * j * xi) * Math.sinh(2 * j * eta);
}
double conformalTangent = Math.sin(xiPrime) / Math.hypot(Math.sinh(etaPrime), Math.cos(xiPrime));
double tau = solveGeodeticTangent(conformalTangent, e, e2);
double phi = Math.atan(tau);
double lambda = centralMeridian + Math.atan2(Math.sinh(etaPrime), Math.cos(xiPrime));
return new GeographicCoordinate(Math.toDegrees(phi), Math.toDegrees(lambda), grid.ellipsoidalHeight());
}
/**
* Evaluates the Krueger forward series and returns {@code {xi, eta}} for a
* given geodetic latitude and longitude offset from the central meridian.
*/
private static double[] tmForwardSeries(double phi, double deltaLambda, double e, double[] alpha) {
double tau = Math.tan(phi);
double sigma = Math.sinh(e * atanh(e * tau / Math.sqrt(1.0 + tau * tau)));
double conformalTangent = tau * Math.sqrt(1.0 + sigma * sigma) - sigma * Math.sqrt(1.0 + tau * tau);
double xiPrime = Math.atan2(conformalTangent, Math.cos(deltaLambda));
double etaPrime = asinh(Math.sin(deltaLambda) / Math.hypot(conformalTangent, Math.cos(deltaLambda)));
double xi = xiPrime;
double eta = etaPrime;
for (int j = 1; j <= 6; j++) {
xi += alpha[j] * Math.sin(2 * j * xiPrime) * Math.cosh(2 * j * etaPrime);
eta += alpha[j] * Math.cos(2 * j * xiPrime) * Math.sinh(2 * j * etaPrime);
}
return new double[]{xi, eta};
}
/**
* Recovers the geodetic latitude tangent {@code tan(phi)} from the conformal
* latitude tangent using Karney's Newton iteration.
*/
private static double solveGeodeticTangent(double conformalTangent, double e, double e2) {
double tau = conformalTangent;
for (int i = 0; i < MAX_ITERATIONS; i++) {
double sigma = Math.sinh(e * atanh(e * tau / Math.sqrt(1.0 + tau * tau)));
double approximated = tau * Math.sqrt(1.0 + sigma * sigma) - sigma * Math.sqrt(1.0 + tau * tau);
double delta = (conformalTangent - approximated) * (1.0 + (1.0 - e2) * tau * tau)
/ ((1.0 - e2) * Math.sqrt(1.0 + tau * tau) * Math.sqrt(1.0 + approximated * approximated));
tau += delta;
if (Math.abs(delta) < TOLERANCE) {
break;
}
}
return tau;
}
private static double rectifyingRadius(double a, double n) {
double n2 = n * n;
double n4 = n2 * n2;
double n6 = n4 * n2;
return a / (1.0 + n) * (1.0 + n2 / 4.0 + n4 / 64.0 + n6 / 256.0);
}
private static double[] kruegerAlpha(double n) {
double n2 = n * n;
double n3 = n2 * n;
double n4 = n3 * n;
double n5 = n4 * n;
double n6 = n5 * n;
double[] alpha = new double[7];
alpha[1] = 1.0 / 2.0 * n - 2.0 / 3.0 * n2 + 5.0 / 16.0 * n3 + 41.0 / 180.0 * n4 - 127.0 / 288.0 * n5 + 7891.0 / 37800.0 * n6;
alpha[2] = 13.0 / 48.0 * n2 - 3.0 / 5.0 * n3 + 557.0 / 1440.0 * n4 + 281.0 / 630.0 * n5 - 1983433.0 / 1935360.0 * n6;
alpha[3] = 61.0 / 240.0 * n3 - 103.0 / 140.0 * n4 + 15061.0 / 26880.0 * n5 + 167603.0 / 181440.0 * n6;
alpha[4] = 49561.0 / 161280.0 * n4 - 179.0 / 168.0 * n5 + 6601661.0 / 7257600.0 * n6;
alpha[5] = 34729.0 / 80640.0 * n5 - 3418889.0 / 1995840.0 * n6;
alpha[6] = 212378941.0 / 319334400.0 * n6;
return alpha;
}
private static double[] kruegerBeta(double n) {
double n2 = n * n;
double n3 = n2 * n;
double n4 = n3 * n;
double n5 = n4 * n;
double n6 = n5 * n;
double[] beta = new double[7];
beta[1] = 1.0 / 2.0 * n - 2.0 / 3.0 * n2 + 37.0 / 96.0 * n3 - 1.0 / 360.0 * n4 - 81.0 / 512.0 * n5 + 96199.0 / 604800.0 * n6;
beta[2] = 1.0 / 48.0 * n2 + 1.0 / 15.0 * n3 - 437.0 / 1440.0 * n4 + 46.0 / 105.0 * n5 - 1118711.0 / 3870720.0 * n6;
beta[3] = 17.0 / 480.0 * n3 - 37.0 / 840.0 * n4 - 209.0 / 4480.0 * n5 + 5569.0 / 90720.0 * n6;
beta[4] = 4397.0 / 161280.0 * n4 - 11.0 / 504.0 * n5 - 830251.0 / 7257600.0 * n6;
beta[5] = 4583.0 / 161280.0 * n5 - 108847.0 / 3991680.0 * n6;
beta[6] = 20648693.0 / 638668800.0 * n6;
return beta;
}
private static int utmZone(double longitude) {
int zone = (int) Math.floor((longitude + 180.0) / 6.0) + 1;
return Math.max(1, Math.min(60, zone));
}
// ------------------------------------------------------------------
// Swiss oblique Mercator (rigorous swisstopo algorithm)
// ------------------------------------------------------------------
private static GridCoordinate swissObliqueMercatorForward(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
SwissConstants c = new SwissConstants(datum, grid);
double phi = Math.toRadians(geographic.latitude());
double lambda = Math.toRadians(geographic.longitude());
double s = c.alpha * Math.log(Math.tan(Math.PI / 4.0 + phi / 2.0))
- c.alpha * c.e / 2.0 * Math.log((1.0 + c.e * Math.sin(phi)) / (1.0 - c.e * Math.sin(phi)))
+ c.k;
double sphereLatitude = 2.0 * (Math.atan(Math.exp(s)) - Math.PI / 4.0);
double sphereLongitude = c.alpha * (lambda - c.lambda0);
double sinPseudoLatitude = Math.cos(c.b0) * Math.sin(sphereLatitude)
- Math.sin(c.b0) * Math.cos(sphereLatitude) * Math.cos(sphereLongitude);
double pseudoLongitude = Math.atan2(
Math.cos(sphereLatitude) * Math.sin(sphereLongitude),
Math.sin(c.b0) * Math.sin(sphereLatitude) + Math.cos(c.b0) * Math.cos(sphereLatitude) * Math.cos(sphereLongitude));
double easting = c.r * pseudoLongitude + grid.falseEasting();
double northing = c.r * atanh(sinPseudoLatitude) + grid.falseNorthing();
return new GridCoordinate(easting, northing, geographic.ellipsoidalHeight());
}
private static GeographicCoordinate swissObliqueMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
SwissConstants c = new SwissConstants(datum, mapGrid);
double pseudoLongitude = (grid.easting() - mapGrid.falseEasting()) / c.r;
double pseudoLatitude = 2.0 * (Math.atan(Math.exp((grid.northing() - mapGrid.falseNorthing()) / c.r)) - Math.PI / 4.0);
double sinSphereLatitude = Math.cos(c.b0) * Math.sin(pseudoLatitude)
+ Math.sin(c.b0) * Math.cos(pseudoLatitude) * Math.cos(pseudoLongitude);
double sphereLatitude = Math.asin(sinSphereLatitude);
double sphereLongitude = Math.atan2(
Math.cos(pseudoLatitude) * Math.sin(pseudoLongitude),
Math.cos(c.b0) * Math.cos(pseudoLatitude) * Math.cos(pseudoLongitude) - Math.sin(c.b0) * Math.sin(pseudoLatitude));
double lambda = c.lambda0 + sphereLongitude / c.alpha;
double isometric = (Math.log(Math.tan(Math.PI / 4.0 + sphereLatitude / 2.0)) - c.k) / c.alpha;
double phi = 2.0 * Math.atan(Math.exp(isometric)) - Math.PI / 2.0;
for (int i = 0; i < MAX_ITERATIONS; i++) {
double phiNext = 2.0 * Math.atan(Math.exp(isometric + c.e * atanh(c.e * Math.sin(phi)))) - Math.PI / 2.0;
if (Math.abs(phiNext - phi) < TOLERANCE) {
phi = phiNext;
break;
}
phi = phiNext;
}
return new GeographicCoordinate(Math.toDegrees(phi), Math.toDegrees(lambda), grid.ellipsoidalHeight());
}
/**
* Derived constants of the Swiss oblique Mercator projection (projection
* sphere radius, sphere distortion factor, fundamental point on the sphere
* and the integration constant), computed from the datum's ellipsoid and
* the grid's fundamental point.
*/
private static final class SwissConstants {
private final double e;
private final double r;
private final double alpha;
private final double b0;
private final double k;
private final double lambda0;
private SwissConstants(MapDatum datum, MapGrid grid) {
double a = datum.semiMajorAxis();
double e2 = datum.firstEccentricitySquared();
this.e = datum.firstEccentricity();
double phi0 = Math.toRadians(grid.latitudeOfOrigin());
this.lambda0 = Math.toRadians(grid.centralMeridian());
double sinPhi0 = Math.sin(phi0);
this.r = a * Math.sqrt(1.0 - e2) / (1.0 - e2 * sinPhi0 * sinPhi0);
this.alpha = Math.sqrt(1.0 + e2 / (1.0 - e2) * Math.pow(Math.cos(phi0), 4));
this.b0 = Math.asin(sinPhi0 / alpha);
this.k = Math.log(Math.tan(Math.PI / 4.0 + b0 / 2.0))
- alpha * Math.log(Math.tan(Math.PI / 4.0 + phi0 / 2.0))
+ alpha * e / 2.0 * Math.log((1.0 + e * sinPhi0) / (1.0 - e * sinPhi0));
}
}
// ------------------------------------------------------------------
// Math helpers (not provided by java.lang.Math)
// ------------------------------------------------------------------
private static double asinh(double x) {
return Math.log(x + Math.sqrt(x * x + 1.0));
}
private static double atanh(double x) {
return 0.5 * Math.log((1.0 + x) / (1.0 - x));
}
}
@@ -0,0 +1,36 @@
package dev.coph.flightscore.backend.map;
import lombok.Getter;
import lombok.experimental.Accessors;
/**
* A geographic coordinate (longitude / latitude) with an ellipsoidal height,
* expressed on a specific {@link MapDatum}.
*/
@Getter
@Accessors(fluent = true)
public class GeographicCoordinate {
/** Latitude in decimal degrees, positive north. */
private final double latitude;
/** Longitude in decimal degrees, positive east. */
private final double longitude;
/** Height above the reference ellipsoid in metres. */
private final double ellipsoidalHeight;
public GeographicCoordinate(double latitude, double longitude, double ellipsoidalHeight) {
this.latitude = latitude;
this.longitude = longitude;
this.ellipsoidalHeight = ellipsoidalHeight;
}
public GeographicCoordinate(double latitude, double longitude) {
this(latitude, longitude, 0.0);
}
@Override
public String toString() {
return "GeographicCoordinate{latitude=" + latitude + ", longitude=" + longitude
+ ", ellipsoidalHeight=" + ellipsoidalHeight + '}';
}
}
@@ -0,0 +1,62 @@
package dev.coph.flightscore.backend.map;
import lombok.Getter;
import lombok.experimental.Accessors;
/**
* A projected grid coordinate (easting / northing) with an ellipsoidal height.
* <p>
* For zone dependent grids (UTM) {@link #utmZone()} and {@link #southernHemisphere()}
* carry the information that cannot be recovered from easting/northing alone.
* For all other grids {@code utmZone} is {@code null} and {@code southernHemisphere}
* is {@code false}.
* <p>
* For the {@link MapGrid#LAT_LONG} grid the easting holds the longitude and the
* northing holds the latitude, both in decimal degrees.
*/
@Getter
@Accessors(fluent = true)
public class GridCoordinate {
/** Easting in metres (or longitude in degrees for {@link MapGrid#LAT_LONG}). */
private final double easting;
/** Northing in metres (or latitude in degrees for {@link MapGrid#LAT_LONG}). */
private final double northing;
/** Height above the reference ellipsoid in metres. */
private final double ellipsoidalHeight;
/** UTM zone number (1..60) for zone dependent grids, otherwise {@code null}. */
private final Integer utmZone;
/** Whether the coordinate lies on the southern hemisphere (relevant for UTM false northing). */
private final boolean southernHemisphere;
public GridCoordinate(double easting, double northing, double ellipsoidalHeight,
Integer utmZone, boolean southernHemisphere) {
this.easting = easting;
this.northing = northing;
this.ellipsoidalHeight = ellipsoidalHeight;
this.utmZone = utmZone;
this.southernHemisphere = southernHemisphere;
}
public GridCoordinate(double easting, double northing, double ellipsoidalHeight) {
this(easting, northing, ellipsoidalHeight, null, false);
}
public GridCoordinate(double easting, double northing) {
this(easting, northing, 0.0, null, false);
}
/** UTM grid coordinate carrying its zone and hemisphere. */
public static GridCoordinate utm(double easting, double northing, double ellipsoidalHeight,
int zone, boolean southernHemisphere) {
return new GridCoordinate(easting, northing, ellipsoidalHeight, zone, southernHemisphere);
}
@Override
public String toString() {
return "GridCoordinate{easting=" + easting + ", northing=" + northing
+ ", ellipsoidalHeight=" + ellipsoidalHeight
+ (utmZone != null ? ", utmZone=" + utmZone + ", southernHemisphere=" + southernHemisphere : "")
+ '}';
}
}
@@ -1,31 +1,54 @@
package dev.coph.flightscore.backend.map;
/**
* Geodetic datums together with their reference ellipsoid and the 7-parameter
* Helmert (Bursa-Wolf) transformation towards {@code WGS84}.
* <p>
* The Helmert parameters are given as {@code datum -> WGS84} using the
* <b>Position Vector</b> rotation convention. Translations are in metres,
* rotations in arc-seconds and the scale difference in parts per million (ppm).
* The values are the standard published mean parameters for each datum; they
* are accurate to roughly a metre. For survey-grade results a region specific
* datum-shift grid would be required.
*/
public enum MapDatum {
WGS84("World Geodetic System 1984", 6378137.0, 298.257223563, 0.0, 0.0),
ETRS89("European Terrestrial Reference System 1989", 6378137.0, 298.257222101, 0.0, 0.0),
NAD83("North American Datum 1983", 6378137.0, 298.257222101, 0.0, 0.0),
CH1903("Swiss MN03", 6377397.155, 299.1528128, 600000.0, 200000.0),
CH1903_PLUS("Swiss MN95", 6377397.155, 299.1528128, 2600000.0, 1200000.0),
OGB1936("OSGB 1936 / British National Grid", 6377563.396, 299.3249646, 400000.0, -100000.0),
LKS94("Lithuania 1994", 6378137.0, 298.257222101, 500000.0, 0.0),
SWEREF99("Swedish Reference Frame 1999", 6378137.0, 298.257222101, 150000.0, 0.0),
TWD67("Taiwan Datum 1967", 6378160.0, 298.25, 250000.0, 0.0),
TWD97("Taiwan Datum 1997", 6378137.0, 298.257222101, 250000.0, 0.0),
AGD1966("Australian Geodetic Datum 1966", 6378160.0, 298.25, 500000.0, 10000000.0),
ED50("European Datum 1950", 6378388.0, 297.0, 500000.0, 0.0);
// a 1/f FE FN dX dY dZ rX rY rZ ds(ppm)
WGS84("World Geodetic System 1984", 6378137.0, 298.257223563, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
ETRS89("European Terrestrial Reference System 1989", 6378137.0, 298.257222101, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
NAD83("North American Datum 1983", 6378137.0, 298.257222101, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
CH1903("Swiss MN03", 6377397.155, 299.1528128, 674.374, 15.056, 405.346, 0.0, 0.0, 0.0, 0.0),
CH1903_PLUS("Swiss MN95", 6377397.155, 299.1528128, 674.374, 15.056, 405.346, 0.0, 0.0, 0.0, 0.0),
OGB1936("OSGB 1936 / British National Grid", 6377563.396, 299.3249646, -446.448, 125.157, -542.060, -0.1502, -0.2470, -0.8421, 20.4894),
LKS94("Lithuania 1994", 6378137.0, 298.257222101, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
SWEREF99("Swedish Reference Frame 1999", 6378137.0, 298.257222101, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
TWD67("Taiwan Datum 1967", 6378160.0, 298.25, -752.130, -358.989, -179.110, 0.0, 0.0, 0.0, 0.0),
TWD97("Taiwan Datum 1997", 6378137.0, 298.257222101, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
AGD1966("Australian Geodetic Datum 1966", 6378160.0, 298.25, -133.0, -48.0, 148.0, 0.0, 0.0, 0.0, 0.0),
ED50("European Datum 1950", 6378388.0, 297.0, -87.0, -98.0, -121.0, 0.0, 0.0, 0.0, 0.0);
private final String displayName;
private final double semiMajorAxis;
private final double inverseFlattening;
private final double falseEasting;
private final double falseNorthing;
private final double deltaX;
private final double deltaY;
private final double deltaZ;
private final double rotationX;
private final double rotationY;
private final double rotationZ;
private final double scaleDifference;
MapDatum(String displayName, double semiMajorAxis, double inverseFlattening, double falseEasting, double falseNorthing) {
MapDatum(String displayName, double semiMajorAxis, double inverseFlattening,
double deltaX, double deltaY, double deltaZ, double rotationX, double rotationY, double rotationZ, double scaleDifference) {
this.displayName = displayName;
this.semiMajorAxis = semiMajorAxis;
this.inverseFlattening = inverseFlattening;
this.falseEasting = falseEasting;
this.falseNorthing = falseNorthing;
this.deltaX = deltaX;
this.deltaY = deltaY;
this.deltaZ = deltaZ;
this.rotationX = rotationX;
this.rotationY = rotationY;
this.rotationZ = rotationZ;
this.scaleDifference = scaleDifference;
}
public String displayName() {
@@ -40,14 +63,6 @@ public enum MapDatum {
return inverseFlattening;
}
public double falseEasting() {
return falseEasting;
}
public double falseNorthing() {
return falseNorthing;
}
public double flattening() {
return 1.0 / inverseFlattening;
}
@@ -55,4 +70,54 @@ public enum MapDatum {
public double semiMinorAxis() {
return semiMajorAxis * (1.0 - flattening());
}
}
/**
* First eccentricity squared {@code e^2 = f * (2 - f)} of the reference ellipsoid.
*/
public double firstEccentricitySquared() {
double f = flattening();
return f * (2.0 - f);
}
/**
* First eccentricity {@code e} of the reference ellipsoid.
*/
public double firstEccentricity() {
return Math.sqrt(firstEccentricitySquared());
}
/** Helmert translation along the X axis towards WGS84, in metres. */
public double deltaX() {
return deltaX;
}
/** Helmert translation along the Y axis towards WGS84, in metres. */
public double deltaY() {
return deltaY;
}
/** Helmert translation along the Z axis towards WGS84, in metres. */
public double deltaZ() {
return deltaZ;
}
/** Helmert rotation about the X axis towards WGS84, in arc-seconds (Position Vector convention). */
public double rotationX() {
return rotationX;
}
/** Helmert rotation about the Y axis towards WGS84, in arc-seconds (Position Vector convention). */
public double rotationY() {
return rotationY;
}
/** Helmert rotation about the Z axis towards WGS84, in arc-seconds (Position Vector convention). */
public double rotationZ() {
return rotationZ;
}
/** Helmert scale difference towards WGS84, in parts per million. */
public double scaleDifference() {
return scaleDifference;
}
}
@@ -1,23 +1,47 @@
package dev.coph.flightscore.backend.map;
/**
* Projected grid systems together with the parameters required to evaluate
* their projection.
* <p>
* Angular parameters ({@link #centralMeridian()}, {@link #latitudeOfOrigin()})
* are stored in decimal degrees. A {@link Double#NaN} central meridian marks a
* zone dependent grid (UTM): for such grids the central meridian is derived
* from the UTM zone instead.
*/
public enum MapGrid {
LAT_LONG("Lat/Long", "Geographic"),
UTM("Universal Transverse Mercator", "Transverse Mercator"),
SWISS_GRID("Swiss Grid", "Oblique Mercator"),
BRITISH_NATIONAL_GRID("British National Grid", "Transverse Mercator"),
LITHUANIA_TM("Lithuania TM", "Transverse Mercator"),
SWEREF_99_TM("SWEREF 99 TM", "Transverse Mercator"),
TWD67_TM2("TWD67 TM2", "Transverse Mercator"),
TWD97_TM2("TWD97 TM2", "Transverse Mercator"),
BNG_OGB_7("BNG OGB-7", "Transverse Mercator"),
BNG_OGB_M("BNG OGB-M", "Transverse Mercator");
// method centralMeridian latOfOrigin k0 FE FN
LAT_LONG("Lat/Long", "Geographic", ProjectionMethod.NONE, 0.0, 0.0, 1.0, 0.0, 0.0),
UTM("Universal Transverse Mercator", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, Double.NaN, 0.0, 0.9996, 500000.0, 0.0),
SWISS_GRID_LV03("Swiss Grid LV03 / MN03", "Oblique Mercator", ProjectionMethod.SWISS_OBLIQUE_MERCATOR, 7.439583333333333, 46.95240555555556, 1.0, 600000.0, 200000.0),
SWISS_GRID_LV95("Swiss Grid LV95 / MN95", "Oblique Mercator", ProjectionMethod.SWISS_OBLIQUE_MERCATOR, 7.439583333333333, 46.95240555555556, 1.0, 2600000.0, 1200000.0),
BRITISH_NATIONAL_GRID("British National Grid", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, -2.0, 49.0, 0.9996012717, 400000.0, -100000.0),
LITHUANIA_TM("Lithuania TM", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, 24.0, 0.0, 0.9998, 500000.0, 0.0),
SWEREF_99_TM("SWEREF 99 TM", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, 15.0, 0.0, 0.9996, 500000.0, 0.0),
TWD67_TM2("TWD67 TM2", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, 121.0, 0.0, 0.9999, 250000.0, 0.0),
TWD97_TM2("TWD97 TM2", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, 121.0, 0.0, 0.9999, 250000.0, 0.0),
BNG_OGB_7("BNG OGB-7", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, -2.0, 49.0, 0.9996012717, 400000.0, -100000.0),
BNG_OGB_M("BNG OGB-M", "Transverse Mercator", ProjectionMethod.TRANSVERSE_MERCATOR, -2.0, 49.0, 0.9996012717, 400000.0, -100000.0);
private final String displayName;
private final String projectionType;
private final ProjectionMethod projectionMethod;
private final double centralMeridian;
private final double latitudeOfOrigin;
private final double scaleFactor;
private final double falseEasting;
private final double falseNorthing;
MapGrid(String displayName, String projectionType) {
MapGrid(String displayName, String projectionType, ProjectionMethod projectionMethod,
double centralMeridian, double latitudeOfOrigin, double scaleFactor, double falseEasting, double falseNorthing) {
this.displayName = displayName;
this.projectionType = projectionType;
this.projectionMethod = projectionMethod;
this.centralMeridian = centralMeridian;
this.latitudeOfOrigin = latitudeOfOrigin;
this.scaleFactor = scaleFactor;
this.falseEasting = falseEasting;
this.falseNorthing = falseNorthing;
}
public String displayName() {
@@ -27,4 +51,41 @@ public enum MapGrid {
public String projectionType() {
return projectionType;
}
}
public ProjectionMethod projectionMethod() {
return projectionMethod;
}
/** Central meridian (longitude of origin) in decimal degrees, or {@link Double#NaN} for zone dependent grids. */
public double centralMeridian() {
return centralMeridian;
}
/** Latitude of origin in decimal degrees. */
public double latitudeOfOrigin() {
return latitudeOfOrigin;
}
/** Scale factor on the central meridian / origin. */
public double scaleFactor() {
return scaleFactor;
}
/** False easting added to the projected easting, in metres. */
public double falseEasting() {
return falseEasting;
}
/** False northing added to the projected northing, in metres. */
public double falseNorthing() {
return falseNorthing;
}
/**
* Whether the central meridian depends on a (UTM) zone and therefore has to
* be derived from the coordinate rather than read from {@link #centralMeridian()}.
*/
public boolean isZoneDependent() {
return Double.isNaN(centralMeridian);
}
}
@@ -0,0 +1,14 @@
package dev.coph.flightscore.backend.map;
/**
* Mathematical projection method used by a {@link MapGrid} to map geographic
* coordinates onto a planar grid.
*/
public enum ProjectionMethod {
/** No projection: the grid coordinates are geographic longitude/latitude in degrees. */
NONE,
/** Transverse Mercator, evaluated with the Kr&uuml;ger series (UTM, BNG, national TM grids). */
TRANSVERSE_MERCATOR,
/** Swiss oblique Mercator (rigorous ellipsoid &rarr; sphere &rarr; oblique cylinder algorithm). */
SWISS_OBLIQUE_MERCATOR
}
@@ -1,12 +1,12 @@
package dev.coph.flightscore.backend.track.parser;
import dev.coph.flightscore.backend.coordinate.Coordinate;
import dev.coph.flightscore.backend.track.Track;
public class BalloonLiveParser implements TrackParser{
public class BalloonLiveParser implements TrackParser {
@Override
public Track parse(String[] lines) {
public Track parse(String[] lines, Coordinate referenceCoordinate) {
return null;
}
}