001 //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.3_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 }