001 //$HeadURL: https://sushibar/svn/deegree/base/trunk/src/org/deegree/model/csct/ct/CoordinateTransformationFactory.java $
002 /*---------------- FILE HEADER ------------------------------------------
003 This file is part of deegree.
004 Copyright (C) 2001-2008 by:
005 Department of Geography, University of Bonn
006 http://www.giub.uni-bonn.de/deegree/
007 lat/lon GmbH
008 http://www.lat-lon.de
009
010 This library is free software; you can redistribute it and/or
011 modify it under the terms of the GNU Lesser General Public
012 License as published by the Free Software Foundation; either
013 version 2.1 of the License, or (at your option) any later version.
014 This library is distributed in the hope that it will be useful,
015 but WITHOUT ANY WARRANTY; without even the implied warranty of
016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
017 Lesser General Public License for more details.
018 You should have received a copy of the GNU Lesser General Public
019 License along with this library; if not, write to the Free Software
020 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
021 Contact:
022
023 Andreas Poth
024 lat/lon GmbH
025 Aennchenstr. 19
026 53177 Bonn
027 Germany
028 E-Mail: poth@lat-lon.de
029
030 Prof. Dr. Klaus Greve
031 Department of Geography
032 University of Bonn
033 Meckenheimer Allee 166
034 53115 Bonn
035 Germany
036 E-Mail: greve@giub.uni-bonn.de
037 ---------------------------------------------------------------------------*/
038 package org.deegree.crs.transformations;
039
040 // OpenGIS dependencies
041 import javax.vecmath.GMatrix;
042 import javax.vecmath.Matrix3d;
043 import javax.vecmath.Matrix4d;
044
045 import org.deegree.crs.components.Axis;
046 import org.deegree.crs.components.Ellipsoid;
047 import org.deegree.crs.components.GeodeticDatum;
048 import org.deegree.crs.components.PrimeMeridian;
049 import org.deegree.crs.components.Unit;
050 import org.deegree.crs.coordinatesystems.CoordinateSystem;
051 import org.deegree.crs.coordinatesystems.GeocentricCRS;
052 import org.deegree.crs.coordinatesystems.GeographicCRS;
053 import org.deegree.crs.coordinatesystems.ProjectedCRS;
054 import org.deegree.crs.exceptions.TransformationException;
055 import org.deegree.crs.projections.ProjectionUtils;
056 import org.deegree.crs.utilities.Matrix;
057 import org.deegree.framework.log.ILogger;
058 import org.deegree.framework.log.LoggerFactory;
059
060 /**
061 * The <code>GeocentricTransform</code> class is the central access point for all transformations between different
062 * crs's.
063 * <p>
064 * It creates a transformation chain for two given CoordinateSystems by considering their type. For example the
065 * Transformation chain from EPSG:31466 ( a projected crs with underlying geographic crs epsg:4314 using the DHDN datum
066 * and the TransverseMercator Projection) to EPSG:28992 (another projected crs with underlying geographic crs epsg:4289
067 * using the 'new Amersfoort Datum' and the StereographicAzimuthal Projection) would result in following Transformation
068 * Chain:
069 * <ol>
070 * <li>Inverse projection - thus getting the coordinates in lat/lon for geographic crs epsg:4314</li>
071 * <li>Geodetic transformation - thus getting x-y-z coordinates for geographic crs epsg:4314</li>
072 * <li>WGS84 transformation -thus getting the x-y-z coordinates for the WGS84 datum </li>
073 * <li>Inverse WGS84 transformation -thus getting the x-y-z coordinates for the geodetic from epsg:4289</li>
074 * <li>Inverse geodetic - thus getting the lat/lon for epsg:4289</li>
075 * <li>projection - getting the coordinates (in meters) for epsg:28992</li>
076 * </ol>
077 *
078 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
079 *
080 * @author last edited by: $Author:$
081 *
082 * @version $Revision:$, $Date:$
083 *
084 */
085 public class TransformationFactory {
086 private static ILogger LOG = LoggerFactory.getLogger( TransformationFactory.class );
087
088 /**
089 * The default coordinate transformation factory. Will be constructed only when first needed.
090 */
091 private static TransformationFactory DEFAULT_INSTANCE = null;
092
093 private TransformationFactory() {
094 // nottin
095 }
096
097 /**
098 * @return the default coordinate transformation factory.
099 */
100 public static synchronized TransformationFactory getInstance() {
101 if ( DEFAULT_INSTANCE == null ) {
102 // DEFAULT_INSTANCE = new CoordinateTransformationFactory( MathTransformFactory.getDefault() );
103 DEFAULT_INSTANCE = new TransformationFactory();
104 }
105 return DEFAULT_INSTANCE;
106 }
107
108 /**
109 * Creates a transformation between two coordinate systems. This method will examine the coordinate systems in order
110 * to construct a transformation between them. This method may fail if no path between the coordinate systems is
111 * found.
112 *
113 * @param sourceCRS
114 * Input coordinate system.
115 * @param targetCRS
116 * Output coordinate system.
117 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
118 * @throws TransformationException
119 * @throws TransformationException
120 * if no transformation path has been found.
121 *
122 */
123 public CRSTransformation createFromCoordinateSystems( final CoordinateSystem sourceCRS,
124 final CoordinateSystem targetCRS )
125 throws TransformationException {
126 if ( sourceCRS == null ) {
127 throw new IllegalArgumentException( "The source CRS may not be null" );
128 }
129 if ( targetCRS == null ) {
130 throw new IllegalArgumentException( "The target CRS may not be null" );
131 }
132 if ( sourceCRS.equals( targetCRS ) ) {
133 LOG.logDebug( "Source crs and target crs are equal, no transformation needed (returning identity matrix)." );
134 final Matrix matrix = new Matrix( sourceCRS.getDimension() + 1 );
135 matrix.setIdentity();
136 return createAffineTransform( sourceCRS, targetCRS, matrix );
137 }
138 if ( ( sourceCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS && sourceCRS.getType() != CoordinateSystem.PROJECTED_CRS && sourceCRS.getType() != CoordinateSystem.GEOCENTRIC_CRS ) || ( targetCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS && targetCRS.getType() != CoordinateSystem.PROJECTED_CRS && targetCRS.getType() != CoordinateSystem.GEOCENTRIC_CRS ) ) {
139 throw new TransformationException( sourceCRS,
140 targetCRS,
141 "Either the target crs type or the source crs type was unknown" );
142 }
143 CRSTransformation result = null;
144 if ( sourceCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
145 /**
146 * Geographic --> Geographic, Projected or Geocentric
147 */
148 final GeographicCRS source = (GeographicCRS) sourceCRS;
149 if ( targetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
150 result = createTransformationStep( source, (GeographicCRS) targetCRS );
151 }
152 if ( targetCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
153 result = createTransformationStep( source, (ProjectedCRS) targetCRS );
154 }
155 if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
156 result = createTransformationStep( source, (GeocentricCRS) targetCRS );
157 }
158 } else if ( sourceCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
159 /**
160 * Projected --> Projected, Geographic or Geocentric
161 */
162 final ProjectedCRS source = (ProjectedCRS) sourceCRS;
163 if ( targetCRS.getType() == CoordinateSystem.PROJECTED_CRS ) {
164 result = createTransformationStep( source, (ProjectedCRS) targetCRS );
165 }
166 if ( targetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS ) {
167 result = createTransformationStep( source, (GeographicCRS) targetCRS );
168 }
169 if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
170 result = createTransformationStep( source, (GeocentricCRS) targetCRS );
171 }
172 } else if ( sourceCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
173 /**
174 * Geocentric --> Geocentric, Horizontal or Compound
175 */
176 final GeocentricCRS source = (GeocentricCRS) sourceCRS;
177 if ( targetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS ) {
178 result = createTransformationStep( source, (GeocentricCRS) targetCRS );
179 }
180 }
181 if ( result == null ) {
182 LOG.logDebug( "The resulting transformation was null, returning an identity matrix." );
183 final Matrix matrix = new Matrix( sourceCRS.getDimension() + 1 );
184 matrix.setIdentity();
185 result = createAffineTransform( sourceCRS, targetCRS, matrix );
186 }
187 if ( LOG.getLevel() == ILogger.LOG_DEBUG ) {
188 StringBuilder output = new StringBuilder( "The resulting transformation chain: \n" );
189 output = result.getTransformationPath( output );
190 LOG.logDebug( output.toString() );
191 }
192 return result;
193
194 }
195
196 /**
197 * Creates a transformation between two geographic coordinate systems. This method is automatically invoked by
198 * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust
199 * axis order and orientation (e.g. transforming from <code>(NORTH,WEST)</code> to <code>(EAST,NORTH)</code>),
200 * performs units conversion and apply Bursa Wolf transformation if needed.
201 *
202 * @param sourceCS
203 * Input coordinate system.
204 * @param targetCS
205 * Output coordinate system.
206 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
207 * @throws TransformationException
208 * if no transformation path has been found.
209 */
210 private CRSTransformation createTransformationStep( final GeographicCRS sourceCS, final GeographicCRS targetCS )
211 throws TransformationException {
212 if ( sourceCS.getGeodeticDatum().equals( targetCS.getGeodeticDatum() ) ) {
213 LOG.logDebug( "The datums of geographic (source): " + sourceCS.getIdAndName()
214 + " equals geographic (target): "
215 + targetCS.getIdAndName()
216 + " returning null" );
217 return null;
218 }
219 LOG.logDebug( "Creating geographic ->geographic transformation: from (source): " + sourceCS.getIdAndName()
220 + " to(target): "
221 + targetCS.getIdAndName() );
222 final GeodeticDatum sourceDatum = sourceCS.getGeodeticDatum();
223 final GeodeticDatum targetDatum = targetCS.getGeodeticDatum();
224 final Ellipsoid sourceEllipsoid = sourceDatum.getEllipsoid();
225 final Ellipsoid targetEllipsoid = targetDatum.getEllipsoid();
226 if ( sourceEllipsoid != null && !sourceEllipsoid.equals( targetEllipsoid ) ) {
227 /*
228 * If the two geographic coordinate systems use differents ellipsoid, convert from the source to target
229 * ellipsoid through the geocentric coordinate system.
230 */
231 final GeocentricCRS sourceGCS = new GeocentricCRS( sourceDatum,
232 sourceCS.getIdentifier(),
233 sourceCS.getName() + "Geocentric" );
234 final GeocentricCRS targetGCS = new GeocentricCRS( targetDatum,
235 targetCS.getIdentifier(),
236 targetCS.getName() + "Geocentric" );
237 CRSTransformation step1 = createTransformationStep( sourceCS, sourceGCS );
238 LOG.logDebug( "Step 1 from geographic to geographic:\n" + step1 );
239 CRSTransformation step2 = createTransformationStep( sourceGCS, targetGCS );
240 LOG.logDebug( "Step 2 from geographic to geographic:\n" + step2 );
241 CRSTransformation step3 = createTransformationStep( targetCS, targetGCS );
242 LOG.logDebug( "Step 3 from geographic to geographic:\n" + step3 );
243 if ( step3 != null ) {
244 step3.inverse();// call inversetransform from step 3.
245 }
246 return concatenate( step1, step2, step3 );
247 }
248
249 /*
250 * Swap axis order, and rotate the longitude coordinate if prime meridians are different.
251 */
252 final Matrix matrix = swapAndScaleGeoAxis( sourceCS, targetCS );
253 // TODO: We should ensure that longitude is in range [-180..+180°].
254 return createAffineTransform( sourceCS, targetCS, matrix );
255 // return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION, transform );
256 }
257
258 /**
259 * Creates a transformation between a geographic and a projected coordinate systems. This method is automatically
260 * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
261 *
262 * @param sourceCRS
263 * Input coordinate system.
264 * @param targetCRS
265 * Output coordinate system.
266 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
267 * @throws TransformationException
268 * if no transformation path has been found.
269 */
270 private CRSTransformation createTransformationStep( final GeographicCRS sourceCRS, final ProjectedCRS targetCRS )
271 throws TransformationException {
272
273 LOG.logDebug( "Creating geographic->projected transformation: from (source): " + sourceCRS.getIdAndName()
274 + " to(target): "
275 + targetCRS.getIdAndName() );
276 final GeographicCRS stepGeoCS = targetCRS.getGeographicCRS();
277
278 // should perform a datum shift
279 final CRSTransformation step1 = createTransformationStep( sourceCRS, stepGeoCS );
280 final CRSTransformation step2 = new ProjectionTransform( targetCRS );
281 return concatenate( step1, step2 );
282 }
283
284 /**
285 * Creates a transformation between a geographic and a geocentric coordinate systems. Since the source coordinate
286 * systems doesn't have a vertical axis, height above the ellipsoid is assumed equals to zero everywhere. This
287 * method is automatically invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
288 *
289 * @param sourceCS
290 * Input geographic coordinate system.
291 * @param targetCS
292 * Output coordinate system.
293 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
294 * @throws TransformationException
295 * if no transformation path has been found.
296 */
297 private CRSTransformation createTransformationStep( final GeographicCRS sourceCS, final GeocentricCRS targetCS )
298 throws TransformationException {
299 LOG.logDebug( "Creating geographic ->geocentric transformation: from (source): " + sourceCS.getIdAndName()
300 + " to(target): "
301 + targetCS.getIdAndName() );
302 if ( !PrimeMeridian.GREENWICH.equals( targetCS.getGeodeticDatum().getPrimeMeridian() ) ) {
303 throw new TransformationException( "Rotation of prime meridian not yet implemented" );
304 }
305 final CRSTransformation step1 = createAffineTransform( sourceCS,
306 targetCS,
307 swapAndScaleGeoAxis( sourceCS, GeographicCRS.WGS84 ) );
308
309 // TODO maybe something for a pool of transformations?
310 final CRSTransformation step2 = new GeocentricTransform( sourceCS, targetCS );
311 return createConcatenatedTransform( step1, step2 );
312 }
313
314 /**
315 * Creates a transformation between two projected coordinate systems. This method is automatically invoked by
316 * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust
317 * axis order and orientation. It also performs units conversion if it is the only extra change needed. Otherwise,
318 * it performs three steps:
319 *
320 * <ol>
321 * <li>Unproject <code>sourceCS</code>.</li>
322 * <li>Transform from <code>sourceCS.geographicCS</code> to <code>targetCS.geographicCS</code>.</li>
323 * <li>Project <code>targetCS</code>.</li>
324 * </ol>
325 *
326 * @param sourceCS
327 * Input coordinate system.
328 * @param targetCS
329 * Output coordinate system.
330 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
331 * @throws TransformationException
332 * if no transformation path has been found.
333 */
334 private CRSTransformation createTransformationStep( final ProjectedCRS sourceCS, final ProjectedCRS targetCS )
335 throws TransformationException {
336 LOG.logDebug( "Creating projected -> projected transformation: from (source): " + sourceCS.getIdAndName()
337 + " to(target): "
338 + targetCS.getIdAndName() );
339 if ( sourceCS.getProjection().equals( targetCS.getProjection() ) ) {
340 return null;
341 }
342 final GeographicCRS sourceGeo = sourceCS.getGeographicCRS();
343 final GeographicCRS targetGeo = targetCS.getGeographicCRS();
344 final CRSTransformation step1 = createTransformationStep( sourceCS, sourceGeo );
345 LOG.logDebug( "Step 1: " + step1 );
346 final CRSTransformation step2 = createTransformationStep( sourceGeo, targetGeo );
347 LOG.logDebug( "Step 2: " + step2 );
348 final CRSTransformation step3 = createTransformationStep( targetGeo, targetCS );
349 LOG.logDebug( "Step 3: " + step3 );
350 return concatenate( step1, step2, step3 );
351 }
352
353 /**
354 * Creates a transformation between a projected and a geocentric coordinate systems. This method is automatically
355 * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. This method doesn't need to be
356 * public since its decomposition in two step should be general enough.
357 *
358 * @param sourceCS
359 * Input projected coordinate system.
360 * @param targetCS
361 * Output coordinate system.
362 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
363 * @throws TransformationException
364 * if no transformation path has been found.
365 */
366 private CRSTransformation createTransformationStep( final ProjectedCRS sourceCS, final GeocentricCRS targetCS )
367 throws TransformationException {
368 LOG.logDebug( "Creating projected -> geocentric transformation: from (source): " + sourceCS.getIdAndName()
369 + " to(target): "
370 + targetCS.getIdAndName() );
371 final GeographicCRS sourceGCS = sourceCS.getGeographicCRS();
372
373 final CRSTransformation step1 = createTransformationStep( sourceCS, sourceGCS );
374 final CRSTransformation step2 = createTransformationStep( sourceGCS, targetCS );
375 return concatenate( step1, step2 );
376 }
377
378 /**
379 * Creates a transformation between a projected and a geographic coordinate systems. This method is automatically
380 * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation
381 * returns <code>{@link #createTransformationStep(GeographicCRS, ProjectedCRS)
382 * createTransformationStep}(targetCS, sourceCS) inverse)</code>.
383 *
384 * @param sourceCRS
385 * Input coordinate system.
386 * @param targetCRS
387 * Output coordinate system.
388 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code> or <code>null</code>
389 * if {@link ProjectedCRS#getGeographicCRS()}.equals targetCRS.
390 * @throws TransformationException
391 * if no transformation path has been found.
392 */
393 private CRSTransformation createTransformationStep( final ProjectedCRS sourceCRS, final GeographicCRS targetCRS )
394 throws TransformationException {
395 LOG.logDebug( "Creating projected->geographic transformation: from (source): " + sourceCRS.getIdAndName()
396 + " to(target): "
397 + targetCRS.getIdAndName() );
398 CRSTransformation result = createTransformationStep( targetCRS, sourceCRS );
399 result.inverse();
400 return result;
401
402 }
403
404 /**
405 * Creates a transformation between two geocentric coordinate systems. This method is automatically invoked by
406 * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust for
407 * axis order and orientation, adjust for prime meridian, performs units conversion and apply Bursa Wolf
408 * transformation if needed.
409 *
410 * @param sourceCS
411 * Input coordinate system.
412 * @param targetCS
413 * Output coordinate system.
414 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
415 * @throws TransformationException
416 * if no transformation path has been found.
417 */
418 private CRSTransformation createTransformationStep( final GeocentricCRS sourceCS, final GeocentricCRS targetCS )
419 throws TransformationException {
420 LOG.logDebug( "Creating geocentric->geocetric transformation: from (source): " + sourceCS.getIdAndName()
421 + " to(target): "
422 + targetCS.getIdAndName() );
423 final GeodeticDatum sourceHD = sourceCS.getGeodeticDatum();
424 final GeodeticDatum targetHD = targetCS.getGeodeticDatum();
425 if ( sourceHD.equals( targetHD ) ) {
426 return null;
427 // final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
428 // return createAffineTransform( sourceCS, targetCS, matrix );
429 }
430 if ( !PrimeMeridian.GREENWICH.equals( sourceHD.getPrimeMeridian() ) || !PrimeMeridian.GREENWICH.equals( targetHD.getPrimeMeridian() ) ) {
431 throw new TransformationException( "Rotation of prime meridian not yet implemented" );
432 }
433 // Transform between differents ellipsoids
434 // using Bursa Wolf parameters.
435 Matrix tmp = swapAndScaleAxis( sourceCS, GeocentricCRS.WGS84 );
436 Matrix4d step1 = null;
437 if ( tmp != null ) {
438 step1 = new Matrix4d();
439 tmp.get( step1 );
440 }
441 final Matrix4d step2 = getWGS84Parameters( sourceHD );
442 final Matrix4d step3 = getWGS84Parameters( targetHD );
443 tmp = swapAndScaleAxis( GeocentricCRS.WGS84, targetCS );
444 Matrix4d step4 = null;
445 if ( tmp != null ) {
446 step4 = new Matrix4d();
447 tmp.get( step4 );
448 }
449 if ( step1 == null && step2 == null && step3 == null && step4 == null ) {
450 LOG.logDebug( "The given geocentric crs's do not need a helmert transformation (but they are not equal), returning identity" );
451 step4 = new Matrix4d();
452 step4.setIdentity();
453 } else {
454 LOG.logDebug( "step1 matrix: \n " + step1 );
455 LOG.logDebug( "step2 matrix: \n " + step2 );
456 LOG.logDebug( "step3 matrix: \n " + step3 );
457 LOG.logDebug( "step4 matrix: \n " + step4 );
458 if ( step3 != null ) {
459 step3.invert(); // Invert in place.
460 LOG.logDebug( "step3 inverted matrix: \n " + step3 );
461 }
462 if ( step4 != null ) {
463 if ( step3 != null ) {
464 step4.mul( step3 ); // step4 = step4*step3
465 LOG.logDebug( "step4 (after mul3): \n " + step4 );
466 }
467 if ( step2 != null ) {
468 step4.mul( step2 ); // step4 = step4*step3*step2
469 LOG.logDebug( "step4 (after mul2): \n " + step4 );
470 }
471 if ( step1 != null ) {
472 step4.mul( step1 ); // step4 = step4*step3*step2*step1
473 }
474 } else if ( step3 != null ) {
475 step4 = step3;
476 if ( step2 != null ) {
477 step4.mul( step2 ); // step4 = step3*step2*step1
478 LOG.logDebug( "step4 (after mul2): \n " + step4 );
479 }
480 if ( step1 != null ) {
481 step4.mul( step1 ); // step4 = step3*step2*step1
482 }
483 } else if ( step2 != null ) {
484 step4 = step2;
485 if ( step1 != null ) {
486 step4.mul( step1 ); // step4 = step2*step1
487 }
488 } else {
489 step4 = step1;
490 }
491 }
492
493 LOG.logDebug( "The resulting geo-geo matrix: from( " + sourceCS.getIdentifier()
494 + ") to("
495 + targetCS.getIdentifier()
496 + ")\n "
497 + step4 );
498 return new MatrixTransform( sourceCS, targetCS, step4, CRSTransformation.createFromTo( sourceCS.getIdentifier(), targetCS.getIdentifier() ), "Helmert-Transformation" );
499 }
500
501 /**
502 * Concatenates two existing transforms.
503 *
504 * @param tr1
505 * The first transform to apply to points.
506 * @param tr2
507 * The second transform to apply to points.
508 * @return The concatenated transform.
509 *
510 */
511 private CRSTransformation createConcatenatedTransform( CRSTransformation tr1, CRSTransformation tr2 ) {
512 if ( tr1 == null ) {
513 return tr2;
514 }
515 if ( tr2 == null ) {
516 return tr1;
517 }
518 // if one of the two is an identity transformation, just return the other.
519 if ( tr1.isIdentity() ) {
520 return tr2;
521 }
522 if ( tr2.isIdentity() ) {
523 return tr1;
524 }
525
526 /*
527 * If one transform is the inverse of the other, just return an identitiy transform.
528 */
529 if ( tr1.areInverse( tr2 ) ) {
530 LOG.logDebug( "Transformation1 and Transformation2 are inverse operations, returning null" );
531 return null;
532 }
533
534 /*
535 * If both transforms use matrix, then we can create a single transform using the concatened matrix.
536 */
537 if ( tr1 instanceof MatrixTransform && tr2 instanceof MatrixTransform ) {
538 GMatrix m1 = ( (MatrixTransform) tr1 ).getMatrix();
539 GMatrix m2 = ( (MatrixTransform) tr2 ).getMatrix();
540 if ( m1 == null ) {
541 if ( m2 == null ) {
542 // both matrices are null, just return the identiy matrix.
543 return new MatrixTransform( tr1.getSourceCRS(),
544 tr1.getTargetCRS(),
545 new GMatrix( tr2.getTargetDimension() + 1, tr1.getSourceDimension() + 1 ) );
546 }
547 return tr2;
548 } else if ( m2 == null ) {
549 return tr1;
550 }
551 m2.mul( m1 );
552 return new MatrixTransform( tr1.getSourceCRS(), tr2.getTargetCRS(), m2 );
553 }
554
555 /*
556 * If one or both math transform are instance of {@link ConcatenatedTransform}, then concatenate <code>tr1</code>
557 * or <code>tr2</code> with one of step transforms.
558 */
559 if ( tr1 instanceof ConcatenatedTransform ) {
560 final ConcatenatedTransform ctr = (ConcatenatedTransform) tr1;
561 tr1 = ctr.getFirstTransform();
562 tr2 = createConcatenatedTransform( ctr.getSecondTransform(), tr2 );
563 }
564 if ( tr2 instanceof ConcatenatedTransform ) {
565 final ConcatenatedTransform ctr = (ConcatenatedTransform) tr2;
566 tr1 = createConcatenatedTransform( tr1, ctr.getFirstTransform() );
567 tr2 = ctr.getSecondTransform();
568 }
569
570 return new ConcatenatedTransform( tr1, tr2 );
571 }
572
573 /**
574 * Creates an affine transform from a matrix.
575 *
576 * @param matrix
577 * The matrix used to define the affine transform.
578 * @return The affine transform.
579 * @throws TransformationException
580 * if the matrix is not affine.
581 *
582 */
583 private MatrixTransform createAffineTransform( CoordinateSystem source, CoordinateSystem target, final Matrix matrix )
584 throws TransformationException {
585 if ( matrix == null ) {
586 return null;
587 }
588 if ( matrix.isAffine() ) {// Affine transform are square.
589 if ( matrix.getNumRow() == 3 && matrix.getNumCol() == 3 && !matrix.isIdentity() ) {
590 LOG.logDebug( "Creating affineTransform of incoming matrix:\n" + matrix );
591 Matrix3d tmp = matrix.toAffineTransform();
592 LOG.logDebug( "Resulting AffineTransform is:\n" + tmp );
593 return new MatrixTransform( source, target, tmp );
594 }
595 return new MatrixTransform( source, target, matrix );
596 }
597 throw new TransformationException( "Given matrix is not affine, cannot continue" );
598 }
599
600 /**
601 * Returns the WGS84 parameters as an affine transform, or <code>null</code> if not available.
602 */
603 private Matrix4d getWGS84Parameters( final GeodeticDatum datum ) {
604 final WGS84ConversionInfo info = datum.getWGS84Conversion();
605 if ( info != null && !( new WGS84ConversionInfo( "test" ) ).equals( info ) ) {
606 return info.getAsAffineTransform();
607 }
608 return null;
609 }
610
611 /**
612 * Concatenate two transformation steps.
613 *
614 * @param step1
615 * The first step, or <code>null</code> for the identity transform.
616 * @param step2
617 * The second step, or <code>null</code> for the identity transform.
618 * @return A concatenated transform, or <code>null</code> if all arguments was nul.
619 */
620 private CRSTransformation concatenate( final CRSTransformation step1, final CRSTransformation step2 ) {
621 if ( step1 == null )
622 return step2;
623 if ( step2 == null )
624 return step1;
625 // if ( !step1.getTargetCRS().equals( step2.getSourceCRS() ) ) {
626 // throw new IllegalArgumentException( "The given crsTranformation can not be concatenated because the first tranformations targetCRS: " + step1.getTargetCRS()
627 // .getIdentifier()
628 // + " does not match the second transformations sourceCRS: "
629 // + step2.getSourceCRS().getIdentifier() );
630 // }
631 return createConcatenatedTransform( step1, step2 );
632 }
633
634 /**
635 * Concatenate three transformation steps.
636 *
637 * @param step1
638 * The first step, or <code>null</code> for the identity transform.
639 * @param step2
640 * The second step, or <code>null</code> for the identity transform.
641 * @param step3
642 * The third step, or <code>null</code> for the identity transform.
643 * @return A concatenated transform, or <code>null</code> if all arguments were <code>null</code>.
644 */
645 private CRSTransformation concatenate( final CRSTransformation step1, final CRSTransformation step2,
646 final CRSTransformation step3 ) {
647 if ( step1 == null ){
648 return concatenate( step2, step3 );
649 }
650 if ( step2 == null ){
651 return concatenate( step1, step3 );
652 }
653 if ( step3 == null ) {
654 return concatenate( step1, step2 );
655 }
656 return createConcatenatedTransform( step1, createConcatenatedTransform( step2, step3 ) );
657 }
658
659 /**
660 * Returns an affine transform between two coordinate systems. Only units and axis order (e.g. transforming from
661 * (NORTH,WEST) to (EAST,NORTH)) are taken in account. Other attributes (especially the datum) must be checked
662 * before invoking this method.
663 *
664 * @param sourceCS
665 * The source coordinate system. If <code>null</code>, then (x,y,z,t) axis order is assumed.
666 * @param targetCS
667 * The target coordinate system. If <code>null</code>, then (x,y,z,t) axis order is assumed.
668 */
669 private Matrix swapAndScaleAxis( final CoordinateSystem sourceCS, final CoordinateSystem targetCS )
670 throws TransformationException {
671 LOG.logDebug( "Creating swap matrix from: " + sourceCS.getIdAndName() + " to: " + targetCS.getIdAndName() );
672 final Axis[] sourceAxis = sourceCS.getAxis();
673 final Axis[] targetAxis = targetCS.getAxis();
674 final Matrix matrix;
675 try {
676 matrix = Matrix.createAffineTransform( sourceAxis, targetAxis );
677 } catch ( RuntimeException e ) {
678 throw new TransformationException( sourceCS, targetCS, e.getMessage() );
679 }
680 // Convert units if necessary
681 final int dimension = matrix.getNumRow() - 1;
682 for ( int i = 0; i < dimension; i++ ) {
683 final Unit sourceUnit = sourceAxis[i].getUnits();
684 final Unit targetUnit = targetAxis[i].getUnits();
685 double offset = 0;
686 double scale = 1;
687 if ( !sourceUnit.equals( targetUnit ) ) {
688 offset = targetUnit.convert( 0, sourceUnit );
689 scale = targetUnit.convert( 1, sourceUnit ) - offset;
690 }
691 matrix.setElement( i, i, scale * matrix.getElement( i, i ) );
692 matrix.setElement( i, dimension, scale * matrix.getElement( i, dimension ) + offset );
693 }
694 if ( matrix.isIdentity() ) {
695 return null;
696 }
697 return matrix;
698 }
699
700 /**
701 * Returns an affine transform between two geographic coordinate systems. Only units, axis order (e.g. transforming
702 * from (NORTH,WEST) to (EAST,NORTH)) and prime meridian are taken in account. Other attributes (especially the
703 * datum) must be checked before invoking this method.
704 *
705 * @param sourceCS
706 * The source coordinate system.
707 * @param targetCS
708 * The target coordinate system.
709 */
710 private Matrix swapAndScaleGeoAxis( final GeographicCRS sourceCS, final GeographicCRS targetCS )
711 throws TransformationException {
712 LOG.logDebug( "Creating geo swap matrix from: " + sourceCS.getIdAndName() + " to: " + targetCS.getIdAndName() );
713 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
714 if ( matrix != null ) {
715 Axis[] targetAxis = targetCS.getAxis();
716 final int lastMatrixColumn = matrix.getNumCol() - 1;
717 for ( int i = targetCS.getDimension(); --i >= 0; ) {
718 // Find longitude, and apply a translation if prime meridians are different.
719 final int orientation = targetAxis[i].getOrientation();
720 if ( Axis.AO_EAST == Math.abs( orientation ) ) {
721 final Unit unit = targetAxis[i].getUnits();
722 final double sourceLongitude = sourceCS.getGeodeticDatum().getPrimeMeridian().getLongitude( unit );
723 final double targetLongitude = targetCS.getGeodeticDatum().getPrimeMeridian().getLongitude( unit );
724 if ( Math.abs( sourceLongitude - targetLongitude ) > ProjectionUtils.EPS11 ) {
725 double translation = targetLongitude - sourceLongitude;
726 if ( Axis.AO_WEST == orientation ) {
727 translation = -translation;
728 }
729 // add the translation to the matrix translate element of this axis
730 matrix.setElement( i, lastMatrixColumn, matrix.getElement( i, lastMatrixColumn ) - translation );
731 }
732 }
733 }
734 }
735 return matrix;
736 }
737 }