diff --git a/src/main/java/dev/coph/flightscore/backend/coordinate/Coordinate.java b/src/main/java/dev/coph/flightscore/backend/coordinate/Coordinate.java index da21d0f..de5dbbb 100644 --- a/src/main/java/dev/coph/flightscore/backend/coordinate/Coordinate.java +++ b/src/main/java/dev/coph/flightscore/backend/coordinate/Coordinate.java @@ -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; + } } diff --git a/src/main/java/dev/coph/flightscore/backend/map/CoordinateConverter.java b/src/main/java/dev/coph/flightscore/backend/map/CoordinateConverter.java new file mode 100644 index 0000000..0fcb407 --- /dev/null +++ b/src/main/java/dev/coph/flightscore/backend/map/CoordinateConverter.java @@ -0,0 +1,502 @@ +package dev.coph.flightscore.backend.map; + +/** + * High precision conversion utility between any combination of {@link MapDatum} + * and {@link MapGrid}. + *

+ * The conversion of a full {@code (datum, grid)} pair into any other + * {@code (datum, grid)} pair runs through three stages: + *

    + *
  1. inverse projection of the source grid onto geographic coordinates + * (on the source datum's ellipsoid),
    2. + *
  2. a geocentric 7-parameter Helmert datum transformation + * (source datum → WGS84 → target datum), carried out in 3D so + * the ellipsoidal height is transformed rigorously,
    3. + *
  3. forward projection onto the target grid (on the target datum's + * ellipsoid).
    4. + *
+ *

Accuracy of the formulas

+ * + * 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)); + } +} diff --git a/src/main/java/dev/coph/flightscore/backend/map/GeographicCoordinate.java b/src/main/java/dev/coph/flightscore/backend/map/GeographicCoordinate.java new file mode 100644 index 0000000..032c7b0 --- /dev/null +++ b/src/main/java/dev/coph/flightscore/backend/map/GeographicCoordinate.java @@ -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 + '}'; + } +} diff --git a/src/main/java/dev/coph/flightscore/backend/map/GridCoordinate.java b/src/main/java/dev/coph/flightscore/backend/map/GridCoordinate.java new file mode 100644 index 0000000..24089a9 --- /dev/null +++ b/src/main/java/dev/coph/flightscore/backend/map/GridCoordinate.java @@ -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. + *

+ * 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}. + *

+ * 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 : "") + + '}'; + } +} diff --git a/src/main/java/dev/coph/flightscore/backend/map/MapDatum.java b/src/main/java/dev/coph/flightscore/backend/map/MapDatum.java index baf2394..19bd81b 100644 --- a/src/main/java/dev/coph/flightscore/backend/map/MapDatum.java +++ b/src/main/java/dev/coph/flightscore/backend/map/MapDatum.java @@ -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}. + *

+ * The Helmert parameters are given as {@code datum -> WGS84} using the + * Position Vector 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()); } -} \ No newline at end of file + + /** + * 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; + } +} diff --git a/src/main/java/dev/coph/flightscore/backend/map/MapGrid.java b/src/main/java/dev/coph/flightscore/backend/map/MapGrid.java index ed05d1b..1ac0268 100644 --- a/src/main/java/dev/coph/flightscore/backend/map/MapGrid.java +++ b/src/main/java/dev/coph/flightscore/backend/map/MapGrid.java @@ -1,23 +1,47 @@ package dev.coph.flightscore.backend.map; +/** + * Projected grid systems together with the parameters required to evaluate + * their projection. + *

+ * 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; } -} \ No newline at end of file + + 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); + } +} diff --git a/src/main/java/dev/coph/flightscore/backend/map/ProjectionMethod.java b/src/main/java/dev/coph/flightscore/backend/map/ProjectionMethod.java new file mode 100644 index 0000000..6174fb4 --- /dev/null +++ b/src/main/java/dev/coph/flightscore/backend/map/ProjectionMethod.java @@ -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üger series (UTM, BNG, national TM grids). */ + TRANSVERSE_MERCATOR, + /** Swiss oblique Mercator (rigorous ellipsoid → sphere → oblique cylinder algorithm). */ + SWISS_OBLIQUE_MERCATOR +} diff --git a/src/main/java/dev/coph/flightscore/backend/track/parser/BalloonLiveParser.java b/src/main/java/dev/coph/flightscore/backend/track/parser/BalloonLiveParser.java index e92a09d..164acd9 100644 --- a/src/main/java/dev/coph/flightscore/backend/track/parser/BalloonLiveParser.java +++ b/src/main/java/dev/coph/flightscore/backend/track/parser/BalloonLiveParser.java @@ -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; } }