    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_testing/src/org/deegree/crs/transformations/TransformationFactory.java $
    002    /*----------------------------------------------------------------------------
    003     This file is part of deegree, http://deegree.org/
    004     Copyright (C) 2001-2009 by:
    005     Department of Geography, University of Bonn
    006     and
    007     lat/lon GmbH
    008
    009     This library is free software; you can redistribute it and/or modify it under
    010     the terms of the GNU Lesser General Public License as published by the Free
    011     Software Foundation; either version 2.1 of the License, or (at your option)
    012     any later version.
    013     This library is distributed in the hope that it will be useful, but WITHOUT
    014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
    015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
    016     details.
    017     You should have received a copy of the GNU Lesser General Public License
    018     along with this library; if not, write to the Free Software Foundation, Inc.,
    019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
    020
    021     Contact information:
    022
    023     lat/lon GmbH
    024     Aennchenstr. 19, 53177 Bonn
    025     Germany
    026     http://lat-lon.de/
    027
    028     Department of Geography, University of Bonn
    029     Prof. Dr. Klaus Greve
    030     Postfach 1147, 53001 Bonn
    031     Germany
    032     http://www.geographie.uni-bonn.de/deegree/
    033
    034     e-mail: info@deegree.org
    035     ----------------------------------------------------------------------------*/
    036    package org.deegree.crs.transformations;
    037
    038    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
    039
    040    import java.util.Arrays;
    041
    042    import javax.vecmath.GMatrix;
    043    import javax.vecmath.Matrix3d;
    044    import javax.vecmath.Matrix4d;
    045
    046    import org.deegree.crs.Identifiable;
    047    import org.deegree.crs.components.Axis;
    048    import org.deegree.crs.components.Ellipsoid;
    049    import org.deegree.crs.components.GeodeticDatum;
    050    import org.deegree.crs.components.Unit;
    051    import org.deegree.crs.coordinatesystems.CompoundCRS;
    052    import org.deegree.crs.coordinatesystems.CoordinateSystem;
    053    import org.deegree.crs.coordinatesystems.GeocentricCRS;
    054    import org.deegree.crs.coordinatesystems.GeographicCRS;
    055    import org.deegree.crs.coordinatesystems.ProjectedCRS;
    056    import org.deegree.crs.exceptions.TransformationException;
    057    import org.deegree.crs.transformations.coordinate.CRSTransformation;
    058    import org.deegree.crs.transformations.coordinate.ConcatenatedTransform;
    059    import org.deegree.crs.transformations.coordinate.DirectTransform;
    060    import org.deegree.crs.transformations.coordinate.GeocentricTransform;
    061    import org.deegree.crs.transformations.coordinate.MatrixTransform;
    062    import org.deegree.crs.transformations.coordinate.ProjectionTransform;
    063    import org.deegree.crs.transformations.helmert.Helmert;
    064    import org.deegree.crs.transformations.polynomial.PolynomialTransformation;
    065    import org.deegree.crs.utilities.Matrix;
    066    import org.deegree.framework.log.ILogger;
    067    import org.deegree.framework.log.LoggerFactory;
    068    import org.deegree.i18n.Messages;
    069
    070    /**
    071     * The <code>TransformationFactory</code> class is the central access point for all transformations between different
    072     * crs's.
    073     * <p>
    074     * It creates a transformation chain for two given CoordinateSystems by considering their type. For example the
    075     * Transformation chain from EPSG:31466 ( a projected crs with underlying geographic crs epsg:4314 using the DHDN datum
    076     * and the TransverseMercator Projection) to EPSG:28992 (another projected crs with underlying geographic crs epsg:4289
    077     * using the 'new Amersfoort Datum' and the StereographicAzimuthal Projection) would result in following Transformation
    078     * Chain:
    079     * <ol>
    080     * <li>Inverse projection - thus getting the coordinates in lat/lon for geographic crs epsg:4314</li>
    081     * <li>Geodetic transformation - thus getting x-y-z coordinates for geographic crs epsg:4314</li>
    082     * <li>WGS84 transformation -thus getting the x-y-z coordinates for the WGS84 datum</li>
    083     * <li>Inverse WGS84 transformation -thus getting the x-y-z coordinates for the geodetic from epsg:4289</li>
    084     * <li>Inverse geodetic - thus getting the lat/lon for epsg:4289</li>
    085     * <li>projection - getting the coordinates (in meters) for epsg:28992</li>
    086     * </ol>
    087     *
    088     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
    089     *
    090     * @author last edited by: $Author: rbezema $
    091     *
    092     * @version $Revision: 19653 $, $Date: 2009-09-15 14:56:30 +0200 (Di, 15. Sep 2009) $
    093     *
    094     */
    095    public class TransformationFactory {
    096        private static ILogger LOG = LoggerFactory.getLogger( TransformationFactory.class );
    097
    098        /**
    099         * The default coordinate transformation factory. Will be constructed only when first needed.
    100         */
    101        private static TransformationFactory DEFAULT_INSTANCE = null;
    102
    103        private TransformationFactory() {
    104            // nottin
    105        }
    106
    107        /**
    108         * @return the default coordinate transformation factory.
    109         */
    110        public static synchronized TransformationFactory getInstance() {
    111            if ( DEFAULT_INSTANCE == null ) {
    112                DEFAULT_INSTANCE = new TransformationFactory();
    113            }
    114            return DEFAULT_INSTANCE;
    115        }
    116
    117        /**
    118         * Creates a transformation between two coordinate systems. This method will examine the coordinate systems in order
    119         * to construct a transformation between them.
    120         *
    121         * @param sourceCRS
    122         *            Input coordinate system.
    123         * @param targetCRS
    124         *            Output coordinate system.
    125         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    126         * @throws TransformationException
    127         * @throws TransformationException
    128         *             if no transformation path has been found.
    129         * @throws IllegalArgumentException
    130         *             if the sourceCRS or targetCRS are <code>null</code>.
    131         *
    132         */
    133        public Transformation createFromCoordinateSystems( final CoordinateSystem sourceCRS,
    134                                                           final CoordinateSystem targetCRS )
    135                                throws TransformationException, IllegalArgumentException {
    136            if ( sourceCRS == null ) {
    137                throw new IllegalArgumentException( "The source CRS may not be null" );
    138            }
    139            if ( targetCRS == null ) {
    140                throw new IllegalArgumentException( "The target CRS may not be null" );
    141            }
    142            if ( ( sourceCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS
    143                   && sourceCRS.getType() != CoordinateSystem.COMPOUND_CRS
    144                   && sourceCRS.getType() != CoordinateSystem.PROJECTED_CRS && sourceCRS.getType() != CoordinateSystem.GEOCENTRIC_CRS )
    145                 || ( targetCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS
    146                      && targetCRS.getType() != CoordinateSystem.COMPOUND_CRS
    147                      && targetCRS.getType() != CoordinateSystem.PROJECTED_CRS && targetCRS.getType() != CoordinateSystem.GEOCENTRIC_CRS ) ) {
    148                throw new TransformationException( sourceCRS, targetCRS,
    149                                                   "Either the target crs type or the source crs type was unknown" );
    150            }
    151
    152            if ( sourceCRS.equals( targetCRS ) ) {
    153                LOG.logDebug( "Source crs and target crs are equal, no transformation needed (returning identity matrix)." );
    154                final Matrix matrix = new Matrix( sourceCRS.getDimension() + 1 );
    155                matrix.setIdentity();
    156                return createMatrixTransform( sourceCRS, targetCRS, matrix );
    157            }
    158
    159            Transformation result = null;
    160            // check if the source crs has an alternative transformation for the given target, if so use it
    161            if ( sourceCRS.hasDirectTransformation( targetCRS ) ) {
    162                PolynomialTransformation direct = sourceCRS.getDirectTransformation( targetCRS );
    163                LOG.logDebug( "Using direct (polynomial) transformation instead of a helmert transformation: "
    164                              + direct.getImplementationName() );
    165                result = new DirectTransform( direct, sourceCRS, new Identifiable( new String[] { direct.getIdentifier()
    166                                                                                                  + "-CRSTransformation" } ) );
    167            } else {
    168                if ( sourceCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
    169                    /**
    170                     * Geographic --> Geographic, Projected, Geocentric or Compound
    171                     */
    172                    final GeographicCRS source = (GeographicCRS) sourceCRS;
    173                    if ( targetCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
    174                        result = createTransformation( source, (ProjectedCRS) targetCRS );
    175                    } else if ( targetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
    176                        result = createTransformation( source, (GeographicCRS) targetCRS );
    177                    } else if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
    178                        result = createTransformation( source, (GeocentricCRS) targetCRS );
    179                    } else if ( targetCRS.getType() == CoordinateSystem.COMPOUND_CRS ) {
    180                        CompoundCRS target = (CompoundCRS) targetCRS;
    181                        CompoundCRS sTmp = new CompoundCRS( target.getHeightAxis(), source, target.getDefaultHeight(),
    182                                                            new Identifiable( new String[] { source.getIdentifier()
    183                                                                                             + "_compound" } ) );
    184                        result = createTransformation( sTmp, target );
    185                    }
    186                } else if ( sourceCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
    187                    /**
    188                     * Projected --> Projected, Geographic, Geocentric or Compound
    189                     */
    190                    final ProjectedCRS source = (ProjectedCRS) sourceCRS;
    191                    if ( targetCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
    192                        result = createTransformation( source, (ProjectedCRS) targetCRS );
    193                    } else if ( targetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
    194                        result = createTransformation( source, (GeographicCRS) targetCRS );
    195                    } else if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
    196                        result = createTransformation( source, (GeocentricCRS) targetCRS );
    197                    } else if ( targetCRS.getType() == CoordinateSystem.COMPOUND_CRS ) {
    198                        CompoundCRS target = (CompoundCRS) targetCRS;
    199                        CompoundCRS sTmp = new CompoundCRS( target.getHeightAxis(), source, target.getDefaultHeight(),
    200                                                            new Identifiable( new String[] { source.getIdentifier()
    201                                                                                             + "_compound" } ) );
    202                        result = createTransformation( sTmp, target );
    203                    }
    204                } else if ( sourceCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
    205                    /**
    206                     * Geocentric --> Projected, Geographic, Geocentric or Compound
    207                     */
    208                    final GeocentricCRS source = (GeocentricCRS) sourceCRS;
    209                    if ( targetCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
    210                        result = createTransformation( source, (ProjectedCRS) targetCRS );
    211                    } else if ( targetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
    212                        result = createTransformation( source, (GeographicCRS) targetCRS );
    213                    } else if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
    214                        result = createTransformation( source, (GeocentricCRS) targetCRS );
    215                    } else if ( targetCRS.getType() == CoordinateSystem.COMPOUND_CRS ) {
    216                        CompoundCRS target = (CompoundCRS) targetCRS;
    217                        CompoundCRS sTmp = new CompoundCRS( target.getHeightAxis(), source, target.getDefaultHeight(),
    218                                                            new Identifiable( new String[] { source.getIdentifier()
    219                                                                                             + "_compound" } ) );
    220                        result = createTransformation( sTmp, target );
    221                    }
    222                } else if ( sourceCRS.getType() == CoordinateSystem.COMPOUND_CRS ) {
    223                    /**
    224                     * Compound --> Projected, Geographic, Geocentric or Compound
    225                     */
    226                    final CompoundCRS source = (CompoundCRS) sourceCRS;
    227                    CompoundCRS target = null;
    228                    if ( targetCRS.getType() != CoordinateSystem.COMPOUND_CRS ) {
    229                        target = new CompoundCRS(
    230                                                  source.getHeightAxis(),
    231                                                  targetCRS,
    232                                                  source.getDefaultHeight(),
    233                                                  new Identifiable(
    234                                                                    new String[] { targetCRS.getIdentifier() + "_compound" } ) );
    235                    } else {
    236                        target = (CompoundCRS) targetCRS;
    237                    }
    238                    result = createTransformation( source, target );
    239                }
    240            }
    241            if ( result == null ) {
    242                LOG.logDebug( "The resulting transformation was null, returning an identity matrix." );
    243                final Matrix matrix = new Matrix( sourceCRS.getDimension() + 1 );
    244                matrix.setIdentity();
    245                result = createMatrixTransform( sourceCRS, targetCRS, matrix );
    246            } else {
    247                LOG.logDebug( "Concatenating the result, with the conversion matrices." );
    248                result = concatenate(
    249                                      createMatrixTransform( sourceCRS, sourceCRS, toStandardizedValues( sourceCRS, false ) ),
    250                                      result, createMatrixTransform( targetCRS, targetCRS, toStandardizedValues( targetCRS,
    251                                                                                                                 true ) ) );
    252            }
    253            if ( LOG.getLevel() == ILogger.LOG_DEBUG ) {
    254                StringBuilder output = new StringBuilder( "The resulting transformation chain: \n" );
    255                output = result.getTransformationPath( output );
    256                LOG.logDebug( output.toString() );
    257
    258                if ( result instanceof MatrixTransform ) {
    259                    LOG.logDebug( "Resulting matrix transform:\n" + ( (MatrixTransform) result ).getMatrix() );
    260                }
    261
    262            }
    263            return result;
    264
    265        }
    266
    267        /**
    268         * Creates a matrix, with which incoming values will be transformed to a standardized form. This means, to radians
    269         * and meters.
    270         *
    271         * @param sourceCRS
    272         *            to create the matrix for.
    273         * @param invert
    274         *            the values. Using the inverted scale, i.e. going from standard to crs specific.
    275         * @return the standardized matrix.
    276         * @throws TransformationException
    277         *             if the unit of one of the axis could not be transformed to one of the base units.
    278         */
    279        private Matrix toStandardizedValues( CoordinateSystem sourceCRS, boolean invert )
    280                                throws TransformationException {
    281            final int dim = sourceCRS.getDimension();
    282            Matrix result = null;
    283            Axis[] allAxis = sourceCRS.getAxis();
    284            for ( int i = 0; i < allAxis.length; ++i ) {
    285                Axis targetAxis = allAxis[i];
    286                if ( targetAxis != null ) {
    287                    Unit targetUnit = targetAxis.getUnits();
    288                    if ( !( Unit.RADIAN.equals( targetUnit ) || Unit.METRE.equals( targetUnit ) ) ) {
    289                        if ( !( targetUnit.canConvert( Unit.RADIAN ) || targetUnit.canConvert( Unit.METRE ) ) ) {
    290                            throw new TransformationException(
    291                                                               Messages.getMessage(
    292                                                                                    "CRS_TRANSFORMATION_NO_APLLICABLE_UNIT",
    293                                                                                    targetUnit ) );
    294                        }
    295                        // lazy instantiation
    296                        if ( result == null ) {
    297                            result = new Matrix( dim + 1 );
    298                            result.setIdentity();
    299                        }
    300                        result.setElement( i, i, invert ? 1d / targetUnit.getScale() : targetUnit.getScale() );
    301                    }
    302                }
    303            }
    304            return result;
    305        }
    306
    307        /**
    308         * @param sourceCRS
    309         * @param targetCRS
    310         * @return the transformation chain or <code>null</code> if the transformation operation is the identity.
    311         * @throws TransformationException
    312         */
    313        private Transformation createTransformation( GeocentricCRS sourceCRS, GeographicCRS targetCRS )
    314                                throws TransformationException {
    315            final Transformation result = createTransformation( targetCRS, sourceCRS );
    316            if ( result != null ) {
    317                result.inverse();
    318            }
    319            return result;
    320        }
    321
    322        /**
    323         * @param sourceCRS
    324         * @param targetCRS
    325         * @return the transformation chain or <code>null</code> if the transformation operation is the identity.
    326         * @throws TransformationException
    327         */
    328        private Transformation createTransformation( GeocentricCRS sourceCRS, ProjectedCRS targetCRS )
    329                                throws TransformationException {
    330            final Transformation result = createTransformation( targetCRS, sourceCRS );
    331            if ( result != null ) {
    332                result.inverse();
    333            }
    334            return result;
    335        }
    336
    337        /**
    338         * This method is valid for all transformations which use a compound crs, because the extra heightvalues need to be
    339         * considered throughout the transformation.
    340         *
    341         * @param sourceCRS
    342         * @param targetCRS
    343         * @return the transformation chain or <code>null</code> if the transformation operation is the identity.
    344         * @throws TransformationException
    345         */
    346        private Transformation createTransformation( CompoundCRS sourceCRS, CompoundCRS targetCRS )
    347                                throws TransformationException {
    348            if ( sourceCRS.getUnderlyingCRS().equals( targetCRS.getUnderlyingCRS() ) ) {
    349                return null;
    350            }
    351            LOG.logDebug( "Creating compound( " + sourceCRS.getUnderlyingCRS().getIdentifier()
    352                          + ") ->compound transformation( " + targetCRS.getUnderlyingCRS().getIdentifier()
    353                          + "): from (source): " + sourceCRS.getIdentifier() + " to(target): " + targetCRS.getIdentifier() );
    354            final int sourceType = sourceCRS.getUnderlyingCRS().getType();
    355            final int targetType = targetCRS.getUnderlyingCRS().getType();
    356            Transformation result = null;
    357            // basic check for simple (invert) projections
    358            if ( sourceType == CoordinateSystem.PROJECTED_CRS && targetType == CoordinateSystem.GEOGRAPHIC_CRS ) {
    359                if ( ( ( (ProjectedCRS) sourceCRS.getUnderlyingCRS() ).getGeographicCRS() ).equals( targetCRS.getUnderlyingCRS() ) ) {
    360                    result = new ProjectionTransform( (ProjectedCRS) sourceCRS.getUnderlyingCRS() );
    361                    result.inverse();
    362                }
    363            }
    364            if ( sourceType == CoordinateSystem.GEOGRAPHIC_CRS && targetType == CoordinateSystem.PROJECTED_CRS ) {
    365                if ( ( ( (ProjectedCRS) targetCRS.getUnderlyingCRS() ).getGeographicCRS() ).equals( sourceCRS.getUnderlyingCRS() ) ) {
    366                    result = new ProjectionTransform( (ProjectedCRS) targetCRS.getUnderlyingCRS() );
    367                }
    368            }
    369            if ( result == null ) {
    370                GeocentricCRS sourceGeocentric = null;
    371                if ( sourceType == CoordinateSystem.GEOCENTRIC_CRS ) {
    372                    sourceGeocentric = (GeocentricCRS) sourceCRS.getUnderlyingCRS();
    373                } else {
    374                    sourceGeocentric = new GeocentricCRS( sourceCRS.getGeodeticDatum(), "tmp_" + sourceCRS.getIdentifier()
    375                                                                                        + "_geocentric",
    376                                                          sourceCRS.getName() + "_Geocentric" );
    377                }
    378                GeocentricCRS targetGeocentric = null;
    379                if ( targetType == CoordinateSystem.GEOCENTRIC_CRS ) {
    380                    targetGeocentric = (GeocentricCRS) targetCRS.getUnderlyingCRS();
    381                } else {
    382                    targetGeocentric = new GeocentricCRS( targetCRS.getGeodeticDatum(), "tmp_" + targetCRS.getIdentifier()
    383                                                                                        + "_geocentric",
    384                                                          targetCRS.getName() + "_Geocentric" );
    385                }
    386                Transformation helmertTransformation = createTransformation( sourceGeocentric, targetGeocentric );
    387
    388                Transformation sourceTransformationChain = null;
    389                Transformation targetTransformationChain = null;
    390                GeographicCRS sourceGeographic = null;
    391                GeographicCRS targetGeographic = null;
    392                switch ( sourceType ) {
    393                case CoordinateSystem.GEOCENTRIC_CRS:
    394                    break;
    395                case CoordinateSystem.PROJECTED_CRS:
    396                    sourceTransformationChain = new ProjectionTransform( (ProjectedCRS) sourceCRS.getUnderlyingCRS() );
    397                    sourceTransformationChain.inverse();
    398                    sourceGeographic = ( (ProjectedCRS) sourceCRS.getUnderlyingCRS() ).getGeographicCRS();
    399                case CoordinateSystem.GEOGRAPHIC_CRS:
    400                    if ( sourceGeographic == null ) {
    401                        sourceGeographic = (GeographicCRS) sourceCRS.getUnderlyingCRS();
    402                    }
    403                    /*
    404                     * Only create a geocentric transform if the helmert transformation != null, e.g. the datums and
    405                     * ellipsoids are not equal.
    406                     */
    407                    if ( helmertTransformation != null ) {
    408                        // create a 2d->3d mapping.
    409                        final CRSTransformation axisAligned = createMatrixTransform( sourceGeographic, sourceGeocentric,
    410                                                                                     swapAxis( sourceGeographic,
    411                                                                                               GeographicCRS.WGS84 ) );
    412                        if ( LOG.isDebug() ) {
    413                            StringBuilder sb = new StringBuilder(
    414                                                                  "Resulting axis alignment between source geographic and source geocentric is:" );
    415                            if ( axisAligned == null ) {
    416                                sb.append( " not necessary" );
    417                            } else {
    418                                sb.append( "\n" ).append( ( (MatrixTransform) axisAligned ).getMatrix() );
    419                            }
    420                            LOG.logDebug( sb.toString() );
    421                        }
    422                        final CRSTransformation geoCentricTransform = new GeocentricTransform( sourceCRS, sourceGeocentric );
    423                        // concatenate the possible projection with the axis alignment and the geocentric transform.
    424                        sourceTransformationChain = concatenate( sourceTransformationChain, axisAligned,
    425                                                                 geoCentricTransform );
    426                    }
    427                    break;
    428                }
    429                switch ( targetType ) {
    430                case CoordinateSystem.GEOCENTRIC_CRS:
    431                    break;
    432                case CoordinateSystem.PROJECTED_CRS:
    433                    targetTransformationChain = new ProjectionTransform( (ProjectedCRS) targetCRS.getUnderlyingCRS() );
    434                    targetGeographic = ( (ProjectedCRS) targetCRS.getUnderlyingCRS() ).getGeographicCRS();
    435                case CoordinateSystem.GEOGRAPHIC_CRS:
    436                    if ( targetGeographic == null ) {
    437                        targetGeographic = (GeographicCRS) targetCRS.getUnderlyingCRS();
    438                    }
    439                    /*
    440                     * Only create a geocentric transform if the helmert transformation != null, e.g. the datums and
    441                     * ellipsoids are not equal.
    442                     */
    443                    if ( helmertTransformation != null ) {
    444                        // create a 2d->3d mapping.
    445                        final CRSTransformation axisAligned = createMatrixTransform( targetGeocentric, targetGeographic,
    446                                                                                     swapAxis( GeographicCRS.WGS84,
    447                                                                                               targetGeographic ) );
    448                        final CRSTransformation geoCentricTransform = new GeocentricTransform( targetCRS, targetGeocentric );
    449                        geoCentricTransform.inverse();
    450                        // concatenate the possible projection with the axis alignment and the geocentric transform.
    451                        targetTransformationChain = concatenate( geoCentricTransform, axisAligned,
    452                                                                 targetTransformationChain );
    453                    }
    454                    break;
    455                }
    456                result = concatenate( sourceTransformationChain, helmertTransformation, targetTransformationChain );
    457            }
    458            return result;
    459        }
    460
    461        /**
    462         * Creates a transformation between two geographic coordinate systems. This method is automatically invoked by
    463         * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust axis
    464         * order and orientation (e.g. transforming from <code>(NORTH,WEST)</code> to <code>(EAST,NORTH)</code>), performs
    465         * units conversion and apply Bursa Wolf transformation if needed.
    466         *
    467         * @param sourceCRS
    468         *            Input coordinate system.
    469         * @param targetCRS
    470         *            Output coordinate system.
    471         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    472         * @throws TransformationException
    473         *             if no transformation path has been found.
    474         */
    475        private Transformation createTransformation( final GeographicCRS sourceCRS, final GeographicCRS targetCRS )
    476                                throws TransformationException {
    477            final GeodeticDatum sourceDatum = sourceCRS.getGeodeticDatum();
    478            final GeodeticDatum targetDatum = targetCRS.getGeodeticDatum();
    479            if ( sourceDatum.equals( targetDatum ) ) {
    480                LOG.logDebug( "The datums of geographic (source): " + sourceCRS.getIdentifier()
    481                              + " equals geographic (target): " + targetCRS.getIdentifier() + " returning null" );
    482                return null;
    483            }
    484            LOG.logDebug( "Creating geographic ->geographic transformation: from (source): " + sourceCRS.getIdentifier()
    485                          + " to(target): " + targetCRS.getIdentifier() );
    486            final Ellipsoid sourceEllipsoid = sourceDatum.getEllipsoid();
    487            final Ellipsoid targetEllipsoid = targetDatum.getEllipsoid();
    488            /*
    489             * If the two geographic coordinate systems use different ellipsoid, convert from the source to target ellipsoid
    490             * through the geocentric coordinate system.
    491             */
    492            if ( sourceEllipsoid != null
    493                 && !sourceEllipsoid.equals( targetEllipsoid )
    494                 && ( ( sourceDatum.getWGS84Conversion() != null && sourceDatum.getWGS84Conversion().hasValues() ) || ( targetDatum.getWGS84Conversion() != null && targetDatum.getWGS84Conversion().hasValues() ) ) ) {
    495
    496                Transformation step1 = null;
    497                CRSTransformation step2 = null;
    498                Transformation step3 = null;
    499                // use the WGS84 Geocentric transform if no toWGS84 parameters are given and the datums ellipsoid is
    500                // actually a sphere.
    501                final GeocentricCRS sourceGCS = ( sourceEllipsoid.isSphere() && ( sourceDatum.getWGS84Conversion() == null || !sourceDatum.getWGS84Conversion().hasValues() ) ) ? GeocentricCRS.WGS84
    502                                                                                                                                                                               : new GeocentricCRS(
    503                                                                                                                                                                                                    sourceDatum,
    504                                                                                                                                                                                                    sourceCRS.getIdentifier(),
    505                                                                                                                                                                                                    sourceCRS.getName()
    506                                                                                                                                                                                                                            + "_Geocentric" );
    507                final GeocentricCRS targetGCS = ( targetEllipsoid.isSphere() && ( targetDatum.getWGS84Conversion() == null || !targetDatum.getWGS84Conversion().hasValues() ) ) ? GeocentricCRS.WGS84
    508                                                                                                                                                                               : new GeocentricCRS(
    509                                                                                                                                                                                                    targetDatum,
    510                                                                                                                                                                                                    targetCRS.getIdentifier(),
    511                                                                                                                                                                                                    targetCRS.getName()
    512                                                                                                                                                                                                                            + "_Geocentric" );
    513
    514                // geographic->geocentric
    515                step1 = createTransformation( sourceCRS, sourceGCS );
    516                // helmert->inv_helmert
    517                step2 = createTransformation( sourceGCS, targetGCS );
    518                // geocentric->geographic
    519                step3 = createTransformation( targetCRS, targetGCS );
    520
    521                if ( step3 != null ) {
    522                    step3.inverse();// call inverseTransform from step 3.
    523                }
    524
    525                return concatenate( step1, step2, step3 );
    526            }
    527
    528            /*
    529             * Swap axis order, and rotate the longitude coordinate if prime meridians are different.
    530             */
    531            final Matrix matrix = swapAndRotateGeoAxis( sourceCRS, targetCRS );
    532            CRSTransformation result = createMatrixTransform( sourceCRS, targetCRS, matrix );
    533            if ( LOG.isDebug() ) {
    534                StringBuilder sb = new StringBuilder(
    535                                                      "Resulting axis alignment between source geographic and target geographic is:" );
    536                if ( result == null ) {
    537                    sb.append( " not necessary" );
    538                } else {
    539                    sb.append( "\n" ).append( ( (MatrixTransform) result ).getMatrix() );
    540                }
    541                LOG.logDebug( sb.toString() );
    542            }
    543            return result;
    544        }
    545
    546        /**
    547         * Creates a transformation between a geographic and a projected coordinate systems. This method is automatically
    548         * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
    549         *
    550         * @param sourceCRS
    551         *            Input coordinate system.
    552         * @param targetCRS
    553         *            Output coordinate system.
    554         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    555         * @throws TransformationException
    556         *             if no transformation path has been found.
    557         */
    558        private Transformation createTransformation( final GeographicCRS sourceCRS, final ProjectedCRS targetCRS )
    559                                throws TransformationException {
    560
    561            LOG.logDebug( "Creating geographic->projected transformation: from (source): " + sourceCRS.getIdentifier()
    562                          + " to(target): " + targetCRS.getIdentifier() );
    563            final GeographicCRS stepGeoCS = targetCRS.getGeographicCRS();
    564
    565            final Transformation geo2geo = createTransformation( sourceCRS, stepGeoCS );
    566            final CRSTransformation swap = createMatrixTransform( stepGeoCS, targetCRS, swapAxis( stepGeoCS, targetCRS ) );
    567            if ( LOG.isDebug() ) {
    568                StringBuilder sb = new StringBuilder(
    569                                                      "Resulting axis alignment between target geographic and target projected is:" );
    570                if ( swap == null ) {
    571                    sb.append( " not necessary" );
    572                } else {
    573                    sb.append( "\n" ).append( ( (MatrixTransform) swap ).getMatrix() );
    574                }
    575                LOG.logDebug( sb.toString() );
    576            }
    577            final CRSTransformation projection = new ProjectionTransform( targetCRS );
    578            return concatenate( geo2geo, swap, projection );
    579        }
    580
    581        /**
    582         * Creates a transformation between a geographic and a geocentric coordinate systems. Since the source coordinate
    583         * systems doesn't have a vertical axis, height above the ellipsoid is assumed equals to zero everywhere. This
    584         * method is automatically invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
    585         *
    586         * @param sourceCRS
    587         *            Input geographic coordinate system.
    588         * @param targetCRS
    589         *            Output coordinate system.
    590         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    591         * @throws TransformationException
    592         *             if no transformation path has been found.
    593         */
    594        private Transformation createTransformation( final GeographicCRS sourceCRS, final GeocentricCRS targetCRS )
    595                                throws TransformationException {
    596            LOG.logDebug( "Creating geographic -> geocentric (helmert) transformation: from (source): "
    597                          + sourceCRS.getIdentifier() + " to(target): " + targetCRS.getIdentifier() );
    598            /*
    599             * if ( !PrimeMeridian.GREENWICH.equals( targetCRS.getGeodeticDatum().getPrimeMeridian() ) ) { throw new
    600             * TransformationException( "The rotation from " +
    601             * targetCRS.getGeodeticDatum().getPrimeMeridian().getIdentifier() + " to Greenwich prime meridian is not yet
    602             * implemented" ); }
    603             */
    604            GeocentricCRS sourceGeocentric = new GeocentricCRS( sourceCRS.getGeodeticDatum(), "tmp_"
    605                                                                                              + sourceCRS.getIdentifier()
    606                                                                                              + "_geocentric",
    607                                                                sourceCRS.getName() + "_geocentric" );
    608            CRSTransformation helmertTransformation = createTransformation( sourceGeocentric, targetCRS );
    609            // if no helmert transformation is needed, the targetCRS equals the source-geocentric.
    610            if ( helmertTransformation == null ) {
    611                sourceGeocentric = targetCRS;
    612            }
    613            final CRSTransformation axisAlign = createMatrixTransform(
    614                                                                       sourceCRS,
    615                                                                       sourceGeocentric,
    616                                                                       swapAndRotateGeoAxis( sourceCRS, GeographicCRS.WGS84 ) );
    617            if ( LOG.isDebug() ) {
    618                StringBuilder sb = new StringBuilder(
    619                                                      "Resulting axis alignment between source geographic and target geocentric is:" );
    620                if ( axisAlign == null ) {
    621                    sb.append( " not necessary" );
    622                } else {
    623                    sb.append( "\n" ).append( ( (MatrixTransform) axisAlign ).getMatrix() );
    624                }
    625                LOG.logDebug( sb.toString() );
    626            }
    627            final CRSTransformation geocentric = new GeocentricTransform( sourceCRS, sourceGeocentric );
    628            return concatenate( axisAlign, geocentric, helmertTransformation );
    629        }
    630
    631        /**
    632         * Creates a transformation between two projected coordinate systems. This method is automatically invoked by
    633         * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust axis
    634         * order and orientation. It also performs units conversion if it is the only extra change needed. Otherwise, it
    635         * performs three steps:
    636         *
    637         * <ol>
    638         * <li>Unproject <code>sourceCRS</code>.</li>
    639         * <li>Transform from <code>sourceCRS.geographicCS</code> to <code>
    640         * targetCRS.geographicCS</code>.</li>
    641         * <li>Project <code>targetCRS</code>.</li>
    642         * </ol>
    643         *
    644         * @param sourceCRS
    645         *            Input coordinate system.
    646         * @param targetCRS
    647         *            Output coordinate system.
    648         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    649         * @throws TransformationException
    650         *             if no transformation path has been found.
    651         */
    652        private Transformation createTransformation( final ProjectedCRS sourceCRS, final ProjectedCRS targetCRS )
    653                                throws TransformationException {
    654            LOG.logDebug( "Creating projected -> projected transformation: from (source): " + sourceCRS.getIdentifier()
    655                          + " to(target): " + targetCRS.getIdentifier() );
    656            if ( sourceCRS.getProjection().equals( targetCRS.getProjection() ) ) {
    657                return null;
    658            }
    659            final GeographicCRS sourceGeo = sourceCRS.getGeographicCRS();
    660            final GeographicCRS targetGeo = targetCRS.getGeographicCRS();
    661            final Transformation inverseProjection = createTransformation( sourceCRS, sourceGeo );
    662            final Transformation geo2geo = createTransformation( sourceGeo, targetGeo );
    663            final Transformation projection = createTransformation( targetGeo, targetCRS );
    664            return concatenate( inverseProjection, geo2geo, projection );
    665        }
    666
    667        /**
    668         * Creates a transformation between a projected and a geocentric coordinate systems. This method is automatically
    669         * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. This method doesn't need to be
    670         * public since its decomposition in two step should be general enough.
    671         *
    672         * @param sourceCRS
    673         *            Input projected coordinate system.
    674         * @param targetCRS
    675         *            Output coordinate system.
    676         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    677         * @throws TransformationException
    678         *             if no transformation path has been found.
    679         */
    680        private Transformation createTransformation( final ProjectedCRS sourceCRS, final GeocentricCRS targetCRS )
    681                                throws TransformationException {
    682            LOG.logDebug( "Creating projected -> geocentric transformation: from (source): " + sourceCRS.getIdentifier()
    683                          + " to(target): " + targetCRS.getIdentifier() );
    684            final GeographicCRS sourceGCS = sourceCRS.getGeographicCRS();
    685
    686            final Transformation inverseProjection = createTransformation( sourceCRS, sourceGCS );
    687            final Transformation geocentric = createTransformation( sourceGCS, targetCRS );
    688            return concatenate( inverseProjection, geocentric );
    689        }
    690
    691        /**
    692         * Creates a transformation between a projected and a geographic coordinate systems. This method is automatically
    693         * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation
    694         * returns <code>{@link #createTransformation(GeographicCRS, ProjectedCRS)} createTransformation}(targetCRS,
    695         * sourceCRS) inverse)</code>.
    696         *
    697         * @param sourceCRS
    698         *            Input coordinate system.
    699         * @param targetCRS
    700         *            Output coordinate system.
    701         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code> or <code>null</code> if
    702         *         {@link ProjectedCRS#getGeographicCRS()}.equals targetCRS.
    703         * @throws TransformationException
    704         *             if no transformation path has been found.
    705         */
    706        private Transformation createTransformation( final ProjectedCRS sourceCRS, final GeographicCRS targetCRS )
    707                                throws TransformationException {
    708            LOG.logDebug( "Creating projected->geographic transformation: from (source): " + sourceCRS.getIdentifier()
    709                          + " to(target): " + targetCRS.getIdentifier() );
    710            Transformation result = createTransformation( targetCRS, sourceCRS );
    711            if ( result != null ) {
    712                result.inverse();
    713            }
    714            return result;
    715
    716        }
    717
    718        /**
    719         * Creates a transformation between two geocentric coordinate systems. This method is automatically invoked by
    720         * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust for
    721         * axis order and orientation, adjust for prime meridian, performs units conversion and apply Bursa Wolf
    722         * transformation if needed.
    723         *
    724         * @param sourceCRS
    725         *            Input coordinate system.
    726         * @param targetCRS
    727         *            Output coordinate system.
    728         * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
    729         * @throws TransformationException
    730         *             if no transformation path has been found.
    731         */
    732        private CRSTransformation createTransformation( final GeocentricCRS sourceCRS, final GeocentricCRS targetCRS )
    733                                throws TransformationException {
    734            LOG.logDebug( "Creating geocentric->geocetric transformation: from (source): " + sourceCRS.getIdentifier()
    735                          + " to(target): " + targetCRS.getIdentifier() );
    736            final GeodeticDatum sourceDatum = sourceCRS.getGeodeticDatum();
    737            final GeodeticDatum targetDatum = targetCRS.getGeodeticDatum();
    738            /*
    739             * if ( !PrimeMeridian.GREENWICH.equals( sourceDatum.getPrimeMeridian() ) || !PrimeMeridian.GREENWICH.equals(
    740             * targetDatum.getPrimeMeridian() ) ) { throw new TransformationException( "Rotation of prime meridian not yet
    741             * implemented" ); }
    742             */
    743            CRSTransformation result = null;
    744            if ( !sourceDatum.equals( targetDatum ) ) {
    745                final Ellipsoid sourceEllipsoid = sourceDatum.getEllipsoid();
    746                final Ellipsoid targetEllipsoid = targetDatum.getEllipsoid();
    747                /*
    748                 * If the two coordinate systems use different ellipsoid, convert from the source to target ellipsoid
    749                 * through the geocentric coordinate system.
    750                 */
    751                if ( sourceEllipsoid != null
    752                     && !sourceEllipsoid.equals( targetEllipsoid )
    753                     && ( ( sourceDatum.getWGS84Conversion() != null && sourceDatum.getWGS84Conversion().hasValues() ) || ( targetDatum.getWGS84Conversion() != null && targetDatum.getWGS84Conversion().hasValues() ) ) ) {
    754                    LOG.logDebug( "Creating helmert transformation: source(" + sourceCRS.getIdentifier() + ")->target("
    755                                  + targetCRS.getIdentifier() + ")." );
    756
    757                    // Transform between different ellipsoids using Bursa Wolf parameters.
    758                    Matrix tmp = swapAxis( sourceCRS, GeocentricCRS.WGS84 );
    759                    Matrix4d forwardAxisAlign = null;
    760                    if ( tmp != null ) {
    761                        forwardAxisAlign = new Matrix4d();
    762                        tmp.get( forwardAxisAlign );
    763                    }
    764                    final Matrix4d forwardToWGS = getWGS84Parameters( sourceDatum );
    765                    final Matrix4d inverseToWGS = getWGS84Parameters( targetDatum );
    766                    tmp = swapAxis( GeocentricCRS.WGS84, targetCRS );
    767                    Matrix4d resultMatrix = null;
    768                    if ( tmp != null ) {
    769                        resultMatrix = new Matrix4d();
    770                        tmp.get( resultMatrix );
    771                    }
    772                    if ( forwardAxisAlign == null && forwardToWGS == null && inverseToWGS == null && resultMatrix == null ) {
    773                        LOG.logDebug( "The given geocentric crs's do not need a helmert transformation (but they are not equal), returning identity" );
    774                        resultMatrix = new Matrix4d();
    775                        resultMatrix.setIdentity();
    776                    } else {
    777                        LOG.logDebug( "step1 matrix: \n " + forwardAxisAlign );
    778                        LOG.logDebug( "step2 matrix: \n " + forwardToWGS );
    779                        LOG.logDebug( "step3 matrix: \n " + inverseToWGS );
    780                        LOG.logDebug( "step4 matrix: \n " + resultMatrix );
    781                        if ( inverseToWGS != null ) {
    782                            inverseToWGS.invert(); // Invert in place.
    783                            LOG.logDebug( "inverseToWGS inverted matrix: \n " + inverseToWGS );
    784                        }
    785                        if ( resultMatrix != null ) {
    786                            if ( inverseToWGS != null ) {
    787                                resultMatrix.mul( inverseToWGS ); // step4 = step4*step3
    788                                LOG.logDebug( "resultMatrix (after mul with inverseToWGS): \n " + resultMatrix );
    789                            }
    790                            if ( forwardToWGS != null ) {
    791                                resultMatrix.mul( forwardToWGS ); // step4 = step4*step3*step2
    792                                LOG.logDebug( "resultMatrix (after mul with forwardToWGS2): \n " + resultMatrix );
    793                            }
    794                            if ( forwardAxisAlign != null ) {
    795                                resultMatrix.mul( forwardAxisAlign ); // step4 = step4*step3*step2*step1
    796                            }
    797                        } else if ( inverseToWGS != null ) {
    798                            resultMatrix = inverseToWGS;
    799                            if ( forwardToWGS != null ) {
    800                                resultMatrix.mul( forwardToWGS ); // step4 = step3*step2*step1
    801                                LOG.logDebug( "resultMatrix (after mul with forwardToWGS2): \n " + resultMatrix );
    802                            }
    803                            if ( forwardAxisAlign != null ) {
    804                                resultMatrix.mul( forwardAxisAlign ); // step4 = step3*step2*step1
    805                            }
    806                        } else if ( forwardToWGS != null ) {
    807                            resultMatrix = forwardToWGS;
    808                            if ( forwardAxisAlign != null ) {
    809                                resultMatrix.mul( forwardAxisAlign ); // step4 = step2*step1
    810                            }
    811                        } else {
    812                            resultMatrix = forwardAxisAlign;
    813                        }
    814                    }
    815
    816                    LOG.logDebug( "The resulting helmert transformation matrix: from( " + sourceCRS.getIdentifier()
    817                                  + ") to(" + targetCRS.getIdentifier() + ")\n " + resultMatrix );
    818                    result = new MatrixTransform( sourceCRS, targetCRS, resultMatrix, "Helmert-Transformation" );
    819
    820                }
    821            }
    822            if ( result == null ) {
    823                /*
    824                 * Swap axis order, and rotate the longitude coordinate if prime meridians are different.
    825                 */
    826                final Matrix matrix = swapAxis( sourceCRS, targetCRS );
    827                result = createMatrixTransform( sourceCRS, targetCRS, matrix );
    828                if ( LOG.isDebug() ) {
    829                    StringBuilder sb = new StringBuilder(
    830                                                          "Resulting axis alignment between source geocentric and target geocentric is:" );
    831                    if ( result == null ) {
    832                        sb.append( " not necessary" );
    833                    } else {
    834                        sb.append( "\n" ).append( ( (MatrixTransform) result ).getMatrix() );
    835                    }
    836                    LOG.logDebug( sb.toString() );
    837                }
    838            }
    839            return result;
    840        }
    841
    842        /**
    843         * Concatenates two existing transforms.
    844         *
    845         * @param first
    846         *            The first transform to apply to points.
    847         * @param second
    848         *            The second transform to apply to points.
    849         * @return The concatenated transform.
    850         *
    851         */
    852        private Transformation concatenateTransformations( Transformation first, Transformation second ) {
    853            if ( first == null ) {
    854                LOG.logDebug( "Concatenate: the first transform  is null." );
    855                return second;
    856            }
    857            if ( second == null ) {
    858                LOG.logDebug( "Concatenate: the second transform  is null." );
    859                return first;
    860            }
    861            // if one of the two is an identity transformation, just return the other.
    862            if ( first.isIdentity() ) {
    863                LOG.logDebug( "Concatenate:  the first transform  is the identity." );
    864                return second;
    865            }
    866            if ( second.isIdentity() ) {
    867                LOG.logDebug( "Concatenate:  the second transform is the identity." );
    868                return first;
    869            }
    870
    871            /*
    872             * If one transform is the inverse of the other, just return an identitiy transform.
    873             */
    874            if ( first.areInverse( second ) ) {
    875                LOG.logDebug( "Transformation1 and Transformation2 are inverse operations, returning null" );
    876                return null;
    877            }
    878
    879            /*
    880             * If both transforms use matrix, then we can create a single transform using the concatened matrix.
    881             */
    882            if ( first instanceof MatrixTransform && second instanceof MatrixTransform ) {
    883                GMatrix m1 = ( (MatrixTransform) first ).getMatrix();
    884                GMatrix m2 = ( (MatrixTransform) second ).getMatrix();
    885                if ( m1 == null ) {
    886                    if ( m2 == null ) {
    887                        // both matrices are null, just return the identiy matrix.
    888                        return new MatrixTransform( first.getSourceCRS(), first.getTargetCRS(),
    889                                                    new GMatrix( second.getTargetDimension() + 1,
    890                                                                 first.getSourceDimension() + 1 ) );
    891                    }
    892                    return second;
    893                } else if ( m2 == null ) {
    894                    return first;
    895                }
    896
    897                m2.mul( m1 );
    898                // m1.mul( m2 );
    899                LOG.logDebug( "Concatenate: both transforms are matrices, resulting multiply:\n" + m2 );
    900                MatrixTransform result = new MatrixTransform( first.getSourceCRS(), second.getTargetCRS(), m2 );
    901                // if ( result.isIdentity() ) {
    902                // return null;
    903                // }
    904                return result;
    905            }
    906
    907            /*
    908             * If one or both math transform are instance of {@link ConcatenatedTransform}, then concatenate
    909             * <code>tr1</code> or <code>tr2</code> with one of step transforms.
    910             */
    911            if ( first instanceof ConcatenatedTransform ) {
    912                final ConcatenatedTransform ctr = (ConcatenatedTransform) first;
    913                first = ctr.getFirstTransform();
    914                second = concatenateTransformations( ctr.getSecondTransform(), second );
    915            } else if ( second instanceof ConcatenatedTransform ) {
    916                final ConcatenatedTransform ctr = (ConcatenatedTransform) second;
    917                first = concatenateTransformations( first, ctr.getFirstTransform() );
    918                second = ctr.getSecondTransform();
    919            }
    920            // because of the concatenation one of the transformations may be null.
    921            if ( first == null ) {
    922                return second;
    923            }
    924            if ( second == null ) {
    925                return first;
    926            }
    927
    928            return new ConcatenatedTransform( first, second );
    929        }
    930
    931        /**
    932         * Creates an affine transform from a matrix.
    933         *
    934         * @param matrix
    935         *            The matrix used to define the affine transform.
    936         * @return The affine transform.
    937         * @throws TransformationException
    938         *             if the matrix is not affine.
    939         *
    940         */
    941        private MatrixTransform createMatrixTransform( CoordinateSystem sourceCRS, CoordinateSystem targetCRS,
    942                                                       final Matrix matrix )
    943                                throws TransformationException {
    944            if ( matrix == null ) {
    945                return null;
    946            }
    947            if ( matrix.isAffine() ) {// Affine transform are square.
    948                if ( matrix.getNumRow() == 3 && matrix.getNumCol() == 3 && !matrix.isIdentity() ) {
    949                    Matrix3d tmp = matrix.toAffineTransform();
    950                    return new MatrixTransform( sourceCRS, targetCRS, tmp );
    951                }
    952                return new MatrixTransform( sourceCRS, targetCRS, matrix );
    953            }
    954            throw new TransformationException( "Given matrix is not affine, cannot continue" );
    955        }
    956
    957        /**
    958         * @return the WGS84 parameters as an affine transform, or <code>null</code> if not available.
    959         */
    960        private Matrix4d getWGS84Parameters( final GeodeticDatum datum ) {
    961            final Helmert info = datum.getWGS84Conversion();
    962            if ( info != null && info.hasValues() ) {
    963                return info.getAsAffineTransform();
    964            }
    965            return null;
    966        }
    967
    968        /**
    969         * Concatenate two transformation steps.
    970         *
    971         * @param step1
    972         *            The first step, or <code>null</code> for the identity transform.
    973         * @param step2
    974         *            The second step, or <code>null</code> for the identity transform.
    975         * @return A concatenated transform, or <code>null</code> if all arguments was nul.
    976         */
    977        private Transformation concatenate( final Transformation step1, final Transformation step2 ) {
    978            if ( step1 == null )
    979                return step2;
    980            if ( step2 == null )
    981                return step1;
    982            return concatenateTransformations( step1, step2 );
    983        }
    984
    985        /**
    986         * Concatenate three transformation steps.
    987         *
    988         * @param step1
    989         *            The first step, or <code>null</code> for the identity transform.
    990         * @param step2
    991         *            The second step, or <code>null</code> for the identity transform.
    992         * @param step3
    993         *            The third step, or <code>null</code> for the identity transform.
    994         * @return A concatenated transform, or <code>null</code> if all arguments were <code>null</code>.
    995         */
    996        private Transformation concatenate( final Transformation step1, final Transformation step2,
    997                                            final Transformation step3 ) {
    998            if ( step1 == null ) {
    999                return concatenate( step2, step3 );
    1000            }
    1001            if ( step2 == null ) {
    1002                return concatenate( step1, step3 );
    1003            }
    1004            if ( step3 == null ) {
    1005                return concatenate( step1, step2 );
    1006            }
    1007            return concatenateTransformations( step1, concatenateTransformations( step2, step3 ) );
    1008        }
    1009
    1010        /**
    1011         * @return an affine transform between two coordinate systems. Only units and axis order (e.g. transforming from
    1012         *         (NORTH,WEST) to (EAST,NORTH)) are taken in account. Other attributes (especially the datum) must be
    1013         *         checked before invoking this method.
    1014         *
    1015         * @param sourceCRS
    1016         *            The source coordinate system.
    1017         * @param targetCRS
    1018         *            The target coordinate system.
    1019         */
    1020        private Matrix swapAxis( final CoordinateSystem sourceCRS, final CoordinateSystem targetCRS )
    1021                                throws TransformationException {
    1022            if ( LOG.isDebug() ) {
    1023                LOG.logDebug( "Creating swap matrix from: " + sourceCRS.getIdentifier() + " to: "
    1024                              + targetCRS.getIdentifier() );
    1025                LOG.logDebug( "Source Axis:\n" + Arrays.toString( sourceCRS.getAxis() ) );
    1026                LOG.logDebug( "Target Axis:\n" + Arrays.toString( targetCRS.getAxis() ) );
    1027            }
    1028            final Matrix matrix;
    1029            try {
    1030                matrix = new Matrix( sourceCRS.getAxis(), targetCRS.getAxis() );
    1031            } catch ( RuntimeException e ) {
    1032                throw new TransformationException( sourceCRS, targetCRS, e.getMessage() );
    1033            }
    1034            return matrix.isIdentity() ? null : matrix;
    1035        }
    1036
    1037        /**
    1038         * @return an affine transform between two geographic coordinate systems. Only units, axis order (e.g. transforming
    1039         *         from (NORTH,WEST) to (EAST,NORTH)) and prime meridian are taken in account. Other attributes (especially
    1040         *         the datum) must be checked before invoking this method.
    1041         *
    1042         * @param sourceCRS
    1043         *            The source coordinate system.
    1044         * @param targetCRS
    1045         *            The target coordinate system.
    1046         */
    1047        private Matrix swapAndRotateGeoAxis( final GeographicCRS sourceCRS, final GeographicCRS targetCRS )
    1048                                throws TransformationException {
    1049            if ( LOG.isDebug() ) {
    1050                LOG.logDebug( "Creating geo swap/rotate matrix from: " + sourceCRS.getIdentifier() + " to: "
    1051                              + targetCRS.getIdentifier() );
    1052            }
    1053            Matrix matrix = swapAxis( sourceCRS, targetCRS );
    1054            if ( !sourceCRS.getGeodeticDatum().getPrimeMeridian().equals( targetCRS.getGeodeticDatum().getPrimeMeridian() ) ) {
    1055                if ( matrix == null ) {
    1056                    matrix = new Matrix( sourceCRS.getDimension() + 1 );
    1057                }
    1058                Axis[] targetAxis = targetCRS.getAxis();
    1059                final int lastMatrixColumn = matrix.getNumCol() - 1;
    1060                for ( int i = 0; i < targetAxis.length; ++i ) {
    1061                    // Find longitude, and apply a translation if prime meridians are different.
    1062                    final int orientation = targetAxis[i].getOrientation();
    1063                    if ( Axis.AO_WEST == Math.abs( orientation ) ) {
    1064                        LOG.logDebug( "Adding prime-meridian translation to axis:" + targetAxis[i] );
    1065                        final double sourceLongitude = sourceCRS.getGeodeticDatum().getPrimeMeridian().getLongitudeAsRadian();
    1066                        final double targetLongitude = targetCRS.getGeodeticDatum().getPrimeMeridian().getLongitudeAsRadian();
    1067                        if ( Math.abs( sourceLongitude - targetLongitude ) > EPS11 ) {
    1068                            double translation = targetLongitude - sourceLongitude;
    1069                            if ( Axis.AO_WEST == orientation ) {
    1070                                translation = -translation;
    1071                            }
    1072                            // add the translation to the matrix translate element of this axis
    1073                            matrix.setElement( i, lastMatrixColumn, matrix.getElement( i, lastMatrixColumn ) - translation );
    1074
    1075                        }
    1076                    }
    1077                }
    1078            } else {
    1079                LOG.logDebug( "The primemeridians of the geographic crs's are equal, so no translation needed." );
    1080            }
    1081            return matrix;
    1082        }
    1083    }