|
|
|
@@ -1,6 +1,5 @@
|
|
|
|
|
package dev.coph.flightscore.backend.map;
|
|
|
|
|
|
|
|
|
|
import dev.coph.flightscore.backend.competition.AltitudeSource;
|
|
|
|
|
import dev.coph.flightscore.backend.coordinate.Altitude;
|
|
|
|
|
import dev.coph.flightscore.backend.coordinate.Coordinate;
|
|
|
|
|
|
|
|
|
@@ -8,37 +7,39 @@ import dev.coph.flightscore.backend.coordinate.Coordinate;
|
|
|
|
|
* 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
|
|
|
|
|
* The full conversion of a {@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 → WGS84 → target datum), carried out in 3D so
|
|
|
|
|
* the ellipsoidal height is transformed rigorously,</li>
|
|
|
|
|
* (source datum → WGS84 → target datum),</li>
|
|
|
|
|
* <li>forward projection onto the target grid (on the target datum's
|
|
|
|
|
* ellipsoid).</li>
|
|
|
|
|
* </ol>
|
|
|
|
|
* <h2>Altitudes</h2>
|
|
|
|
|
* The {@link Coordinate} class carries two independent altitude readings
|
|
|
|
|
* (barometric and GPS). Pure projections ({@link #toGrid} / {@link #toGeographic})
|
|
|
|
|
* leave both untouched. {@link #transformDatum} and the variants of
|
|
|
|
|
* {@link #convert} that include a datum change run the Helmert transform once
|
|
|
|
|
* per altitude that is actually set; {@code null} altitudes stay {@code null}
|
|
|
|
|
* and the height is excluded from the geocentric step (a height of {@code 0}
|
|
|
|
|
* <i>is</i> a valid altitude and is transformed normally).
|
|
|
|
|
* <h2>Accuracy of the formulas</h2>
|
|
|
|
|
* <ul>
|
|
|
|
|
* <li>Transverse Mercator is evaluated with the Krüger series to order
|
|
|
|
|
* {@code n^6} (Karney's algorithm), which is accurate to well below a
|
|
|
|
|
* millimetre.</li>
|
|
|
|
|
* {@code n^6} (Karney's algorithm).</li>
|
|
|
|
|
* <li>The Swiss grid uses the rigorous swisstopo algorithm (ellipsoid
|
|
|
|
|
* → sphere → oblique cylinder), accurate to a few millimetres.</li>
|
|
|
|
|
* → sphere → oblique cylinder).</li>
|
|
|
|
|
* <li>Geographic ↔ 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).
|
|
|
|
|
* {@link MapDatum} Helmert parameters.
|
|
|
|
|
*/
|
|
|
|
|
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;
|
|
|
|
|
|
|
|
|
@@ -52,74 +53,38 @@ public final class CoordinateConverter {
|
|
|
|
|
/**
|
|
|
|
|
* 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);
|
|
|
|
|
Coordinate geographic = toGeographic(source, sourceDatum, sourceGrid);
|
|
|
|
|
Coordinate shifted = transformDatum(geographic, sourceDatum, targetDatum);
|
|
|
|
|
return toGrid(shifted, 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}.
|
|
|
|
|
* Converts a geographic {@link Coordinate} from its datum onto any target
|
|
|
|
|
* {@code (datum, grid)} combination. Natural entry point for plain GPS
|
|
|
|
|
* track points (WGS84 lon/lat).
|
|
|
|
|
*/
|
|
|
|
|
public static GridCoordinate convert(GeographicCoordinate source, MapDatum sourceDatum,
|
|
|
|
|
public static GridCoordinate convert(Coordinate source, MapDatum sourceDatum,
|
|
|
|
|
MapDatum targetDatum, MapGrid targetGrid) {
|
|
|
|
|
return toGrid(transformDatum(source, sourceDatum, targetDatum), targetDatum, targetGrid);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Projects a {@link Coordinate} (referenced to {@code datum}) onto {@code grid}.
|
|
|
|
|
* <p>
|
|
|
|
|
* The {@link Coordinate#gpsAltitude() GPS altitude} is taken as the ellipsoidal
|
|
|
|
|
* height since it is the GNSS-derived height and therefore the closest match;
|
|
|
|
|
* the barometric altitude is used as a fallback and {@code 0} if neither is set.
|
|
|
|
|
* Both altitudes are carried through unchanged.
|
|
|
|
|
*/
|
|
|
|
|
public static GridCoordinate toGrid(Coordinate coordinate, AltitudeSource altitudeSource, MapDatum datum, MapGrid grid) {
|
|
|
|
|
return toGrid(toGeographicCoordinate(coordinate, altitudeSource), datum, grid);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Inverse projection of a grid coordinate onto a {@link Coordinate}. The
|
|
|
|
|
* resulting ellipsoidal height is stored as the coordinate's altitude.
|
|
|
|
|
*/
|
|
|
|
|
public static Coordinate toCoordinate(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
GeographicCoordinate geographic = toGeographic(grid, datum, mapGrid);
|
|
|
|
|
return new Coordinate(geographic.latitude(), geographic.longitude(),
|
|
|
|
|
Altitude.fromMeters(geographic.ellipsoidalHeight()));
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static GeographicCoordinate toGeographicCoordinate(Coordinate coordinate, AltitudeSource source) {
|
|
|
|
|
Altitude altitude = switch (source) {
|
|
|
|
|
case GPS -> coordinate.gpsAltitude();
|
|
|
|
|
case BAROMETRIC -> coordinate.barometricAltitude();
|
|
|
|
|
};
|
|
|
|
|
double height = altitude != null ? altitude.meters() : 0.0;
|
|
|
|
|
return new GeographicCoordinate(coordinate.latitude(), coordinate.longitude(), height);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Projects a geographic coordinate (referenced to {@code datum}) onto {@code grid}.
|
|
|
|
|
*/
|
|
|
|
|
public static GridCoordinate toGrid(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
|
|
|
|
|
public static GridCoordinate toGrid(Coordinate coordinate, MapDatum datum, MapGrid grid) {
|
|
|
|
|
switch (grid.projectionMethod()) {
|
|
|
|
|
case NONE:
|
|
|
|
|
return new GridCoordinate(geographic.longitude(), geographic.latitude(), geographic.ellipsoidalHeight());
|
|
|
|
|
return new GridCoordinate(coordinate.longitude(), coordinate.latitude(),
|
|
|
|
|
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, false);
|
|
|
|
|
case TRANSVERSE_MERCATOR:
|
|
|
|
|
return transverseMercatorForward(geographic, datum, grid);
|
|
|
|
|
return transverseMercatorForward(coordinate, datum, grid);
|
|
|
|
|
case SWISS_OBLIQUE_MERCATOR:
|
|
|
|
|
return swissObliqueMercatorForward(geographic, datum, grid);
|
|
|
|
|
return swissObliqueMercatorForward(coordinate, datum, grid);
|
|
|
|
|
default:
|
|
|
|
|
throw new IllegalArgumentException("Unsupported projection method: " + grid.projectionMethod());
|
|
|
|
|
}
|
|
|
|
@@ -127,12 +92,14 @@ public final class CoordinateConverter {
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Inverse projection of a grid coordinate (projected with {@code grid},
|
|
|
|
|
* referenced to {@code datum}) onto geographic coordinates.
|
|
|
|
|
* referenced to {@code datum}) onto a geographic {@link Coordinate}. Both
|
|
|
|
|
* altitudes are carried through unchanged.
|
|
|
|
|
*/
|
|
|
|
|
public static GeographicCoordinate toGeographic(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
public static Coordinate toGeographic(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
switch (mapGrid.projectionMethod()) {
|
|
|
|
|
case NONE:
|
|
|
|
|
return new GeographicCoordinate(grid.northing(), grid.easting(), grid.ellipsoidalHeight());
|
|
|
|
|
return new Coordinate(grid.northing(), grid.easting(),
|
|
|
|
|
grid.barometricAltitude(), grid.gpsAltitude());
|
|
|
|
|
case TRANSVERSE_MERCATOR:
|
|
|
|
|
return transverseMercatorInverse(grid, datum, mapGrid);
|
|
|
|
|
case SWISS_OBLIQUE_MERCATOR:
|
|
|
|
@@ -143,41 +110,73 @@ public final class CoordinateConverter {
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Transforms a geographic coordinate from one datum to another using a
|
|
|
|
|
* geocentric 7-parameter Helmert transformation via WGS84.
|
|
|
|
|
* Transforms a {@link Coordinate} from one datum to another using a
|
|
|
|
|
* geocentric 7-parameter Helmert transformation via WGS84. Each set
|
|
|
|
|
* altitude is transformed independently; {@code null} altitudes remain
|
|
|
|
|
* {@code null}.
|
|
|
|
|
*/
|
|
|
|
|
public static GeographicCoordinate transformDatum(GeographicCoordinate geographic, MapDatum from, MapDatum to) {
|
|
|
|
|
public static Coordinate transformDatum(Coordinate coordinate, MapDatum from, MapDatum to) {
|
|
|
|
|
if (from == to) {
|
|
|
|
|
return geographic;
|
|
|
|
|
return coordinate;
|
|
|
|
|
}
|
|
|
|
|
double[] geocentric = geographicToGeocentric(geographic, from);
|
|
|
|
|
geocentric = helmertToWgs84(geocentric, from);
|
|
|
|
|
geocentric = helmertFromWgs84(geocentric, to);
|
|
|
|
|
return geocentricToGeographic(geocentric, to);
|
|
|
|
|
Altitude baro = coordinate.barometricAltitude();
|
|
|
|
|
Altitude gps = coordinate.gpsAltitude();
|
|
|
|
|
|
|
|
|
|
// Use any present altitude to compute the (essentially h-independent) horizontal result.
|
|
|
|
|
double referenceHeight = gps != null ? gps.meters() : (baro != null ? baro.meters() : 0.0);
|
|
|
|
|
double[] reference = transformGeographicPosition(coordinate.latitude(), coordinate.longitude(),
|
|
|
|
|
referenceHeight, from, to);
|
|
|
|
|
|
|
|
|
|
Altitude newBaro = null;
|
|
|
|
|
Altitude newGps = null;
|
|
|
|
|
if (gps != null && baro != null && gps != baro) {
|
|
|
|
|
// Both set, with different values: run the second transform for the other altitude.
|
|
|
|
|
newGps = Altitude.fromMeters(reference[2]);
|
|
|
|
|
double[] secondary = transformGeographicPosition(coordinate.latitude(), coordinate.longitude(),
|
|
|
|
|
baro.meters(), from, to);
|
|
|
|
|
newBaro = Altitude.fromMeters(secondary[2]);
|
|
|
|
|
} else if (gps != null) {
|
|
|
|
|
newGps = Altitude.fromMeters(reference[2]);
|
|
|
|
|
if (baro != null) {
|
|
|
|
|
newBaro = newGps; // shared instance, both pointed to the same value originally
|
|
|
|
|
}
|
|
|
|
|
} else if (baro != null) {
|
|
|
|
|
newBaro = Altitude.fromMeters(reference[2]);
|
|
|
|
|
}
|
|
|
|
|
return new Coordinate(reference[0], reference[1], newBaro, newGps);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// ------------------------------------------------------------------
|
|
|
|
|
// Datum transformation (geocentric Helmert)
|
|
|
|
|
// ------------------------------------------------------------------
|
|
|
|
|
|
|
|
|
|
private static double[] geographicToGeocentric(GeographicCoordinate geographic, MapDatum datum) {
|
|
|
|
|
/** Returns {@code {latDeg, lonDeg, height}}. */
|
|
|
|
|
private static double[] transformGeographicPosition(double latDeg, double lonDeg, double height,
|
|
|
|
|
MapDatum from, MapDatum to) {
|
|
|
|
|
double[] xyz = geographicToGeocentric(latDeg, lonDeg, height, from);
|
|
|
|
|
xyz = helmertToWgs84(xyz, from);
|
|
|
|
|
xyz = helmertFromWgs84(xyz, to);
|
|
|
|
|
return geocentricToGeographic(xyz, to);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static double[] geographicToGeocentric(double latDeg, double lonDeg, double height, 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 phi = Math.toRadians(latDeg);
|
|
|
|
|
double lambda = Math.toRadians(lonDeg);
|
|
|
|
|
|
|
|
|
|
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;
|
|
|
|
|
double x = (primeVerticalRadius + height) * cosPhi * Math.cos(lambda);
|
|
|
|
|
double y = (primeVerticalRadius + height) * cosPhi * Math.sin(lambda);
|
|
|
|
|
double z = (primeVerticalRadius * (1.0 - e2) + height) * sinPhi;
|
|
|
|
|
return new double[]{x, y, z};
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static GeographicCoordinate geocentricToGeographic(double[] xyz, MapDatum datum) {
|
|
|
|
|
/** Returns {@code {latDeg, lonDeg, height}}. */
|
|
|
|
|
private static double[] geocentricToGeographic(double[] xyz, MapDatum datum) {
|
|
|
|
|
double a = datum.semiMajorAxis();
|
|
|
|
|
double e2 = datum.firstEccentricitySquared();
|
|
|
|
|
double x = xyz[0];
|
|
|
|
@@ -187,11 +186,11 @@ public final class CoordinateConverter {
|
|
|
|
|
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
|
|
|
|
|
if (p < 1.0e-11) {
|
|
|
|
|
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);
|
|
|
|
|
return new double[]{Math.toDegrees(phiPole), Math.toDegrees(lambda), height};
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
double phi = Math.atan2(z, p * (1.0 - e2));
|
|
|
|
@@ -207,7 +206,7 @@ public final class CoordinateConverter {
|
|
|
|
|
}
|
|
|
|
|
phi = phiNext;
|
|
|
|
|
}
|
|
|
|
|
return new GeographicCoordinate(Math.toDegrees(phi), Math.toDegrees(lambda), height);
|
|
|
|
|
return new double[]{Math.toDegrees(phi), Math.toDegrees(lambda), height};
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static double[] helmertToWgs84(double[] xyz, MapDatum datum) {
|
|
|
|
@@ -234,10 +233,6 @@ public final class CoordinateConverter {
|
|
|
|
|
};
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* 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;
|
|
|
|
@@ -273,21 +268,21 @@ public final class CoordinateConverter {
|
|
|
|
|
// Transverse Mercator (Krueger series, Karney's algorithm)
|
|
|
|
|
// ------------------------------------------------------------------
|
|
|
|
|
|
|
|
|
|
private static GridCoordinate transverseMercatorForward(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
|
|
|
|
|
private static GridCoordinate transverseMercatorForward(Coordinate coordinate, 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());
|
|
|
|
|
double phi = Math.toRadians(coordinate.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;
|
|
|
|
|
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 {
|
|
|
|
@@ -295,22 +290,20 @@ public final class CoordinateConverter {
|
|
|
|
|
}
|
|
|
|
|
double k0 = grid.scaleFactor();
|
|
|
|
|
double falseEasting = grid.falseEasting();
|
|
|
|
|
double deltaLambda = Math.toRadians(geographic.longitude()) - centralMeridian;
|
|
|
|
|
double deltaLambda = Math.toRadians(coordinate.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);
|
|
|
|
|
double easting = falseEasting + k0 * rectifyingRadius * series[1];
|
|
|
|
|
double northing = falseNorthing + k0 * rectifyingRadius * (series[0] - xiOrigin);
|
|
|
|
|
return new GridCoordinate(easting, northing,
|
|
|
|
|
coordinate.barometricAltitude(), coordinate.gpsAltitude(), zone, southern);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static GeographicCoordinate transverseMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
private static Coordinate transverseMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
double a = datum.semiMajorAxis();
|
|
|
|
|
double f = datum.flattening();
|
|
|
|
|
double e = datum.firstEccentricity();
|
|
|
|
@@ -351,13 +344,10 @@ public final class CoordinateConverter {
|
|
|
|
|
|
|
|
|
|
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());
|
|
|
|
|
return new Coordinate(Math.toDegrees(phi), Math.toDegrees(lambda),
|
|
|
|
|
grid.barometricAltitude(), grid.gpsAltitude());
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* 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)));
|
|
|
|
@@ -375,10 +365,6 @@ public final class CoordinateConverter {
|
|
|
|
|
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++) {
|
|
|
|
@@ -442,11 +428,11 @@ public final class CoordinateConverter {
|
|
|
|
|
// Swiss oblique Mercator (rigorous swisstopo algorithm)
|
|
|
|
|
// ------------------------------------------------------------------
|
|
|
|
|
|
|
|
|
|
private static GridCoordinate swissObliqueMercatorForward(GeographicCoordinate geographic, MapDatum datum, MapGrid grid) {
|
|
|
|
|
private static GridCoordinate swissObliqueMercatorForward(Coordinate coordinate, MapDatum datum, MapGrid grid) {
|
|
|
|
|
SwissConstants c = new SwissConstants(datum, grid);
|
|
|
|
|
|
|
|
|
|
double phi = Math.toRadians(geographic.latitude());
|
|
|
|
|
double lambda = Math.toRadians(geographic.longitude());
|
|
|
|
|
double phi = Math.toRadians(coordinate.latitude());
|
|
|
|
|
double lambda = Math.toRadians(coordinate.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)))
|
|
|
|
@@ -462,10 +448,11 @@ public final class CoordinateConverter {
|
|
|
|
|
|
|
|
|
|
double easting = c.r * pseudoLongitude + grid.falseEasting();
|
|
|
|
|
double northing = c.r * atanh(sinPseudoLatitude) + grid.falseNorthing();
|
|
|
|
|
return new GridCoordinate(easting, northing, geographic.ellipsoidalHeight());
|
|
|
|
|
return new GridCoordinate(easting, northing,
|
|
|
|
|
coordinate.barometricAltitude(), coordinate.gpsAltitude(), null, false);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
private static GeographicCoordinate swissObliqueMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
private static Coordinate swissObliqueMercatorInverse(GridCoordinate grid, MapDatum datum, MapGrid mapGrid) {
|
|
|
|
|
SwissConstants c = new SwissConstants(datum, mapGrid);
|
|
|
|
|
|
|
|
|
|
double pseudoLongitude = (grid.easting() - mapGrid.falseEasting()) / c.r;
|
|
|
|
@@ -490,15 +477,10 @@ public final class CoordinateConverter {
|
|
|
|
|
}
|
|
|
|
|
phi = phiNext;
|
|
|
|
|
}
|
|
|
|
|
return new GeographicCoordinate(Math.toDegrees(phi), Math.toDegrees(lambda), grid.ellipsoidalHeight());
|
|
|
|
|
return new Coordinate(Math.toDegrees(phi), Math.toDegrees(lambda),
|
|
|
|
|
grid.barometricAltitude(), grid.gpsAltitude());
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* 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;
|
|
|
|
|