Examples with DatumShiftMethod org.apache.sis.internal.referencing.provider.DatumShiftMethod used on opensource projects

Example 1 with DatumShiftMethod

use of org.apache.sis.internal.referencing.provider.DatumShiftMethod in project sis by apache.

the class CoordinateOperationFinder method createOperationStep.

/**
 * Creates an operation between two geodetic (geographic or geocentric) coordinate reference systems.
 * The default implementation can:
 *
 * <ul>
 *   <li>adjust axis order and orientation, for example converting from (<cite>North</cite>, <cite>West</cite>)
 *       axes to (<cite>East</cite>, <cite>North</cite>) axes,</li>
 *   <li>apply units conversion if needed,</li>
 *   <li>perform longitude rotation if needed,</li>
 *   <li>perform datum shift if {@linkplain BursaWolfParameters Bursa-Wolf parameters} are available
 *       for the area of interest.</li>
 * </ul>
 *
 * <p>This method returns only <em>one</em> step for a chain of concatenated operations (to be built by the caller).
 * But a list is returned because the same step may be implemented by different operation methods. Only one element
 * in the returned list should be selected (usually the first one).</p>
 *
 * @param  sourceCRS  input coordinate reference system.
 * @param  targetCRS  output coordinate reference system.
 * @return a coordinate operation from {@code sourceCRS} to {@code targetCRS}.
 * @throws FactoryException if the operation can not be constructed.
 */
