blob: 761fc018295b385ca57e50159bff0610253d122b [file] [log] [blame]
/*
* 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.projection;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
import java.util.Objects;
import java.io.Serializable;
import java.lang.reflect.Modifier;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptor;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransform2D;
import org.opengis.referencing.operation.SingleOperation;
import org.opengis.referencing.operation.OperationMethod;
import org.opengis.referencing.operation.TransformException;
import org.opengis.referencing.operation.MathTransformFactory;
import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.opengis.util.FactoryException;
import org.apache.sis.util.Debug;
import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.parameter.Parameters;
import org.apache.sis.parameter.ParameterBuilder;
import org.apache.sis.metadata.iso.citation.Citations;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.operation.transform.AbstractMathTransform;
import org.apache.sis.referencing.operation.transform.AbstractMathTransform2D;
import org.apache.sis.referencing.operation.transform.ContextualParameters;
import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory;
import org.apache.sis.referencing.operation.transform.MathTransformProvider;
import org.apache.sis.internal.referencing.provider.MapProjection;
import org.apache.sis.internal.referencing.CoordinateOperations;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.system.Modules;
import org.apache.sis.internal.system.Loggers;
import org.apache.sis.internal.util.Constants;
import org.apache.sis.internal.util.Numerics;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.logging.Logging;
import static java.lang.Math.*;
import org.opengis.referencing.ReferenceIdentifier;
/**
* Base class for conversion services between ellipsoidal and cartographic projections.
* This conversion works on a normalized spaces, where angles are expressed in radians and
* computations are performed for a sphere having a semi-major axis of 1. More specifically:
*
* <ul class="verbose">
* <li>On input, the {@link #transform(double[], int, double[], int, boolean) transform(…)} method
* expects (<var>longitude</var>, <var>latitude</var>) angles in <strong>radians</strong>,
* sometime pre-multiplied by other projection-specific factors (see point #3 below).
* Longitudes have the <cite>central meridian</cite> (λ₀) removed before the transform method is invoked.
* The conversion from degrees to radians and the longitude rotation are applied by the
* {@linkplain ContextualParameters#normalizeGeographicInputs normalization} affine transform.</li>
*
* <li>On output, the {@link #transform(double[],int,double[],int,boolean) transform(…)} method returns
* (<var>x</var>, <var>y</var>) values on a sphere or ellipse having a semi-major axis length (<var>a</var>) of 1,
* sometime divided by other projection-specific factors (see point #3 below).
* The multiplication by the scale factor (<var>k₀</var>) and the translation by false easting (FE) and false
* northing (FN) are applied by the {@linkplain ContextualParameters#getMatrix denormalization} affine transform.</li>
*
* <li>In addition to above-cited conversions, subclasses may opportunistically concatenate other linear operations
* (scales and translations). They do that by changing the normalization and denormalization matrices shown below.
* When such changes are applied, the {@code transform(…)} inputs are no longer angles in radians but some other
* derived values.</li>
* </ul>
*
* The normalization and denormalization steps are represented below by the matrices immediately on the left and right
* sides of {@code NormalizedProjection} respectively. Those matrices show only the basic parameters common to most projections.
* Some projections will put more elements in those matrices.
*
* <div class="horizontal-flow" style="align-items:center">
* <div>{@include ../transform/formulas.html#SwapAxes}</div>
* <div>→</div>
* <div>{@include ../transform/formulas.html#NormalizeGeographic}</div>
* <div>→</div>
* <div>{@code NormalizedProjection}</div>
* <div>→</div>
* <div>{@include ../transform/formulas.html#DenormalizeCartesian}</div>
* </div>
*
* <div class="note"><b>Note:</b>
* The first matrix on the left side is for {@linkplain org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes
* swapping axes} from (<var>latitude</var>, <var>longitude</var>) to (<var>longitude</var>, <var>latitude</var>) order.
* This matrix is shown here for completeness, but is not managed by this projection package. Axes swapping is managed
* at a {@linkplain org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory#createParameterizedTransform
* higher level}.</div>
*
* {@code NormalizedProjection} does not store the above cited parameters (central meridian, scale factor, <i>etc.</i>)
* on intent (except indirectly), in order to make clear that those parameters are not used by subclasses.
* The ability to recognize two {@code NormalizedProjection}s as {@linkplain #equals(Object, ComparisonMode) equivalent}
* without consideration for the scale factor (among other) allow more efficient concatenation in some cases
* (typically some combinations of inverse projection followed by a direct projection).
*
* <p>All angles (either fields, method parameters or return values) in this class and subclasses are
* in radians. This is the opposite of {@link Parameters} where all angles are in CRS-dependent units,
* typically decimal degrees.</p>
*
* <h2>Serialization</h2>
* Serialization of this class is appropriate for short-term storage or RMI use, but may not be compatible
* with future versions. For long term storage, WKT (Well Know Text) or XML are more appropriate.
*
* @author Martin Desruisseaux (MPO, IRD, Geomatys)
* @author André Gosselin (MPO)
* @author Rueben Schulz (UBC)
* @author Rémi Maréchal (Geomatys)
* @version 1.1
*
* @see ContextualParameters
* @see <a href="http://mathworld.wolfram.com/MapProjection.html">Map projections on MathWorld</a>
*
* @since 0.6
* @module
*/
public abstract class NormalizedProjection extends AbstractMathTransform2D implements Serializable {
/**
* For cross-version compatibility.
*/
private static final long serialVersionUID = -4010883312927645853L;
/**
* Number of input and output dimensions of all {@code NormalizedProjection}.
* We define this constant for clarity only; its value shall not be modified.
*
* @see #getSourceDimensions()
* @see #getTargetDimensions()
*/
static final int DIMENSION = 2;
/**
* Maximum difference allowed when comparing longitudes or latitudes in radians.
* The current value takes the system-wide angular tolerance value (equivalent to
* about 1 cm on Earth) converted to radians.
*
* <p>Some formulas use this tolerance value for testing sines or cosines of an angle.
* In the sine case, this is justified because sin(θ) ≅ θ when θ is small.
* Similar reasoning applies to cosine with cos(θ) ≅ θ + π/2 when θ is small.</p>
*
* <p>Some formulas may use this tolerance value as a <em>linear</em> tolerance on the unit sphere.
* This is okay because the arc length for an angular tolerance θ is r⋅θ, but in this class r=1.</p>
*/
static final double ANGULAR_TOLERANCE = Formulas.ANGULAR_TOLERANCE * (PI/180);
// Note: an alternative way to compute this value could be Formulas.LINEAR_TOLERANCE / AUTHALIC_RADIUS.
// But the later is only 0.07% lower than the current value.
/**
* Desired accuracy for the result of iterative computations, in radians.
* This constant defines the desired accuracy of methods like {@link ConformalProjection#φ(double)}.
*
* <p>The current value is 0.25 time the accuracy derived from {@link Formulas#LINEAR_TOLERANCE}.
* So if the linear tolerance is 1 cm, then the accuracy that we will seek for is 0.25 cm (about
* 4E-10 radians). The 0.25 factor is a safety margin for meeting the 1 cm accuracy.</p>
*/
static final double ITERATION_TOLERANCE = ANGULAR_TOLERANCE * 0.25;
/**
* Maximum number of iterations for iterative computations.
* The iterative methods used in subclasses should converge quickly (in 3 or 4 iterations)
* when used for a planet with an eccentricity similar to Earth. But we allow a high limit
* in case someone uses SIS for some planet with higher eccentricity.
*/
static final int MAXIMUM_ITERATIONS = Formulas.MAXIMUM_ITERATIONS;
/**
* The internal parameter descriptors. Keys are implementation classes. Values are parameter descriptor groups
* containing at least a parameter for the {@link #eccentricity} value, and optionally other internal parameter
* added by some subclasses.
*
* <p>Entries are created only when first needed. Those descriptors are usually never created since they are
* used only by {@link #getParameterDescriptors()}, which is itself invoked mostly for debugging purpose.</p>
*/
@Debug
private static final Map<Class<?>,ParameterDescriptorGroup> DESCRIPTORS = new HashMap<>();
/**
* The parameters used for creating this projection. They are used for formatting <cite>Well Known Text</cite> (WKT)
* and error messages. Subclasses shall not use the values defined in this object for computation purpose, except at
* construction time.
*
* @see #getContextualParameters()
*/
final ContextualParameters context;
/**
* Ellipsoid eccentricity, equals to <code>sqrt({@linkplain #eccentricitySquared})</code>.
* Value 0 means that the ellipsoid is spherical.
*/
protected final double eccentricity;
/**
* The square of eccentricity: ℯ² = (a²-b²)/a² where
* <var>ℯ</var> is the {@linkplain #eccentricity eccentricity},
* <var>a</var> is the <cite>semi-major</cite> axis length and
* <var>b</var> is the <cite>semi-minor</cite> axis length.
*/
protected final double eccentricitySquared;
/**
* The inverse of this map projection.
*
* <div class="note"><b>Note:</b>
* creation of this object is not deferred to the first call to the {@link #inverse()} method because this
* object is lightweight and typically needed soon anyway (may be as soon as {@code ConcatenatedTransform}
* construction time). In addition this field is part of serialization form in order to preserve the
* references graph.</div>
*/
private final MathTransform2D inverse;
/**
* Maps the parameters to be used for initializing {@link NormalizedProjection} and its
* {@linkplain ContextualParameters#getMatrix normalization / denormalization} matrices.
* This is an enumeration of parameters found in almost every map projections, but under different names.
* This enumeration allows {@code NormalizedProjection} subclasses to specify which parameter names, ranges
* and default values should be used by the
* {@linkplain NormalizedProjection#NormalizedProjection(OperationMethod, Parameters, Map) projection constructor}.
*
* <p>{@code NormalizedProjection} subclasses will typically provide values only for the following keys:
* {@link #CENTRAL_MERIDIAN}, {@link #SCALE_FACTOR}, {@link #FALSE_EASTING} and {@link #FALSE_NORTHING}.</p>
*
* @author Martin Desruisseaux (Geomatys)
* @version 0.6
*
* @see NormalizedProjection#NormalizedProjection(OperationMethod, Parameters, Map)
*
* @since 0.6
* @module
*/
protected enum ParameterRole {
/**
* Maps the <cite>semi-major axis length</cite> parameter (symbol: <var>a</var>).
* This value is used for computing {@link NormalizedProjection#eccentricity},
* and is also a multiplication factor for the denormalization matrix.
*
* <p>Unless specified otherwise, this is always mapped to a parameter named {@code "semi_major"}.
* {@code NormalizedProjection} subclasses typically do not need to provide a value for this key.</p>
*/
SEMI_MAJOR,
/**
* Maps the <cite>semi-minor axis length</cite> parameter (symbol: <var>b</var>).
* This value is used for computing {@link NormalizedProjection#eccentricity}.
*
* <p>Unless specified otherwise, this is always mapped to a parameter named {@code "semi_minor"}.
* {@code NormalizedProjection} subclasses typically do not need to provide a value for this key.</p>
*/
SEMI_MINOR,
/**
* Maps the parameter for the latitude where to compute the <cite>radius of conformal sphere</cite>
* (symbol: <var>R</var><sub>c</sub>). If this parameter is provided, then the radius of the conformal
* sphere at latitude φ will be used instead than the semi-major axis length in the denormalisation matrix.
* In other words, if provided then <var>a</var> is replaced by <var>R</var><sub>c</sub> below:
*
* <div style="text-align:center">{@include ../transform/formulas.html#DenormalizeCartesian}</div>
*
* <p>This enumeration shall be used <strong>only</strong> when the user requested explicitly spherical
* formulas, for example the <cite>"Mercator (Spherical)"</cite> projection (EPSG:1026), but the figure
* of the Earth may be an ellipsoid rather than a sphere. In the majority of cases, this enumeration should
* not be used.</p>
*/
LATITUDE_OF_CONFORMAL_SPHERE_RADIUS,
/**
* Maps the <cite>central meridian</cite> parameter (symbol: λ₀).
* This value is subtracted from the longitude values before the map projections.
*
* <p>Some common names for this parameter are:</p>
* <ul>
* <li>Longitude of origin</li>
* <li>Longitude of false origin</li>
* <li>Longitude of natural origin</li>
* <li>Spherical longitude of origin</li>
* <li>Longitude of projection centre</li>
* </ul>
*/
CENTRAL_MERIDIAN,
/**
* Maps the <cite>scale factor</cite> parameter (symbol: <var>k₀</var>).
* This is a multiplication factor for the (<var>x</var>,<var>y</var>) values obtained after map projections.
*
* <p>Some common names for this parameter are:</p>
* <ul>
* <li>Scale factor at natural origin</li>
* <li>Scale factor on initial line</li>
* <li>Scale factor on pseudo standard parallel</li>
* </ul>
*/
SCALE_FACTOR,
/**
* Maps the <cite>false easting</cite> parameter (symbol: <var>FE</var>).
* This is a translation term for the <var>x</var> values obtained after map projections.
*
* <p>Some common names for this parameter are:</p>
* <ul>
* <li>False easting</li>
* <li>Easting at false origin</li>
* <li>Easting at projection centre</li>
* </ul>
*/
FALSE_EASTING,
/**
* Maps the <cite>false westing</cite> parameter (symbol: <var>FW</var>).
* This is the same <var>x</var> translation than {@link #FALSE_EASTING}, but of opposite sign.
*
* <p>Actually, there is usually no parameter named "false westing" in a map projection.
* But some projections like <cite>"Lambert Conic Conformal (West Orientated)"</cite> are
* defined in such a way that their "false easting" parameter is effectively a "false westing".
* This enumeration value can be used for informing {@link NormalizedProjection} about that fact.</p>
*/
FALSE_WESTING,
/**
* Maps the <cite>false northing</cite> parameter (symbol: <var>FN</var>).
* This is a translation term for the <var>y</var> values obtained after map projections.
*
* <p>Some common names for this parameter are:</p>
* <ul>
* <li>False northing</li>
* <li>Northing at false origin</li>
* <li>Northing at projection centre</li>
* </ul>
*/
FALSE_NORTHING,
/**
* Maps the <cite>false southing</cite> parameter (symbol: <var>FS</var>).
* This is the same <var>y</var> translation than {@link #FALSE_NORTHING}, but of opposite sign.
*
* <p>Actually, there is usually no parameter named "false southing" in a map projection.
* But some projections like <cite>"Transverse Mercator (South Orientated)"</cite> are
* defined in such a way that their "false northing" parameter is effectively a "false southing".
* This enumeration value can be used for informing {@link NormalizedProjection} about that fact.</p>
*/
FALSE_SOUTHING
}
/**
* Constructs a new map projection from the supplied parameters.
* This constructor applies the following operations on the
* {@linkplain #getContextualParameters() contextual parameters}:
*
* <ul>
* <li>On the <b>normalization</b> matrix (to be applied before {@code this} transform):
* <ul>
* <li>{@linkplain ContextualParameters#normalizeGeographicInputs(double) Subtract}
* the <cite>central meridian</cite> value.</li>
* <li>Convert from degrees to radians.</li>
* </ul>
* </li>
* <li>On the <b>denormalization</b> matrix (to be applied after {@code this} transform):
* <ul>
* <li>{@linkplain MatrixSIS#convertAfter(int, Number, Number) Scale} by the <cite>semi-major</cite> axis length.</li>
* <li>If a scale factor is present (not all map projections have a scale factor), apply that scale.</li>
* <li>Translate by the <cite>false easting</cite> and <cite>false northing</cite> (after the scale).</li>
* </ul>
* </li>
* <li>On the <b>contextual parameters</b> (not the parameters of {@code this} transform):
* <ul>
* <li>Store the values for <cite>semi-major</cite> axis length, <cite>semi-minor</cite> axis length,
* <cite>scale factor</cite> (if present), <cite>central meridian</cite>,
* <cite>false easting</cite> and <cite>false northing</cite> values.</li>
* </ul>
* </li>
* </ul>
*
* In matrix form, this constructor creates the following matrices (subclasses are free to modify):
* <table class="sis">
* <caption>Initial matrix coefficients after construction</caption>
* <tr>
* <th>Normalization</th>
* <th class="sep">Denormalization</th>
* </tr>
* <tr>
* <td>{@include ../transform/formulas.html#NormalizeGeographic}</td>
* <td class="sep">{@include ../transform/formulas.html#DenormalizeCartesian}</td>
* </tr>
* </table>
*
* <h4>Which parameters are considered</h4>
* The {@code roles} map specifies which parameters to look for <cite>central meridian</cite>,
* <cite>scale factor</cite>, <cite>false easting</cite>, <cite>false northing</cite> and other values.
* All entries in the {@code roles} map are optional.
* All descriptors in the map shall comply to the following constraints:
*
* <ul>
* <li>Descriptors associated to {@link ParameterRole#SEMI_MAJOR}, {@link ParameterRole#SEMI_MINOR SEMI_MINOR},
* {@link ParameterRole#FALSE_EASTING FALSE_EASTING} and {@link ParameterRole#FALSE_NORTHING FALSE_NORTHING}
* shall have the same linear unit of measurement (usually metre).</li>
* <li>Descriptors associated to angular measures ({@link ParameterRole#CENTRAL_MERIDIAN} and
* {@link ParameterRole#LATITUDE_OF_CONFORMAL_SPHERE_RADIUS LATITUDE_OF_CONFORMAL_SPHERE_RADIUS})
* shall use degrees.</li>
* </ul>
*
* Note that users can still use units of their choice in the {@link Parameters} object given in argument to
* this constructor. But those values will be converted to the units of measurement specified by the parameter
* descriptors in the {@code roles} map, which must be the above-cited units.
*
* @param method description of the map projection parameters.
* @param parameters the parameters of the projection to be created.
* @param roles parameters to look for <cite>central meridian</cite>, <cite>scale factor</cite>,
* <cite>false easting</cite>, <cite>false northing</cite> and other values.
*/
protected NormalizedProjection(final OperationMethod method, final Parameters parameters,
final Map<ParameterRole, ? extends ParameterDescriptor<? extends Number>> roles)
{
this(new Initializer(method, parameters, roles, (byte) 0));
}
/**
* Creates a new normalized projection from the parameters computed by the given initializer.
*
* @param initializer the initializer for computing map projection internal parameters.
*/
NormalizedProjection(final Initializer initializer) {
context = initializer.context;
eccentricitySquared = initializer.eccentricitySquared.doubleValue();
eccentricity = sqrt(eccentricitySquared); // DoubleDouble.sqrt() does not make any difference here.
inverse = new Inverse(this);
}
/**
* Creates a new projection initialized to the values of the given one. This constructor may be invoked after
* we determined that the default implementation can be replaced by an other one, for example using spherical
* formulas instead than the ellipsoidal ones. This constructor allows to transfer all parameters to the new
* instance without recomputing them.
*/
NormalizedProjection(final NormalizedProjection other) {
context = other.context;
eccentricity = other.eccentricity;
eccentricitySquared = other.eccentricitySquared;
inverse = new Inverse(this);
}
/**
* Returns {@code true} if the projection specified by the given method has the given keyword or identifier.
* If non-null, the given identifier is presumed in the EPSG namespace and has precedence over the keyword.
*
* <div class="note"><b>Implementation note:</b>
* Since callers usually give a constant string for the {@code regex} argument, it would be more efficient to
* compile the {@link java.util.regex.Pattern} once for all. However the regular expression is used only as a
* fallback if the descriptor does not contain EPSG identifier, which should be rare. Usually, the regular
* expression will never be compiled.</div>
*
* @param method the user-specified projection method.
* @param regex the regular expression to use when using the operation name as the criterion.
* @param identifier the identifier to compare against the operation method name.
* @return {@code true} if the name of the given operation method contains the given keyword
* or has an EPSG identifier equals to the given identifier.
*/
static boolean identMatch(final OperationMethod method, final String regex, final String identifier) {
if (identifier != null) {
for (final ReferenceIdentifier id : method.getIdentifiers()) {
if (Constants.EPSG.equals(id.getCodeSpace())) {
return identifier.equals(id.getCode());
}
}
}
return method.getName().getCode().replace('_',' ').matches(regex);
}
/**
* Returns the sequence of <cite>normalization</cite> → {@code this} → <cite>denormalization</cite> transforms
* as a whole. The transform returned by this method expects (<var>longitude</var>, <var>latitude</var>)
* coordinates in <em>degrees</em> and returns (<var>x</var>,<var>y</var>) coordinates in <em>metres</em>.
* Conversion to other units and {@linkplain org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes
* changes in axis order} are <strong>not</strong> managed by the returned transform.
*
* <p>The default implementation is as below:</p>
* {@preformat java
* return getContextualParameters().completeTransform(factory, this);
* }
*
* Subclasses can override this method if they wish to use alternative implementations under some circumstances.
* For example many subclasses will replace {@code this} by a specialized implementation if they detect that the
* ellipsoid is actually spherical.
*
* @param factory the factory to use for creating the transform.
* @return the map projection from (λ,φ) to (<var>x</var>,<var>y</var>) coordinates.
* @throws FactoryException if an error occurred while creating a transform.
*
* @see ContextualParameters#completeTransform(MathTransformFactory, MathTransform)
*/
public MathTransform createMapProjection(final MathTransformFactory factory) throws FactoryException {
return context.completeTransform(factory, this);
}
/**
* If this map projection can not handle the parameters given by the user but an other projection could, delegates
* to the other projection. This method can be invoked by some {@link #createMapProjection(MathTransformFactory)}
* implementations when the other projection can be seen as a special case.
*
* <div class="note"><b>Example:</b>
* the {@link ObliqueStereographic} formulas do not work anymore when the latitude of origin is 90°N or 90°S,
* because some internal coefficients become infinite. However the {@link PolarStereographic} implementation
* is designed especially for those special cases. So the {@code ObliqueStereographic.createMapProjection(…)}
* method can redirect to {@code PolarStereographic.createMapProjection(…)} when it detects such cases.</div>
*
* It is caller's responsibility to choose an alternative method that can understand the parameters which were
* given to this original projection.
*
* @param factory the factory given to {@link #createMapProjection(MathTransformFactory)}.
* @param name the name of the alternative map projection to use.
* @return the alternative projection.
* @throws FactoryException if an error occurred while creating the alternative projection.
*
* @since 0.7
*/
final MathTransform delegate(final MathTransformFactory factory, final String name) throws FactoryException {
final OperationMethod method;
if (factory instanceof DefaultMathTransformFactory) {
method = ((DefaultMathTransformFactory) factory).getOperationMethod(name);
} else {
method = CoordinateOperations.getOperationMethod(factory.getAvailableMethods(SingleOperation.class), name);
}
if (method instanceof MathTransformProvider) {
return ((MathTransformProvider) method).createMathTransform(factory, context);
} else {
throw new FactoryException(Errors.format(Errors.Keys.UnsupportedImplementation_1,
(method != null ? method : factory).getClass()));
}
}
/**
* Returns the parameters used for creating the complete map projection. Those parameters describe a sequence of
* <cite>normalize</cite> → {@code this} → <cite>denormalize</cite> transforms, <strong>not</strong> including
* {@linkplain org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes axis swapping}.
* Those parameters are used for formatting <cite>Well Known Text</cite> (WKT) and error messages.
* Subclasses shall not use the values defined in the returned object for computation purpose,
* except at construction time.
*
* @return the parameters values for the sequence of <cite>normalize</cite> → {@code this} → <cite>denormalize</cite>
* transforms, or {@code null} if unspecified.
*/
@Override
protected final ContextualParameters getContextualParameters() {
return context;
}
/**
* Returns a copy of non-linear internal parameter values of this {@code NormalizedProjection}.
* The returned group contains at least the {@link #eccentricity} parameter value.
* Some subclasses add more non-linear parameters, but most of them do not because many parameters
* like the <cite>scale factor</cite> or the <cite>false easting/northing</cite> are handled by the
* {@linkplain ContextualParameters#getMatrix (de)normalization affine transforms} instead.
*
* <div class="note"><b>Note:</b>
* 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}.</div>
*
* @return a copy of the internal parameter values for this normalized projection.
*/
@Debug
@Override
public ParameterValueGroup getParameterValues() {
final ParameterValueGroup group = getParameterDescriptors().createValue();
group.parameter("eccentricity").setValue(eccentricity);
final String[] names = getInternalParameterNames();
final double[] values = getInternalParameterValues();
for (int i=0; i<names.length; i++) {
group.parameter(names[i]).setValue(values[i]);
}
return group;
}
/**
* Returns a description of the non-linear internal parameters of this {@code NormalizedProjection}.
* The returned group contains at least a descriptor for the {@link #eccentricity} parameter.
* Subclasses may add more parameters.
*
* <p>This method is for inspecting the parameter values of this non-linear kernel only,
* not for inspecting the {@linkplain #getContextualParameters() contextual parameters}.
* Inspecting the kernel parameter values is usually for debugging purpose only.</p>
*
* @return a description of the internal parameters.
*/
@Debug
@Override
public ParameterDescriptorGroup getParameterDescriptors() {
Class<?> type = getClass();
while (!Modifier.isPublic(type.getModifiers())) {
type = type.getSuperclass();
}
ParameterDescriptorGroup group;
synchronized (DESCRIPTORS) {
group = DESCRIPTORS.get(type);
if (group == null) {
final ParameterBuilder builder = new ParameterBuilder().setRequired(true);
if (type.getName().startsWith(Modules.CLASSNAME_PREFIX)) {
builder.setCodeSpace(Citations.SIS, Constants.SIS);
}
final String[] names = getInternalParameterNames();
final ParameterDescriptor<?>[] parameters = new ParameterDescriptor<?>[names.length + 1];
parameters[0] = MapProjection.ECCENTRICITY;
for (int i=1; i<parameters.length; i++) {
parameters[i] = builder.addName(names[i-1]).create(Double.class, null);
}
group = builder.addName(CharSequences.camelCaseToSentence(
type.getSimpleName()) + " (radians domain)").createGroup(1, 1, parameters);
DESCRIPTORS.put(type, group);
}
}
return group;
}
/**
* Returns the names of any additional internal parameters (other than {@link #eccentricity})
* that this projection has. The length of this array must be the same than the length of the
* {@link #getInternalParameterValues()} array, if the later is non-null.
*/
String[] getInternalParameterNames() {
return CharSequences.EMPTY_ARRAY;
}
/**
* Returns the values of any additional internal parameters (other than {@link #eccentricity}) that
* this projection has. Those values are also compared by {@link #equals(Object, ComparisonMode)}.
*/
double[] getInternalParameterValues() {
return null;
}
/**
* If the given scaled longitude θ=n⋅λ is outside the [−n⋅π … n⋅π] range, maybe shifts θ to that range.
* This method intentionally does <strong>not</strong> force θ to be inside that range in all cases.
* We avoid explicit wraparounds as much as possible (as opposed to implicit wraparounds performed by
* trigonometric functions) because they tend to introduce discontinuities. We perform wraparounds only
* when necessary for the problem of area spanning the anti-meridian (±180°).
*
* <div class="note"><b>Example:</b>
* a CRS for Alaska may have the central meridian at λ₀=−154° of longitude. If the point to project is
* at λ=177° of longitude, calculations will be performed with Δλ=331° while the correct value that we
* need to use is Δλ=−29°.</div>
*
* In order to avoid wraparound operations as much as possible, we test only the bound where anti-meridian
* problem may happen; no wraparound will be applied for the opposite bound. Furthermore we add or subtract
* 360° only once. Even if the point did many turns around the Earth, the 360° shift will still be applied
* at most once. The desire to apply the minimal amount of shifts is the reason why we do not use
* {@link Math#IEEEremainder(double, double)}.
*
* @param θ the scaled longitude value θ=n⋅λ where <var>n</var> is a projection-dependent factor.
* @param θ_bound minimal (if negative) or maximal (if positive) value of θ before to apply the shift.
* This is computed by <code>{@linkplain Initializer#boundOfScaledLongitude(double)
* Initializer.boundOfScaledLongitude}(n)</code>
* @return θ or shifted θ.
*
* @see Initializer#boundOfScaledLongitude(double)
* @see <a href="https://issues.apache.org/jira/browse/SIS-486">SIS-486</a>
*/
static double wraparoundScaledLongitude(double θ, final double θ_bound) {
if _bound < 0 ? θ < θ_bound : θ > θ_bound) {
θ -= 2_bound;
}
return θ;
}
/**
* Converts a single coordinate in {@code srcPts} at the given offset and stores the result
* in {@code dstPts} at the given offset. In addition, opportunistically computes the
* transform derivative if requested.
*
* <h4>Normalization</h4>
* The input coordinates are (<var>λ</var>,<var>φ</var>) (the variable names for <var>longitude</var> and
* <var>latitude</var> respectively) angles in radians, eventually pre-multiplied by projection-specific factors.
* Input coordinate shall have the <cite>central meridian</cite> removed from the longitude by the caller
* before this method is invoked. After this method is invoked, the caller will need to multiply the output
* coordinate by the global <cite>scale factor</cite>,
* apply the (<cite>false easting</cite>, <cite>false northing</cite>) offset
* and eventually other projection-specific factors.
* This means that projections that implement this method are performed on a sphere or ellipse
* having a semi-major axis length of 1.
*
* <div class="note"><b>Note 1:</b> it is generally not necessary to know the projection-specific additional
* factors applied by subclasses on the input and output values, because {@code NormalizedProjection} should
* never be used directly. {@code NormalizedProjection} instances are used only indirectly as a step in a
* concatenated transform that include the <cite>normalization</cite> and <cite>denormalization</cite>
* matrices documented in this class javadoc.</div>
*
* <div class="note"><b>Note 2:</b> in <a href="https://proj4.org/">Proj.4</a>, the same standardization,
* described above, is handled by {@code pj_fwd.c}, except for the projection-specific additional factors.</div>
*
* <h4>Argument checks</h4>
* The input longitude and latitude are usually (but not always) in the range [-π … π] and [-π/2 … π/2] respectively.
* However values outside those ranges are accepted on the assumption that most implementations use those values
* only in trigonometric functions like {@linkplain Math#sin(double) sine} and {@linkplain Math#cos(double) cosine}.
* If this assumption is not applicable to a particular subclass, then it is implementer responsibility to check
* the range.
*
* @param srcPts the array containing the source point coordinate, as (<var>longitude</var>, <var>latitude</var>)
* angles in <strong>radians</strong>.
* @param srcOff the offset of the single coordinate to be converted in the source array.
* @param dstPts the array into which the converted coordinate is returned (may be the same than {@code srcPts}).
* Coordinates will be expressed in a dimensionless unit, as a linear distance on a unit sphere or ellipse.
* @param dstOff the offset of the location of the converted coordinate that is stored in the destination array.
* @param derivate {@code true} for computing the derivative, or {@code false} if not needed.
* @return the matrix of the projection derivative at the given source position,
* or {@code null} if the {@code derivate} argument is {@code false}.
* @throws ProjectionException if the coordinate can not be converted.
*/
@Override
public abstract Matrix transform(double[] srcPts, int srcOff, double[] dstPts, int dstOff, boolean derivate)
throws ProjectionException;
/**
* Inverse converts the single coordinate in {@code srcPts} at the given offset and stores the result in
* {@code ptDst} at the given offset. The output coordinates are (<var>longitude</var>, <var>latitude</var>)
* angles in radians, usually (but not necessarily) in the range [-π … π] and [-π/2 … π/2] respectively.
*
* <h4>Normalization</h4>
* Input coordinate shall have the (<cite>false easting</cite>, <cite>false northing</cite>) removed
* by the caller and the result divided by the global <cite>scale factor</cite> before this method is invoked.
* After this method is invoked, the caller will need to add the <cite>central meridian</cite> to the longitude
* in the output coordinate. This means that projections that implement this method are performed on a sphere
* or ellipse having a semi-major axis of 1.
* Additional projection-specific factors may also need to be applied (see class javadoc).
*
* <div class="note"><b>Note:</b> in <a href="https://proj4.org/">Proj.4</a>, the same standardization,
* described above, is handled by {@code pj_inv.c}, except for the projection-specific additional factors.</div>
*
* @param srcPts the array containing the source point coordinate, as linear distance on a unit sphere or ellipse.
* @param srcOff the offset of the point to be converted in the source array.
* @param dstPts the array into which the converted point coordinate is returned (may be the same than {@code srcPts}).
* Coordinates will be (<var>longitude</var>, <var>latitude</var>) angles in <strong>radians</strong>.
* @param dstOff the offset of the location of the converted point that is stored in the destination array.
* @throws ProjectionException if the point can not be converted.
*/
protected abstract void inverseTransform(double[] srcPts, int srcOff, double[] dstPts, int dstOff)
throws ProjectionException;
/**
* Returns the inverse of this map projection.
* Subclasses do not need to override this method, as they should override
* {@link #inverseTransform(double[], int, double[], int) inverseTransform(…)} instead.
*
* @return the inverse of this map projection.
*/
@Override
public MathTransform2D inverse() {
return inverse;
}
/**
* Inverse of a normalized map projection.
*
* @author Martin Desruisseaux (Geomatys)
* @version 1.0
* @since 0.6
* @module
*/
private static final class Inverse extends AbstractMathTransform2D.Inverse implements Serializable {
/**
* For cross-version compatibility.
*/
private static final long serialVersionUID = 6014176098150309651L;
/**
* The enclosing transform.
*/
private final NormalizedProjection forward;
/**
* Default constructor.
*/
Inverse(final NormalizedProjection forward) {
this.forward = forward;
}
/**
* Returns the inverse of this math transform.
*/
@Override
public MathTransform2D inverse() {
return forward;
}
/**
* Inverse transforms the specified {@code srcPts} and stores the result in {@code dstPts}.
* If the derivative has been requested, then this method will delegate the derivative
* calculation to the enclosing class and inverts the resulting matrix.
*/
@Override
public Matrix transform(final double[] srcPts, final int srcOff,
double[] dstPts, int dstOff,
final boolean derivate) throws TransformException
{
if (!derivate) {
forward.inverseTransform(srcPts, srcOff, dstPts, dstOff);
return null;
} else {
if (dstPts == null) {
dstPts = new double[DIMENSION];
dstOff = 0;
}
forward.inverseTransform(srcPts, srcOff, dstPts, dstOff);
return Matrices.inverse(forward.transform(dstPts, dstOff, null, 0, true));
}
}
/**
* Inverse transforms an arbitrary amount of coordinates. This method optimizes the
* case where conversions can be applied by a loop with indices in increasing order.
*/
@Override
public void transform(final double[] srcPts, int srcOff,
final double[] dstPts, int dstOff, int numPts) throws TransformException
{
if (srcPts == dstPts && srcOff < dstOff) {
super.transform(srcPts, srcOff, dstPts, dstOff, numPts);
} else while (--numPts >= 0) {
forward.inverseTransform(srcPts, srcOff, dstPts, dstOff);
srcOff += DIMENSION;
dstOff += DIMENSION;
}
}
/**
* Concatenates or pre-concatenates in an optimized way this projection with the given transform, if possible.
* If no optimization is available, returns {@code null}.
*/
@Override
protected MathTransform tryConcatenate(final boolean applyOtherFirst, final MathTransform other,
final MathTransformFactory factory) throws FactoryException
{
final Matrix m = getMiddleMatrix(this, other, applyOtherFirst);
if (m != null) {
/*
* 'projectedSpace' values:
* - false if applyOtherFirst == false since we have (inverse projection) → (affine) → (projection).
* - true if applyOtherFirst == true since we have (projection) → (affine) → (inverse projection).
*/
return forward.tryConcatenate(applyOtherFirst, m, factory);
}
return super.tryConcatenate(applyOtherFirst, other, factory);
}
}
/**
* Concatenates or pre-concatenates in an optimized way this projection with the given transform, if possible.
* If transforms are concatenated in an (inverse projection) → (affine) → (projection) sequence where the
* (projection) and (inverse projection) steps are the {@linkplain #inverse() inverse} of each other,
* then in some particular case the sequence can be replaced by a single affine transform.
* If no such simplification is possible, this method returns {@code null}.
*
* @return the simplified (usually affine) transform, or {@code null} if no such optimization is available.
* @throws FactoryException if an error occurred while combining the transforms.
*
* @since 0.8
*/
@Override
protected MathTransform tryConcatenate(final boolean applyOtherFirst, final MathTransform other,
final MathTransformFactory factory) throws FactoryException
{
final Matrix m = getMiddleMatrix(this, other, applyOtherFirst);
if (m != null) {
/*
* 'projectedSpace' values:
* - false if applyOtherFirst == true since we have (inverse projection) → (affine) → (projection).
* - true if applyOtherFirst == false since we have (projection) → (affine) → (inverse projection).
*/
return tryConcatenate(!applyOtherFirst, m, factory);
}
return super.tryConcatenate(applyOtherFirst, other, factory);
}
/**
* Returns the concatenation of (inverse) → (affine) → (this) transforms, or {@code null} if none.
* The affine transform is applied in the geographic space with angular values in radians.
*
* <p>Above description is for the usual case where {@code projectedSpace} is {@code false}.
* If {@code true} (should be unusual), then the affine transform is rather applied in the
* projected space and the sequence to concatenate is (this) → (affine) → (inverse).</p>
*
* <p>Default implementation returns {@code null}. Subclasses should override if applicable.</p>
*
* @param projectedSpace {@code true} if affine transform is applied in projected instead than geographic space.
* @param affine the affine transform in the middle of (inverse) → (affine) → (this) transform sequence.
* @param factory the factory to use for creating new transform, or {@code null}.
* @return the optimized concatenation, or {@code null} if none.
*/
MathTransform tryConcatenate(boolean projectedSpace, Matrix affine, MathTransformFactory factory)
throws FactoryException
{
return null;
}
/**
* If a sequence of 3 transforms are (inverse projection) → (affine) → (projection) where
* the (projection) and (inverse projection) steps are the inverse of each other, returns
* the matrix of the affine transform step. Otherwise returns {@code null}. This method
* accepts also (projection) → (affine) → (inverse projection) sequence, but such sequences
* should be much more unusual.
*
* @param projection either {@link NormalizedProjection} or {@link Inverse}.
* @param other the arbitrary transforms to be concatenated with the given projection.
* @param applyOtherFirst whether {@code other} is concatenated before {@code projection} or the converse.
* @return the 3×3 matrix of the affine transform step, or {@code null} if none.
*/
private static Matrix getMiddleMatrix(final AbstractMathTransform projection, final MathTransform other,
final boolean applyOtherFirst)
{
final List<MathTransform> steps = MathTransforms.getSteps(other);
if (steps.size() == 2) try {
final int oi = applyOtherFirst ? 0 : 1;
if (projection.equals(steps.get(oi).inverse(), ComparisonMode.IGNORE_METADATA)) {
final Matrix m = MathTransforms.getMatrix(steps.get(oi ^ 1));
if (Matrices.isAffine(m) && m.getNumRow() == DIMENSION+1 && m.getNumCol() == DIMENSION+1) {
return m;
}
}
} catch (NoninvertibleTransformException e) {
Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION),
(projection instanceof NormalizedProjection) ? NormalizedProjection.class : projection.getClass(),
"tryConcatenate", e);
}
return null;
}
/**
* Computes a hash code value for this {@code NormalizedProjection}.
*
* @return the hash code value.
*/
@Override
protected int computeHashCode() {
long c = Double.doubleToLongBits(eccentricity);
final double[] parameters = getInternalParameterValues();
if (parameters != null) {
for (int i=0; i<parameters.length; i++) {
c = c*31 + Double.doubleToLongBits(parameters[i]);
}
}
return super.computeHashCode() ^ Long.hashCode(c);
}
/**
* Compares the given object with this transform for equivalence. The default implementation checks if
* {@code object} is an instance of the same class than {@code this}, then compares the eccentricity.
*
* <p>If this method returns {@code true}, then for any given identical source position, the two compared map
* projections shall compute the same target position. Many of the {@linkplain #getContextualParameters()
* contextual parameters} used for creating the map projections are irrelevant and do not need to be known.
* Those projection parameters will be compared only if the comparison mode is {@link ComparisonMode#STRICT}
* or {@link ComparisonMode#BY_CONTRACT BY_CONTRACT}.</p>
*
* <div class="note"><b>Example:</b>
* a {@linkplain Mercator Mercator} projection can be created in the 2SP case with a <cite>standard parallel</cite>
* value of 60°. The same projection can also be created in the 1SP case with a <cite>scale factor</cite> of 0.5.
* Nevertheless those two map projections applied on a sphere gives identical results. Considering them as
* equivalent allows the referencing module to transform coordinates between those two projections more efficiently.
* </div>
*
* @param object the object to compare with this map projection for equivalence.
* @param mode the strictness level of the comparison. Default to {@link ComparisonMode#STRICT}.
* @return {@code true} if the given object is equivalent to this map projection.
*/
@Override
@SuppressWarnings("fallthrough")
public boolean equals(final Object object, final ComparisonMode mode) {
if (object == this) {
return true;
}
if (!super.equals(object, mode)) {
return false;
}
final NormalizedProjection that = (NormalizedProjection) object;
switch (mode) {
case STRICT:
case BY_CONTRACT: {
if (!Objects.equals(context, that.context)) {
return false;
}
// Fall through for comparing the eccentricity.
}
case IGNORE_METADATA: {
/*
* There is no need to compare both 'eccentricity' and 'eccentricitySquared' since the former
* is computed from the later. We are better to compare 'eccentricitySquared' since it is the
* original value from which the other value is derived.
*/
if (!Numerics.equals(eccentricitySquared, that.eccentricitySquared)) {
return false;
}
break;
}
default: {
/*
* We want to compare the eccentricity with a tolerance threshold corresponding approximately
* to an error of 1 cm on Earth. The eccentricity for an ellipsoid of semi-major axis a=1 is:
*
* ℯ² = 1 - b²
*
* If we add a slight ε error to the semi-minor axis length (where ε will be our linear tolerance
* threshold), we get:
*
* (ℯ + ε′)² = 1 - (b + ε)² ≈ 1 - (b² + 2⋅b⋅ε) assuming ε ≪ b
*
* Replacing 1 - b² by ℯ²:
*
* ℯ² + 2⋅ℯ⋅ε′ ≈ ℯ² - 2⋅b⋅ε
*
* After a few rearrangements:
*
* ε′ ≈ ε⋅(ℯ - 1/ℯ)
*
* Note that ε′ is negative for ℯ < 1 so we actually need to compute ε⋅(1/ℯ - ℯ) instead.
* The result is less than 2E-8 for the eccentricity of the Earth.
*/
final double e = max(eccentricity, that.eccentricity);
if (!Numerics.epsilonEqual(eccentricity, that.eccentricity, ANGULAR_TOLERANCE * (1/e - e))) {
assert (mode != ComparisonMode.DEBUG) : Numerics.messageForDifference(
"eccentricity", eccentricity, that.eccentricity);
return false;
}
break;
}
}
final double[] parameters = getInternalParameterValues();
if (parameters != null) {
/*
* super.equals(…) guarantees that the two objects are of the same class.
* So in SIS implementation, this implies that the arrays have the same length.
*/
final double[] others = that.getInternalParameterValues();
assert others.length == parameters.length;
for (int i=0; i<parameters.length; i++) {
if (!Numerics.epsilonEqual(parameters[i], others[i], mode)) {
assert (mode != ComparisonMode.DEBUG) : Numerics.messageForDifference(
getInternalParameterNames()[i], parameters[i], others[i]);
return false;
}
}
}
return true;
}
}