This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
commit 4956bf7511f6840524685a3c021b96b871b1878b
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Sat Oct 28 23:29:32 2023 +0200
Following the removal of `InterpolatedMolodenskyTransform`, retrofit
`MolodenskyFormula` into `MolodenskyTransform`.
---
.../main/module-info.java | 2 +-
.../transform/AbridgedMolodenskyTransform2D.java | 2 +-
.../operation/transform/MolodenskyFormula.java | 463 ---------------------
.../operation/transform/MolodenskyTransform.java | 418 ++++++++++++++++++-
4 files changed, 401 insertions(+), 484 deletions(-)
diff --git a/endorsed/src/org.apache.sis.referencing/main/module-info.java
b/endorsed/src/org.apache.sis.referencing/main/module-info.java
index 6b668db7a6..af2202a759 100644
--- a/endorsed/src/org.apache.sis.referencing/main/module-info.java
+++ b/endorsed/src/org.apache.sis.referencing/main/module-info.java
@@ -21,7 +21,7 @@
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @author Maxime Gavens (Geomatys)
- * @version 1.4
+ * @version 1.5
* @since 0.3
*/
module org.apache.sis.referencing {
diff --git
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbridgedMolodenskyTransform2D.java
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbridgedMolodenskyTransform2D.java
index d1ca0e7fd7..0eaef9b946 100644
---
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbridgedMolodenskyTransform2D.java
+++
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbridgedMolodenskyTransform2D.java
@@ -26,7 +26,7 @@ import static java.lang.Math.*;
/**
* Two-dimensional abridged Molodensky transform with all translation terms
fixed to zero.
* This implementation performs only a change of ellipsoid. It provides
nothing new compared
- * to {@link MolodenskyFormula}, except performance.
+ * to {@link MolodenskyTransform}, except performance.
*
* <p><b>Note:</b> this transform is yet more abridged than standard "abridged
Molondensky" transform since
* it sets all translation terms to zero. A better class name could be "Very
abridged Molodensky transform".
diff --git
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyFormula.java
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyFormula.java
deleted file mode 100644
index 63c0a9f3e9..0000000000
---
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyFormula.java
+++ /dev/null
@@ -1,463 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package org.apache.sis.referencing.operation.transform;
-
-import javax.measure.Unit;
-import javax.measure.quantity.Length;
-import org.opengis.parameter.ParameterValueGroup;
-import org.opengis.parameter.ParameterDescriptorGroup;
-import org.opengis.referencing.datum.Ellipsoid;
-import org.opengis.referencing.operation.Matrix;
-import org.opengis.referencing.operation.TransformException;
-import org.apache.sis.referencing.operation.matrix.Matrices;
-import org.apache.sis.referencing.datum.DefaultEllipsoid;
-import org.apache.sis.referencing.datum.DatumShiftGrid;
-import org.apache.sis.referencing.operation.provider.Molodensky;
-import org.apache.sis.util.ArgumentChecks;
-import org.apache.sis.util.ComparisonMode;
-import org.apache.sis.util.Debug;
-import org.apache.sis.util.internal.Numerics;
-import org.apache.sis.parameter.Parameters;
-
-import static java.lang.Math.*;
-
-
-/**
- * Implementation of Molodensky formulas.
- *
- * @author Martin Desruisseaux (Geomatys)
- */
-abstract class MolodenskyFormula extends DatumShiftTransform {
- /**
- * Serial number for inter-operability with different versions.
- */
- private static final long serialVersionUID = 7684676923384073055L;
-
- /**
- * The value of 1/sin(1″) multiplied by the conversion factor from
arc-seconds to radians (π/180)/(60⋅60).
- * This is the final multiplication factor for Δλ and Δφ.
- */
- static final double ANGULAR_SCALE = 1.00000000000391744;
-
- /**
- * {@code true} if the source coordinates have a height.
- */
- final boolean isSource3D;
-
- /**
- * {@code true} if the target coordinates have a height.
- */
- final boolean isTarget3D;
-
- /**
- * {@code true} for the abridged formula, or {@code false} for the
complete one.
- */
- final boolean isAbridged;
-
- /**
- * Shift along the geocentric X axis (toward prime meridian)
- * in units of the semi-major axis of the source ellipsoid.
- *
- * @see org.apache.sis.referencing.datum.BursaWolfParameters#tX
- */
- protected final double tX;
-
- /**
- * Shift along the geocentric Y axis (toward 90°E)
- * in units of the semi-major axis of the source ellipsoid.
- *
- * @see org.apache.sis.referencing.datum.BursaWolfParameters#tY
- */
- protected final double tY;
-
- /**
- * Shift along the geocentric Z axis (toward north pole)
- * in units of the semi-major axis of the source ellipsoid.
- *
- * @see org.apache.sis.referencing.datum.BursaWolfParameters#tZ
- */
- protected final double tZ;
-
- /**
- * Difference in the semi-major axes of the target and source ellipsoids:
{@code Δa = target a - source a}.
- *
- * @see DefaultEllipsoid#semiMajorAxisDifference(Ellipsoid)
- */
- final double Δa;
-
- /**
- * Difference between the flattening of the target and source ellipsoids
(Δf), opportunistically modified
- * with additional terms. The value depends on whether this Molodensky
transform is abridged or not:
- *
- * <ul>
- * <li>For Molodensky, this field is set to (b⋅Δf).</li>
- * <li>For Abridged Molodensky, this field is set to (a⋅Δf) +
(f⋅Δa).</li>
- * </ul>
- *
- * where Δf = <var>target flattening</var> - <var>source flattening</var>.
- */
- final double Δfmod;
-
- /**
- * Semi-major axis length (<var>a</var>) of the source ellipsoid.
- */
- protected final double semiMajor;
-
- /**
- * The square of eccentricity of the source ellipsoid.
- * This can be computed by ℯ² = (a²-b²)/a² where
- * <var>a</var> is the <cite>semi-major</cite> axis length and
- * <var>b</var> is the <cite>semi-minor</cite> axis length.
- *
- * @see DefaultEllipsoid#getEccentricitySquared()
- */
- protected final double eccentricitySquared;
-
- /**
- * Constructs the inverse of a Molodensky transform.
- *
- * @param inverse the transform for which to create the inverse.
- * @param source the source ellipsoid of the given {@code inverse}
transform.
- * @param target the target ellipsoid of the given {@code inverse}
transform.
- * @param descriptor the contextual parameter descriptor.
- */
- MolodenskyFormula(final MolodenskyFormula inverse, final Ellipsoid source,
final Ellipsoid target,
- final ParameterDescriptorGroup descriptor)
- {
- this(target, inverse.isTarget3D,
- source, inverse.isSource3D,
- -inverse.tX, -inverse.tY, -inverse.tZ, inverse.grid,
inverse.isAbridged, descriptor);
- }
-
- /**
- * Creates a Molodensky transform from the specified parameters.
- * If a non-null {@code grid} is specified, it is caller's responsibility
to verify its validity.
- *
- * @param source the source ellipsoid.
- * @param isSource3D {@code true} if the source coordinates have a height.
- * @param target the target ellipsoid.
- * @param isTarget3D {@code true} if the target coordinates have a height.
- * @param tX the geocentric <var>X</var> translation in same
units than the source ellipsoid axes.
- * @param tY the geocentric <var>Y</var> translation in same
units than the source ellipsoid axes.
- * @param tZ the geocentric <var>Z</var> translation in same
units than the source ellipsoid axes.
- * @param grid interpolation grid in geocentric coordinates, or
{@code null} if none.
- * @param isAbridged {@code true} for the abridged formula, or {@code
false} for the complete one.
- * @param descriptor the contextual parameter descriptor.
- */
- @SuppressWarnings("OverridableMethodCallInConstructor")
- MolodenskyFormula(final Ellipsoid source, final boolean isSource3D,
- final Ellipsoid target, final boolean isTarget3D,
- final double tX, final double tY, final double tZ,
- final DatumShiftGrid<?,?> grid, final boolean isAbridged,
- final ParameterDescriptorGroup descriptor)
- {
- super(descriptor, isSource3D, isTarget3D, grid);
- ArgumentChecks.ensureNonNull("source", source);
- ArgumentChecks.ensureNonNull("target", target);
- final DefaultEllipsoid src = DefaultEllipsoid.castOrCopy(source);
- this.isSource3D = isSource3D;
- this.isTarget3D = isTarget3D;
- this.isAbridged = isAbridged;
- this.semiMajor = src.getSemiMajorAxis();
- this.Δa = src.semiMajorAxisDifference(target);
- this.tX = tX;
- this.tY = tY;
- this.tZ = tZ;
-
- final double semiMinor = src.getSemiMinorAxis();
- final double Δf = src.flatteningDifference(target);
- eccentricitySquared = src.getEccentricitySquared();
- Δfmod = isAbridged ? (semiMajor * Δf) + (semiMajor - semiMinor) * (Δa
/ semiMajor)
- : (semiMinor * Δf);
- /*
- * Copy parameters to the ContextualParameter. Those parameters are
not used directly
- * by MolodenskyTransform, but we need to store them in case the user
asks for them.
- */
- final Unit<Length> unit = src.getAxisUnit();
- setContextParameters(semiMajor, semiMinor, unit, target);
- completeParameters(context, semiMinor, unit, Δf);
- /*
- * Prepare two affine transforms to be executed before and after the
MolodenskyTransform:
- *
- * - A "normalization" transform for converting degrees to radians,
- * - A "denormalization" transform for for converting radians to
degrees.
- */
- context.normalizeGeographicInputs(0);
- context.denormalizeGeographicOutputs(0);
- }
-
- /**
- * Returns a copy of internal parameter values of this transform.
- * The returned group contains parameters for the source ellipsoid
semi-axis lengths
- * and the differences between source and target ellipsoid parameters.
- *
- * <h4>Usage note</h4>
- * This method is mostly for {@linkplain
org.apache.sis.io.wkt.Convention#INTERNAL debugging purposes}
- * since the isolation of non-linear parameters in this class is highly
implementation dependent.
- * Most GIS applications will instead be interested in the {@linkplain
#getContextualParameters()
- * contextual parameters}.
- *
- * @return a copy of the internal parameter values for this transform.
- */
- @Debug
- @Override
- public ParameterValueGroup getParameterValues() {
- final Unit<?> unit =
context.getOrCreate(Molodensky.SRC_SEMI_MAJOR).getUnit();
- final double semiMinor =
context.getOrCreate(Molodensky.SRC_SEMI_MINOR).doubleValue(unit);
- final Parameters pg =
Parameters.castOrWrap(getParameterDescriptors().createValue());
- pg.getOrCreate(Molodensky.SRC_SEMI_MAJOR).setValue(semiMajor, unit);
- pg.getOrCreate(Molodensky.SRC_SEMI_MINOR).setValue(semiMinor, unit);
- completeParameters(pg, semiMinor, unit, Double.NaN);
- return pg;
- }
-
- /**
- * Sets parameter values in the given group for parameters other than axis
lengths.
- * This method is invoked for both completing contextual parameters
({@code pg == context}) and
- * for completing internal parameters ({@code pg != context}). When this
method is invoked, the
- * following parameters are already set:
- *
- * <ul>
- * <li>{@code "src_semi_major"}</li>
- * <li>{@code "src_semi_minor"}</li>
- * <li>{@code "tgt_semi_major"} (contextual parameters only)</li>
- * <li>{@code "tgt_semi_minor"} (contextual parameters only)</li>
- * </ul>
- *
- * The default implementation sets the {@code "dim"} parameter.
- * Subclasses shall override this method for completing also the following
parameters:
- *
- * <ul>
- * <li><cite>"X-axis translation"</cite> (Molodensky only)</li>
- * <li><cite>"Y-axis translation"</cite> (Molodensky only)</li>
- * <li><cite>"Z-axis translation"</cite> (Molodensky only)</li>
- * <li><cite>"Geocentric translation file"</cite> (Geocentric
interpolations only)</li>
- * <li><cite>"Semi-major axis length difference"</cite> (Always for
Molodensky, internal WKT only for geocentric interpolations)</li>
- * <li><cite>"Flattening difference"</cite> (Always for Molodensky,
internal WKT only for geocentric interpolations)</li>
- * </ul>
- *
- * @param pg where to set the parameters.
- * @param semiMinor the semi minor axis length, in unit of {@code unit}.
- * @param unit the unit of measurement to declare.
- * @param Δf the flattening difference to set, or NaN if this
method should fetch that value itself.
- */
- void completeParameters(final Parameters pg, final double semiMinor, final
Unit<?> unit, final double Δf) {
- /*
- * Unconditionally set the "dim" parameters to the number of source
dimensions (do not check for consistency
- * with the number of target dimensions) because source dimensions
determine the value of ellipsoidal heights,
- * which may change the horizontal numerical values. By contrast, the
number of target dimensions does not have
- * any impact on numerical values (it can just causes a drop of the
third value).
- */
- pg.getOrCreate(Molodensky.DIMENSION).setValue(getSourceDimensions());
- }
-
- /**
- * Gets the dimension of input points.
- *
- * @return the input dimension, which is 2 or 3.
- */
- @Override
- public final int getSourceDimensions() {
- return isSource3D ? 3 : 2;
- }
-
- /**
- * Gets the dimension of output points.
- *
- * @return the output dimension, which is 2 or 3.
- */
- @Override
- public final int getTargetDimensions() {
- return isTarget3D ? 3 : 2;
- }
-
- /**
- * Implementation of {@link #transform(double[], int, double[], int,
boolean)} with possibility
- * to override some field values. In this method signature, parameters
having the same name than
- * fields have the same value.
- *
- * @param λ longitude (radians).
- * @param φ latitude (radians).
- * @param h height above the ellipsoid in unit of semi-major
axis.
- * @param dstPts the array into which the transformed coordinate is
returned, or {@code null}.
- * @param dstOff the offset to the location of the transformed point
that is stored in the destination array.
- * @param tX the {@link #tX} field value (or a slightly
different value during geocentric interpolation).
- * @param tY the {@link #tY} field value (or a slightly
different value during geocentric interpolation).
- * @param tZ the {@link #tZ} field value (or a slightly
different value during geocentric interpolation).
- * @param derivate {@code true} for computing the derivative, or
{@code false} if not needed.
- * @throws TransformException if a point cannot be transformed.
- */
- final Matrix transform(final double λ, final double φ, final double h,
final double[] dstPts, int dstOff,
- double tX, double tY, double tZ, final boolean
derivate)
- throws TransformException
- {
- /*
- * Abridged Molodensky formulas from EPSG guidance note:
- *
- * ν = a / √(1 - ℯ²⋅sin²φ) : radius of
curvature in the prime vertical
- * ρ = a⋅(1 – ℯ²) / (1 – ℯ²⋅sin²φ)^(3/2) : radius of
curvature in the meridian
- * Δλ″ = (-tX⋅sinλ + tY⋅cosλ) / (ν⋅cosφ⋅sin1″)
- * Δφ″ = (-tX⋅sinφ⋅cosλ - tY⋅sinφ⋅sinλ + tZ⋅cosφ + [a⋅Δf +
f⋅Δa]⋅sin(2φ)) / (ρ⋅sin1″)
- * Δh = tX⋅cosφ⋅cosλ + tY⋅cosφ⋅sinλ + tZ⋅sinφ + (a⋅Δf +
f⋅Δa)⋅sin²φ - Δa
- *
- * we set:
- *
- * dfm = (a⋅Δf + f⋅Δa) in abridged case (b⋅Δf in non-abridged
case)
- * sin(2φ) = 2⋅sin(φ)⋅cos(φ)
- */
- final double sinλ = sin(λ);
- final double cosλ = cos(λ);
- final double sinφ = sin(φ);
- final double cosφ = cos(φ);
- final double sin2φ = sinφ * sinφ;
- final double ν2den = 1 - eccentricitySquared*sin2φ; //
Square of the denominator of ν
- final double νden = sqrt(ν2den); //
Denominator of ν
- final double ρden = ν2den * νden; //
Denominator of ρ
- double ρ = semiMajor * (1 - eccentricitySquared) / ρden; //
Other notation: Rm = ρ
- double ν = semiMajor / νden; //
Other notation: Rn = ν
- double t = Δfmod * 2; //
A term in the calculation of Δφ
- if (!isAbridged) {
- ρ += h;
- ν += h;
- t = t*(0.5/νden + 0.5/ρden) // = Δf⋅[ν⋅(b/a) +
ρ⋅(a/b)] (without the +h in ν and ρ)
- + Δa*eccentricitySquared/νden; // = Δa⋅[ℯ²⋅ν/a]
- }
- final double spcλ = tY*sinλ + tX*cosλ; // "spc"
stands for "sin plus cos"
- final double cmsλ = tY*cosλ - tX*sinλ; // "cms"
stands for "cos minus sin"
- final double cmsφ = (tZ + t*sinφ)*cosφ - spcλ*sinφ;
- final double scaleX = ANGULAR_SCALE / (ν*cosφ);
- final double scaleY = ANGULAR_SCALE / ρ;
- final double λt = λ + (cmsλ * scaleX); // The target
geographic coordinates
- final double φt = φ + (cmsφ * scaleY);
- if (dstPts != null) {
- dstPts[dstOff++] = λt;
- dstPts[dstOff++] = φt;
- if (isTarget3D) {
- double t1 = Δfmod * sin2φ; // A term in the
calculation of Δh
- double t2 = Δa;
- if (!isAbridged) {
- t1 /= νden; // = Δf⋅(b/a)⋅ν⋅sin²φ
- t2 *= νden; // = Δa⋅(a/ν)
- }
- dstPts[dstOff++] = h + spcλ*cosφ + tZ*sinφ + t1 - t2;
- }
- }
- if (!derivate) {
- return null;
- }
- /*
- * At this point the (Abridged) Molodensky transformation is finished.
- * Code below this point is only for computing the derivative, if
requested.
- * Note: variable names do not necessarily tell all the terms that
they contain.
- */
- final Matrix matrix = Matrices.createDiagonal(getTargetDimensions(),
getSourceDimensions());
- final double sinφcosφ = sinφ * cosφ;
- final double dν = eccentricitySquared*sinφcosφ / ν2den;
- final double dν3ρ = 3*dν * (1 - eccentricitySquared) / ν2den;
- // double dXdλ = spcλ;
- final double dYdλ = cmsλ * sinφ;
- final double dZdλ = cmsλ * cosφ;
- double dXdφ = dYdλ / cosφ;
- double dYdφ = -tZ*sinφ - cosφ*spcλ + t*(1 - 2*sin2φ);
- double dZdφ = tZ*cosφ - sinφ*spcλ;
- if (isAbridged) {
- /*
- * Δfmod = (a⋅Δf) + (f⋅Δa)
- * t = 2⋅Δfmod
- * dXdh = 0 so no need to set the matrix element.
- * dYdh = 0 so no need to set the matrix element.
- */
- dXdφ -= cmsλ * dν;
- dYdφ -= cmsφ * dν3ρ;
- dZdφ += t*cosφ*sinφ;
- } else {
- /*
- * Δfmod = b⋅Δf
- * t = Δf⋅[ν⋅(b/a) + ρ⋅(a/b)] (real ν and ρ, without +
h)
- * ν is actually ν + h
- * ρ is actually ρ + h
- */
- final double dρ = dν3ρ * νden * (semiMajor / ρ); // Reminder:
that ρ contains a h term.
- dXdφ -= dν * cmsλ * semiMajor / (νden*ν); // Reminder:
that ν contains a h term.
- dYdφ -= dρ * dZdφ - (Δfmod*(dν*2/(1 - eccentricitySquared) + (1 +
1/ν2den)*(dν - dρ))
- + Δa*(dν + 1)*eccentricitySquared) *
sinφcosφ / νden;
- if (isSource3D) {
- final double dXdh = cmsλ / ν;
- final double dYdh = -cmsφ / ρ;
- matrix.setElement(0, 2, -dXdh * scaleX);
- matrix.setElement(1, 2, +dYdh * scaleY);
- }
- final double t1 = Δfmod * (dν*sin2φ + 2*sinφcosφ);
- final double t2 = Δa * dν;
- dZdφ += t1/νden + t2*νden;
- }
- matrix.setElement(0, 0, 1 - spcλ * scaleX);
- matrix.setElement(1, 1, 1 + dYdφ * scaleY);
- matrix.setElement(0, 1, + dXdφ * scaleX);
- matrix.setElement(1, 0, - dYdλ * scaleY);
- if (isTarget3D) {
- matrix.setElement(2, 0, dZdλ);
- matrix.setElement(2, 1, dZdφ);
- }
- return matrix;
- }
-
- /**
- * {@inheritDoc}
- *
- * @return {@inheritDoc}
- */
- @Override
- protected int computeHashCode() {
- int code = super.computeHashCode() + Long.hashCode(
- Double.doubleToLongBits(Δa)
- + Double.doubleToLongBits(Δfmod)
- + 31 * (Double.doubleToLongBits(tX)
- + 31 * (Double.doubleToLongBits(tY)
- + 31 * (Double.doubleToLongBits(tZ)))));
- if (isAbridged) code = ~code;
- return code;
- }
-
- /**
- * Compares the specified object with this math transform for equality.
- *
- * @return {@inheritDoc}
- */
- @Override
- public boolean equals(final Object object, final ComparisonMode mode) {
- if (object == this) {
- // Slight optimization
- return true;
- }
- if (super.equals(object, mode)) {
- final MolodenskyFormula that = (MolodenskyFormula) object;
- return isSource3D == that.isSource3D
- && isTarget3D == that.isTarget3D
- && isAbridged == that.isAbridged
- && Numerics.epsilonEqual(tX, that.tX,
mode)
- && Numerics.epsilonEqual(tY, that.tY,
mode)
- && Numerics.epsilonEqual(tZ, that.tZ,
mode)
- && Numerics.epsilonEqual(Δa, that.Δa,
mode)
- && Numerics.epsilonEqual(Δfmod, that.Δfmod,
mode)
- && Numerics.epsilonEqual(semiMajor, that.semiMajor,
mode)
- && Numerics.epsilonEqual(eccentricitySquared,
that.eccentricitySquared, mode);
- // No need to compare the contextual parameters since this is
done by super-class.
- }
- return false;
- }
-}
diff --git
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
index f8039f182b..7ea6805057 100644
---
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
+++
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
@@ -18,18 +18,26 @@ package org.apache.sis.referencing.operation.transform;
import java.util.Arrays;
import javax.measure.Unit;
+import javax.measure.quantity.Length;
import org.opengis.util.FactoryException;
+import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransformFactory;
import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.datum.DatumShiftGrid;
+import org.apache.sis.referencing.datum.DefaultEllipsoid;
+import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.provider.Molodensky;
import org.apache.sis.referencing.operation.provider.AbridgedMolodensky;
import org.apache.sis.parameter.Parameters;
import org.apache.sis.measure.Units;
+import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.Debug;
+import org.apache.sis.util.internal.Numerics;
import static java.lang.Math.*;
@@ -81,14 +89,14 @@ import static java.lang.Math.*;
* @author Rueben Schulz (UBC)
* @author Martin Desruisseaux (IRD, Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.3
+ * @version 1.5
* @since 0.7
*/
-public class MolodenskyTransform extends MolodenskyFormula {
+public class MolodenskyTransform extends DatumShiftTransform {
/**
* Serial number for inter-operability with different versions.
*/
- private static final long serialVersionUID = 7206439437113286122L;
+ private static final long serialVersionUID = 1741638055238886288L;
/**
* Internal parameter descriptor, used only for debugging purpose.
@@ -99,12 +107,93 @@ public class MolodenskyTransform extends MolodenskyFormula
{
@Debug
private static ParameterDescriptorGroup DESCRIPTOR;
+ /**
+ * The value of 1/sin(1″) multiplied by the conversion factor from
arc-seconds to radians (π/180)/(60⋅60).
+ * This is the final multiplication factor for Δλ and Δφ.
+ */
+ static final double ANGULAR_SCALE = 1.00000000000391744;
+
+ /**
+ * {@code true} if the source coordinates have a height.
+ */
+ private final boolean isSource3D;
+
+ /**
+ * {@code true} if the target coordinates have a height.
+ */
+ private final boolean isTarget3D;
+
+ /**
+ * {@code true} for the abridged formula, or {@code false} for the
complete one.
+ */
+ private final boolean isAbridged;
+
+ /**
+ * Shift along the geocentric X axis (toward prime meridian)
+ * in units of the semi-major axis of the source ellipsoid.
+ *
+ * @see org.apache.sis.referencing.datum.BursaWolfParameters#tX
+ */
+ protected final double tX;
+
+ /**
+ * Shift along the geocentric Y axis (toward 90°E)
+ * in units of the semi-major axis of the source ellipsoid.
+ *
+ * @see org.apache.sis.referencing.datum.BursaWolfParameters#tY
+ */
+ protected final double tY;
+
+ /**
+ * Shift along the geocentric Z axis (toward north pole)
+ * in units of the semi-major axis of the source ellipsoid.
+ *
+ * @see org.apache.sis.referencing.datum.BursaWolfParameters#tZ
+ */
+ protected final double tZ;
+
+ /**
+ * Difference in the semi-major axes of the target and source ellipsoids:
{@code Δa = target a - source a}.
+ *
+ * @see DefaultEllipsoid#semiMajorAxisDifference(Ellipsoid)
+ */
+ private final double Δa;
+
+ /**
+ * Difference between the flattening of the target and source ellipsoids
(Δf), opportunistically modified
+ * with additional terms. The value depends on whether this Molodensky
transform is abridged or not:
+ *
+ * <ul>
+ * <li>For Molodensky, this field is set to (b⋅Δf).</li>
+ * <li>For Abridged Molodensky, this field is set to (a⋅Δf) +
(f⋅Δa).</li>
+ * </ul>
+ *
+ * where Δf = <var>target flattening</var> - <var>source flattening</var>.
+ */
+ final double Δfmod;
+
+ /**
+ * Semi-major axis length (<var>a</var>) of the source ellipsoid.
+ */
+ protected final double semiMajor;
+
+ /**
+ * The square of eccentricity of the source ellipsoid.
+ * This can be computed by ℯ² = (a²-b²)/a² where
+ * <var>a</var> is the <cite>semi-major</cite> axis length and
+ * <var>b</var> is the <cite>semi-minor</cite> axis length.
+ *
+ * @see DefaultEllipsoid#getEccentricitySquared()
+ */
+ protected final double eccentricitySquared;
+
/**
* The inverse of this Molodensky transform.
+ * Shall be considered final.
*
* @see #inverse()
*/
- private final MolodenskyTransform inverse;
+ private MolodenskyTransform inverse;
/**
* Creates a Molodensky transform from the specified parameters.
@@ -148,8 +237,9 @@ public class MolodenskyTransform extends MolodenskyFormula {
final double tX, final double tY, final double
tZ,
final boolean isAbridged)
{
- super(source, isSource3D, target, isTarget3D, tX, tY, tZ, null,
isAbridged,
- isAbridged ? AbridgedMolodensky.PARAMETERS :
Molodensky.PARAMETERS);
+ this(source, isSource3D, target, isTarget3D, tX, tY, tZ, null,
isAbridged,
+ isAbridged ? AbridgedMolodensky.PARAMETERS :
Molodensky.PARAMETERS);
+
if (!isSource3D && !isTarget3D) {
if (isAbridged && tX == 0 && tY == 0 && tZ == 0) {
inverse = new AbridgedMolodenskyTransform2D(this, source,
target);
@@ -169,24 +259,109 @@ public class MolodenskyTransform extends
MolodenskyFormula {
* @param target the target ellipsoid of the given {@code inverse}
transform.
*/
MolodenskyTransform(final MolodenskyTransform inverse, final Ellipsoid
source, final Ellipsoid target) {
- super(inverse, source, target, inverse.context.getDescriptor());
+ this(target, inverse.isTarget3D, source, inverse.isSource3D,
+ -inverse.tX, -inverse.tY, -inverse.tZ, inverse.grid,
+ inverse.isAbridged, inverse.context.getDescriptor());
+
this.inverse = inverse;
}
/**
- * Invoked by constructor and by {@link #getParameterValues()} for setting
all parameters other than axis lengths.
+ * Creates a Molodensky transform from the specified parameters.
+ * If a non-null {@code grid} is specified, it is caller's responsibility
to verify its validity.
*
- * @param pg where to set the parameters.
- * @param semiMinor ignored.
- * @param unit the unit of measurement to declare.
- * @param Δf the flattening difference to set, or NaN if this
method should fetch that value itself.
+ * @param source the source ellipsoid.
+ * @param isSource3D {@code true} if the source coordinates have a height.
+ * @param target the target ellipsoid.
+ * @param isTarget3D {@code true} if the target coordinates have a height.
+ * @param tX the geocentric <var>X</var> translation in same
units than the source ellipsoid axes.
+ * @param tY the geocentric <var>Y</var> translation in same
units than the source ellipsoid axes.
+ * @param tZ the geocentric <var>Z</var> translation in same
units than the source ellipsoid axes.
+ * @param grid interpolation grid in geocentric coordinates, or
{@code null} if none.
+ * @param isAbridged {@code true} for the abridged formula, or {@code
false} for the complete one.
+ * @param descriptor the contextual parameter descriptor.
*/
- @Override
- final void completeParameters(final Parameters pg, final double semiMinor,
final Unit<?> unit, double Δf) {
+ private MolodenskyTransform(final Ellipsoid source, final boolean
isSource3D,
+ final Ellipsoid target, final boolean
isTarget3D,
+ final double tX, final double tY, final double
tZ,
+ final DatumShiftGrid<?,?> grid, final boolean
isAbridged,
+ final ParameterDescriptorGroup descriptor)
+ {
+ super(descriptor, isSource3D, isTarget3D, grid);
+ ArgumentChecks.ensureNonNull("source", source);
+ ArgumentChecks.ensureNonNull("target", target);
+ final DefaultEllipsoid src = DefaultEllipsoid.castOrCopy(source);
+ this.isSource3D = isSource3D;
+ this.isTarget3D = isTarget3D;
+ this.isAbridged = isAbridged;
+ this.semiMajor = src.getSemiMajorAxis();
+ this.Δa = src.semiMajorAxisDifference(target);
+ this.tX = tX;
+ this.tY = tY;
+ this.tZ = tZ;
+
+ final double semiMinor = src.getSemiMinorAxis();
+ final double Δf = src.flatteningDifference(target);
+ eccentricitySquared = src.getEccentricitySquared();
+ Δfmod = isAbridged ? (semiMajor * Δf) + (semiMajor - semiMinor) * (Δa
/ semiMajor)
+ : (semiMinor * Δf);
+ /*
+ * Copy parameters to the ContextualParameter. Those parameters are
not used directly
+ * by MolodenskyTransform, but we need to store them in case the user
asks for them.
+ */
+ final Unit<Length> unit = src.getAxisUnit();
+ setContextParameters(semiMajor, semiMinor, unit, target);
+ completeParameters(context, unit, Δf);
+ /*
+ * Prepare two affine transforms to be executed before and after the
MolodenskyTransform:
+ *
+ * - A "normalization" transform for converting degrees to radians,
+ * - A "denormalization" transform for for converting radians to
degrees.
+ */
+ context.normalizeGeographicInputs(0);
+ context.denormalizeGeographicOutputs(0);
+ }
+
+ /**
+ * Sets parameter values in the given group for parameters other than axis
lengths.
+ * This method is invoked for both completing contextual parameters
({@code pg == context}) and
+ * for completing internal parameters ({@code pg != context}). When this
method is invoked, the
+ * following parameters are already set:
+ *
+ * <ul>
+ * <li>{@code "src_semi_major"}</li>
+ * <li>{@code "src_semi_minor"}</li>
+ * <li>{@code "tgt_semi_major"} (contextual parameters only)</li>
+ * <li>{@code "tgt_semi_minor"} (contextual parameters only)</li>
+ * </ul>
+ *
+ * This method sets the following parameters:
+ *
+ * <ul>
+ * <li><q>dim</q></li>
+ * <li><q>X-axis translation</q> (Molodensky only)</li>
+ * <li><q>Y-axis translation</q> (Molodensky only)</li>
+ * <li><q>Z-axis translation</q> (Molodensky only)</li>
+ * <li><q>Geocentric translation file</q> (Geocentric interpolations
only)</li>
+ * <li><q>Semi-major axis length difference</q> (Always for Molodensky,
internal WKT only for geocentric interpolations)</li>
+ * <li><q>Flattening difference</q> (Always for Molodensky, internal WKT
only for geocentric interpolations)</li>
+ * </ul>
+ *
+ * @param pg where to set the parameters.
+ * @param unit the unit of measurement to declare.
+ * @param Δf the flattening difference to set, or NaN if this method
should fetch that value itself.
+ */
+ private void completeParameters(final Parameters pg, final Unit<?> unit,
double Δf) {
if (Double.isNaN(Δf)) {
Δf = context.doubleValue(Molodensky.FLATTENING_DIFFERENCE);
}
- super.completeParameters(pg, semiMinor, unit, Δf);
+ /*
+ * Unconditionally set the "dim" parameters to the number of source
dimensions (do not check for consistency
+ * with the number of target dimensions) because source dimensions
determine the value of ellipsoidal heights,
+ * which may change the horizontal numerical values. By contrast, the
number of target dimensions does not have
+ * any impact on numerical values (it can just causes a drop of the
third value).
+ */
+ pg.getOrCreate(Molodensky.DIMENSION).setValue(getSourceDimensions());
pg.getOrCreate(Molodensky.TX) .setValue(tX, unit);
pg.getOrCreate(Molodensky.TY) .setValue(tY, unit);
pg.getOrCreate(Molodensky.TZ) .setValue(tZ, unit);
@@ -258,21 +433,156 @@ public class MolodenskyTransform extends
MolodenskyFormula {
return tX == 0 && tY == 0 && tZ == 0 && Δa == 0 && Δfmod == 0 &&
getSourceDimensions() == getTargetDimensions();
}
+ /**
+ * Gets the dimension of input points.
+ *
+ * @return the input dimension, which is 2 or 3.
+ */
+ @Override
+ public final int getSourceDimensions() {
+ return isSource3D ? 3 : 2;
+ }
+
+ /**
+ * Gets the dimension of output points.
+ *
+ * @return the output dimension, which is 2 or 3.
+ */
+ @Override
+ public final int getTargetDimensions() {
+ return isTarget3D ? 3 : 2;
+ }
+
/**
* Transforms the (λ,φ) or (λ,φ,<var>h</var>) coordinates between two
geographic CRS,
* and optionally returns the derivative at that location.
*
+ * @param dstPts the array into which the transformed coordinate is
returned, or {@code null}.
+ * @param srcOff the offset to the point to be transformed in the
source array.
+ * @param dstPts the array into which the transformed coordinates is
returned.
+ * @param dstOff the offset to the location of the transformed point
that is stored in the destination array.
+ * @param derivate {@code true} for computing the derivative, or {@code
false} if not needed.
* @return {@inheritDoc}
* @throws TransformException if the point cannot be transformed or
* if a problem occurred while calculating the derivative.
*/
@Override
public Matrix transform(final double[] srcPts, final int srcOff,
- final double[] dstPts, final int dstOff,
+ final double[] dstPts, int dstOff,
final boolean derivate) throws TransformException
{
- return transform(srcPts[srcOff], srcPts[srcOff+1], isSource3D ?
srcPts[srcOff+2] : 0,
- dstPts, dstOff, tX, tY, tZ, derivate);
+ final double λ = srcPts[srcOff];
+ final double φ = srcPts[srcOff+1];
+ final double h = isSource3D ? srcPts[srcOff+2] : 0;
+ /*
+ * Abridged Molodensky formulas from EPSG guidance note:
+ *
+ * ν = a / √(1 - ℯ²⋅sin²φ) : radius of
curvature in the prime vertical
+ * ρ = a⋅(1 – ℯ²) / (1 – ℯ²⋅sin²φ)^(3/2) : radius of
curvature in the meridian
+ * Δλ″ = (-tX⋅sinλ + tY⋅cosλ) / (ν⋅cosφ⋅sin1″)
+ * Δφ″ = (-tX⋅sinφ⋅cosλ - tY⋅sinφ⋅sinλ + tZ⋅cosφ + [a⋅Δf +
f⋅Δa]⋅sin(2φ)) / (ρ⋅sin1″)
+ * Δh = tX⋅cosφ⋅cosλ + tY⋅cosφ⋅sinλ + tZ⋅sinφ + (a⋅Δf +
f⋅Δa)⋅sin²φ - Δa
+ *
+ * we set:
+ *
+ * dfm = (a⋅Δf + f⋅Δa) in abridged case (b⋅Δf in non-abridged
case)
+ * sin(2φ) = 2⋅sin(φ)⋅cos(φ)
+ */
+ final double sinλ = sin(λ);
+ final double cosλ = cos(λ);
+ final double sinφ = sin(φ);
+ final double cosφ = cos(φ);
+ final double sin2φ = sinφ * sinφ;
+ final double ν2den = 1 - eccentricitySquared*sin2φ; //
Square of the denominator of ν
+ final double νden = sqrt(ν2den); //
Denominator of ν
+ final double ρden = ν2den * νden; //
Denominator of ρ
+ double ρ = semiMajor * (1 - eccentricitySquared) / ρden; //
Other notation: Rm = ρ
+ double ν = semiMajor / νden; //
Other notation: Rn = ν
+ double t = Δfmod * 2; //
A term in the calculation of Δφ
+ if (!isAbridged) {
+ ρ += h;
+ ν += h;
+ t = t*(0.5/νden + 0.5/ρden) // = Δf⋅[ν⋅(b/a) +
ρ⋅(a/b)] (without the +h in ν and ρ)
+ + Δa*eccentricitySquared/νden; // = Δa⋅[ℯ²⋅ν/a]
+ }
+ final double spcλ = tY*sinλ + tX*cosλ; // "spc"
stands for "sin plus cos"
+ final double cmsλ = tY*cosλ - tX*sinλ; // "cms"
stands for "cos minus sin"
+ final double cmsφ = (tZ + t*sinφ)*cosφ - spcλ*sinφ;
+ final double scaleX = ANGULAR_SCALE / (ν*cosφ);
+ final double scaleY = ANGULAR_SCALE / ρ;
+ final double λt = λ + (cmsλ * scaleX); // The target
geographic coordinates
+ final double φt = φ + (cmsφ * scaleY);
+ if (dstPts != null) {
+ dstPts[dstOff++] = λt;
+ dstPts[dstOff++] = φt;
+ if (isTarget3D) {
+ double t1 = Δfmod * sin2φ; // A term in the
calculation of Δh
+ double t2 = Δa;
+ if (!isAbridged) {
+ t1 /= νden; // = Δf⋅(b/a)⋅ν⋅sin²φ
+ t2 *= νden; // = Δa⋅(a/ν)
+ }
+ dstPts[dstOff++] = h + spcλ*cosφ + tZ*sinφ + t1 - t2;
+ }
+ }
+ if (!derivate) {
+ return null;
+ }
+ /*
+ * At this point the (Abridged) Molodensky transformation is finished.
+ * Code below this point is only for computing the derivative, if
requested.
+ * Note: variable names do not necessarily tell all the terms that
they contain.
+ */
+ final Matrix matrix = Matrices.createDiagonal(getTargetDimensions(),
getSourceDimensions());
+ final double sinφcosφ = sinφ * cosφ;
+ final double dν = eccentricitySquared*sinφcosφ / ν2den;
+ final double dν3ρ = 3*dν * (1 - eccentricitySquared) / ν2den;
+ // double dXdλ = spcλ;
+ final double dYdλ = cmsλ * sinφ;
+ final double dZdλ = cmsλ * cosφ;
+ double dXdφ = dYdλ / cosφ;
+ double dYdφ = -tZ*sinφ - cosφ*spcλ + t*(1 - 2*sin2φ);
+ double dZdφ = tZ*cosφ - sinφ*spcλ;
+ if (isAbridged) {
+ /*
+ * Δfmod = (a⋅Δf) + (f⋅Δa)
+ * t = 2⋅Δfmod
+ * dXdh = 0 so no need to set the matrix element.
+ * dYdh = 0 so no need to set the matrix element.
+ */
+ dXdφ -= cmsλ * dν;
+ dYdφ -= cmsφ * dν3ρ;
+ dZdφ += t*cosφ*sinφ;
+ } else {
+ /*
+ * Δfmod = b⋅Δf
+ * t = Δf⋅[ν⋅(b/a) + ρ⋅(a/b)] (real ν and ρ, without +
h)
+ * ν is actually ν + h
+ * ρ is actually ρ + h
+ */
+ final double dρ = dν3ρ * νden * (semiMajor / ρ); // Reminder:
that ρ contains a h term.
+ dXdφ -= dν * cmsλ * semiMajor / (νden*ν); // Reminder:
that ν contains a h term.
+ dYdφ -= dρ * dZdφ - (Δfmod*(dν*2/(1 - eccentricitySquared) + (1 +
1/ν2den)*(dν - dρ))
+ + Δa*(dν + 1)*eccentricitySquared) *
sinφcosφ / νden;
+ if (isSource3D) {
+ final double dXdh = cmsλ / ν;
+ final double dYdh = -cmsφ / ρ;
+ matrix.setElement(0, 2, -dXdh * scaleX);
+ matrix.setElement(1, 2, +dYdh * scaleY);
+ }
+ final double t1 = Δfmod * (dν*sin2φ + 2*sinφcosφ);
+ final double t2 = Δa * dν;
+ dZdφ += t1/νden + t2*νden;
+ }
+ matrix.setElement(0, 0, 1 - spcλ * scaleX);
+ matrix.setElement(1, 1, 1 + dYdφ * scaleY);
+ matrix.setElement(0, 1, + dXdφ * scaleX);
+ matrix.setElement(1, 0, - dYdλ * scaleY);
+ if (isTarget3D) {
+ matrix.setElement(2, 0, dZdλ);
+ matrix.setElement(2, 1, dZdφ);
+ }
+ return matrix;
}
/**
@@ -319,7 +629,7 @@ public class MolodenskyTransform extends MolodenskyFormula {
}
/*
* The code in the following loop is basically a copy-and-paste of the
code in the
- * MolodenskyFormula.transform(λ, φ, h, …) method, but without
derivative matrix
+ * MolodenskyTransform.transform(λ, φ, h, …) method, but without
derivative matrix
* computation and without support for interpolation of (tX,tY,tZ)
values in a grid.
*/
while (--numPts >= 0) {
@@ -407,4 +717,74 @@ public class MolodenskyTransform extends MolodenskyFormula
{
return DESCRIPTOR;
}
}
+
+ /**
+ * Returns a copy of internal parameter values of this transform.
+ * The returned group contains parameters for the source ellipsoid
semi-axis lengths
+ * and the differences between source and target ellipsoid parameters.
+ *
+ * <h4>Usage note</h4>
+ * This method is mostly for {@linkplain
org.apache.sis.io.wkt.Convention#INTERNAL debugging purposes}
+ * since the isolation of non-linear parameters in this class is highly
implementation dependent.
+ * Most GIS applications will instead be interested in the {@linkplain
#getContextualParameters()
+ * contextual parameters}.
+ *
+ * @return a copy of the internal parameter values for this transform.
+ */
+ @Debug
+ @Override
+ public ParameterValueGroup getParameterValues() {
+ final Unit<?> unit =
context.getOrCreate(Molodensky.SRC_SEMI_MAJOR).getUnit();
+ final double semiMinor =
context.getOrCreate(Molodensky.SRC_SEMI_MINOR).doubleValue(unit);
+ final Parameters pg =
Parameters.castOrWrap(getParameterDescriptors().createValue());
+ pg.getOrCreate(Molodensky.SRC_SEMI_MAJOR).setValue(semiMajor, unit);
+ pg.getOrCreate(Molodensky.SRC_SEMI_MINOR).setValue(semiMinor, unit);
+ completeParameters(pg, unit, Double.NaN);
+ return pg;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * @return {@inheritDoc}
+ */
+ @Override
+ protected int computeHashCode() {
+ int code = super.computeHashCode() + Long.hashCode(
+ Double.doubleToLongBits(Δa)
+ + Double.doubleToLongBits(Δfmod)
+ + 31 * (Double.doubleToLongBits(tX)
+ + 31 * (Double.doubleToLongBits(tY)
+ + 31 * (Double.doubleToLongBits(tZ)))));
+ if (isAbridged) code = ~code;
+ return code;
+ }
+
+ /**
+ * Compares the specified object with this math transform for equality.
+ *
+ * @return {@inheritDoc}
+ */
+ @Override
+ public boolean equals(final Object object, final ComparisonMode mode) {
+ if (object == this) {
+ // Slight optimization
+ return true;
+ }
+ if (super.equals(object, mode)) {
+ final MolodenskyTransform that = (MolodenskyTransform) object;
+ return isSource3D == that.isSource3D
+ && isTarget3D == that.isTarget3D
+ && isAbridged == that.isAbridged
+ && Numerics.epsilonEqual(tX, that.tX,
mode)
+ && Numerics.epsilonEqual(tY, that.tY,
mode)
+ && Numerics.epsilonEqual(tZ, that.tZ,
mode)
+ && Numerics.epsilonEqual(Δa, that.Δa,
mode)
+ && Numerics.epsilonEqual(Δfmod, that.Δfmod,
mode)
+ && Numerics.epsilonEqual(semiMajor, that.semiMajor,
mode)
+ && Numerics.epsilonEqual(eccentricitySquared,
that.eccentricitySquared, mode);
+ // No need to compare the contextual parameters since this is
done by super-class.
+ }
+ return false;
+ }
}
