001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_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    }