@SuppressWarnings("null")
protected List<CoordinateOperation> createOperationStep(final GeodeticCRS sourceCRS, final GeodeticCRS targetCRS) throws FactoryException {
    final GeodeticDatum sourceDatum = sourceCRS.getDatum();
    final GeodeticDatum targetDatum = targetCRS.getDatum();
    Matrix datumShift = null;
    /*
         * If the prime meridian is not the same, we will concatenate a longitude rotation before or after datum shift
         * (that concatenation will be performed by the customized DefaultMathTransformFactory.Context created below).
         * Actually we do not know if the longitude rotation should be before or after datum shift. But this ambiguity
         * can usually be ignored because Bursa-Wolf parameters are always used with source and target prime meridians
         * set to Greenwich in EPSG dataset 8.9.  For safety, the SIS's DefaultGeodeticDatum class ensures that if the
         * prime meridians are not the same, then the target meridian must be Greenwich.
         */
    final DefaultMathTransformFactory.Context context = ReferencingUtilities.createTransformContext(sourceCRS, targetCRS, new MathTransformContext(sourceDatum, targetDatum));
    /*
         * If both CRS use the same datum and the same prime meridian, then the coordinate operation is only axis
         * swapping, unit conversion or change of coordinate system type (Ellipsoidal ↔ Cartesian ↔ Spherical).
         * Otherwise (if the datum are not the same), we will need to perform a scale, translation and rotation
         * in Cartesian space using the Bursa-Wolf parameters. If the user does not require the best accuracy,
         * then the Molodensky approximation may be used for avoiding the conversion step to geocentric CRS.
         */
    Identifier identifier;
    boolean isGeographicToGeocentric = false;
    final CoordinateSystem sourceCS = context.getSourceCS();
    final CoordinateSystem targetCS = context.getTargetCS();
    if (equalsIgnoreMetadata(sourceDatum, targetDatum)) {
        final boolean isGeocentricToGeographic;
        isGeographicToGeocentric = (sourceCS instanceof EllipsoidalCS && targetCS instanceof CartesianCS);
        isGeocentricToGeographic = (sourceCS instanceof CartesianCS && targetCS instanceof EllipsoidalCS);
        /*
             * Above booleans should never be true in same time. If it nevertheless happen (we are paranoiac;
             * maybe a lazy user implemented all interfaces in a single class), do not apply any geographic ↔
             * geocentric conversion. Instead do as if the coordinate system types were the same.
             */
        if (isGeocentricToGeographic ^ isGeographicToGeocentric) {
            identifier = GEOCENTRIC_CONVERSION;
        } else {
            identifier = AXIS_CHANGES;
        }
    } else {
        identifier = ELLIPSOID_CHANGE;
        if (sourceDatum instanceof DefaultGeodeticDatum) {
            datumShift = ((DefaultGeodeticDatum) sourceDatum).getPositionVectorTransformation(targetDatum, areaOfInterest);
            if (datumShift != null) {
                identifier = DATUM_SHIFT;
            }
        }
    }
    /*
         * Conceptually, all transformations below could done by first converting from the source coordinate
         * system to geocentric Cartesian coordinates (X,Y,Z), apply an affine transform represented by the
         * datum shift matrix, then convert from the (X′,Y′,Z′) coordinates to the target coordinate system.
         * However there is two exceptions to this path:
         *
         *   1) In the particular where both the source and target CS are ellipsoidal, we may use the
         *      Molodensky approximation as a shortcut (if the desired accuracy allows).
         *
         *   2) Even if we really go through the XYZ coordinates without Molodensky approximation, there is
         *      at least 9 different ways to name this operation depending on whether the source and target
         *      CRS are geocentric or geographic, 2- or 3-dimensional, whether there is a translation or not,
         *      the rotation sign, etc. We try to use the most specific name if we can find one, and fallback
         *      on an arbitrary name only in last resort.
         */
    final DefaultMathTransformFactory mtFactory = factorySIS.getDefaultMathTransformFactory();
    MathTransform before = null, after = null;
    ParameterValueGroup parameters;
    if (identifier == DATUM_SHIFT || identifier == ELLIPSOID_CHANGE) {
        /*
             * If the transform can be represented by a single coordinate operation, returns that operation.
             * Possible operations are:
             *
             *    - Position Vector transformation (in geocentric, geographic-2D or geographic-3D domains)
             *    - Geocentric translation         (in geocentric, geographic-2D or geographic-3D domains)
             *    - [Abridged] Molodensky          (as an approximation of geocentric translation)
             *    - Identity                       (if the desired accuracy is so large than we can skip datum shift)
             *
             * TODO: if both CS are ellipsoidal but with different number of dimensions, then we should use
             * an intermediate 3D geographic CRS in order to enable the use of Molodensky method if desired.
             */
        final DatumShiftMethod preferredMethod = DatumShiftMethod.forAccuracy(desiredAccuracy);
        parameters = GeocentricAffine.createParameters(sourceCS, targetCS, datumShift, preferredMethod);
        if (parameters == null) {
            /*
                 * Failed to select a coordinate operation. Maybe because the coordinate system types are not the same.
                 * Convert unconditionally to XYZ geocentric coordinates and apply the datum shift in that CS space.
                 *
                 * TODO: operation name should not be "Affine" if 'before' or 'after' transforms are not identity.
                 *       Reminder: the parameter group name here determines the OperationMethod later in this method.
                 */
            if (datumShift != null) {
                parameters = TensorParameters.WKT1.createValueGroup(properties(Constants.AFFINE), datumShift);
            } else {
                // Dimension of geocentric CRS.
                parameters = Affine.identity(3);
            }
            final CoordinateSystem normalized = CommonCRS.WGS84.geocentric().getCoordinateSystem();
            before = mtFactory.createCoordinateSystemChange(sourceCS, normalized, sourceDatum.getEllipsoid());
            after = mtFactory.createCoordinateSystemChange(normalized, targetCS, targetDatum.getEllipsoid());
            context.setSource(normalized);
            context.setTarget(normalized);
        }
    } else if (identifier == GEOCENTRIC_CONVERSION) {
        parameters = (isGeographicToGeocentric ? GeographicToGeocentric.PARAMETERS : GeocentricToGeographic.PARAMETERS).createValue();
    } else {
        /*
             * Coordinate system change (including change in the number of dimensions) without datum shift.
             */
        final int sourceDim = sourceCS.getDimension();
        final int targetDim = targetCS.getDimension();
        if (// sourceDim == 2 or 3.
        (sourceDim & ~1) == 2 && // abs(sourceDim - targetDim) == 1.
        (sourceDim ^ targetDim) == 1 && (sourceCS instanceof EllipsoidalCS) && (targetCS instanceof EllipsoidalCS)) {
            parameters = (sourceDim == 2 ? Geographic2Dto3D.PARAMETERS : Geographic3Dto2D.PARAMETERS).createValue();
        } else {
            /*
                 * TODO: instead than creating parameters for an identity operation, we should create the
                 *       CoordinateOperation directly from the MathTransform created by mtFactory below.
                 *       The intent if to get the correct OperationMethod, which should not be "Affine"
                 *       if there is a CS type change.
                 */
            parameters = Affine.identity(targetDim);
            /*
                 * createCoordinateSystemChange(…) needs the ellipsoid associated to the ellipsoidal coordinate system,
                 * if any. If none or both coordinate systems are ellipsoidal, then the ellipsoid will be ignored (see
                 * createCoordinateSystemChange(…) javadoc for the rational) so it does not matter which one we pick.
                 */
            before = mtFactory.createCoordinateSystemChange(sourceCS, targetCS, (sourceCS instanceof EllipsoidalCS ? sourceDatum : targetDatum).getEllipsoid());
            context.setSource(targetCS);
        }
    }
    /*
         * Transform between differents datums using Bursa Wolf parameters. The Bursa Wolf parameters are used
         * with "standard" geocentric CS, i.e. with X axis towards the prime meridian, Y axis towards East and
         * Z axis toward North, unless the Molodensky approximation is used. The following steps are applied:
         *
         *     source CRS                        →
         *     normalized CRS with source datum  →
         *     normalized CRS with target datum  →
         *     target CRS
         *
         * Those steps may be either explicit with the 'before' and 'after' transform, or implicit with the
         * Context parameter.
         */
    MathTransform transform = mtFactory.createParameterizedTransform(parameters, context);
    final OperationMethod method = mtFactory.getLastMethodUsed();
    if (before != null) {
        transform = mtFactory.createConcatenatedTransform(before, transform);
        if (after != null) {
            transform = mtFactory.createConcatenatedTransform(transform, after);
        }
    }
    return asList(createFromMathTransform(properties(identifier), sourceCRS, targetCRS, transform, method, parameters, null));
}