001 //$HeadURL: svn+ssh://rbezema@svn.wald.intevation.org/deegree/base/tags/2.1/src/org/deegree/model/csct/ct/CoordinateTransformationFactory.java $
002 /*---------------- FILE HEADER ------------------------------------------
003
004 This file is part of deegree.
005 Copyright (C) 2001 by:
006 EXSE, Department of Geography, University of Bonn
007 http://www.giub.uni-bonn.de/exse/
008 lat/lon GmbH
009 http://www.lat-lon.de
010
011 It has been implemented within SEAGIS - An OpenSource implementation of OpenGIS specification
012 (C) 2001, Institut de Recherche pour le D�veloppement (http://sourceforge.net/projects/seagis/)
013 SEAGIS Contacts: Surveillance de l'Environnement Assist�e par Satellite
014 Institut de Recherche pour le D�veloppement / US-Espace
015 mailto:seasnet@teledetection.fr
016
017
018 This library is free software; you can redistribute it and/or
019 modify it under the terms of the GNU Lesser General Public
020 License as published by the Free Software Foundation; either
021 version 2.1 of the License, or (at your option) any later version.
022
023 This library is distributed in the hope that it will be useful,
024 but WITHOUT ANY WARRANTY; without even the implied warranty of
025 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
026 Lesser General Public License for more details.
027
028 You should have received a copy of the GNU Lesser General Public
029 License along with this library; if not, write to the Free Software
030 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
031
032 Contact:
033
034 Andreas Poth
035 lat/lon GmbH
036 Aennchenstr. 19
037 53115 Bonn
038 Germany
039 E-Mail: poth@lat-lon.de
040
041 Klaus Greve
042 Department of Geography
043 University of Bonn
044 Meckenheimer Allee 166
045 53115 Bonn
046 Germany
047 E-Mail: klaus.greve@uni-bonn.de
048
049
050 ---------------------------------------------------------------------------*/
051 package org.deegree.model.csct.ct;
052
053 // OpenGIS dependencies
054 import javax.media.jai.ParameterList;
055 import javax.vecmath.SingularMatrixException;
056
057 import org.deegree.model.csct.cs.AxisOrientation;
058 import org.deegree.model.csct.cs.CompoundCoordinateSystem;
059 import org.deegree.model.csct.cs.CoordinateSystem;
060 import org.deegree.model.csct.cs.Ellipsoid;
061 import org.deegree.model.csct.cs.GeocentricCoordinateSystem;
062 import org.deegree.model.csct.cs.GeographicCoordinateSystem;
063 import org.deegree.model.csct.cs.HorizontalCoordinateSystem;
064 import org.deegree.model.csct.cs.HorizontalDatum;
065 import org.deegree.model.csct.cs.PrimeMeridian;
066 import org.deegree.model.csct.cs.ProjectedCoordinateSystem;
067 import org.deegree.model.csct.cs.Projection;
068 import org.deegree.model.csct.cs.TemporalCoordinateSystem;
069 import org.deegree.model.csct.cs.VerticalCoordinateSystem;
070 import org.deegree.model.csct.cs.VerticalDatum;
071 import org.deegree.model.csct.cs.WGS84ConversionInfo;
072 import org.deegree.model.csct.pt.Dimensioned;
073 import org.deegree.model.csct.pt.Matrix;
074 import org.deegree.model.csct.resources.OpenGIS;
075 import org.deegree.model.csct.resources.Utilities;
076 import org.deegree.model.csct.resources.css.ResourceKeys;
077 import org.deegree.model.csct.resources.css.Resources;
078 import org.deegree.model.csct.units.Unit;
079
080 /**
081 * Creates coordinate transformations.
082 *
083 * @version 1.0
084 * @author OpenGIS (www.opengis.org)
085 * @author Martin Desruisseaux
086 *
087 * @see org.opengis.ct.CT_CoordinateTransformationFactory
088 */
089 public class CoordinateTransformationFactory {
090 /**
091 * The default coordinate transformation factory.
092 * Will be constructed only when first needed.
093 */
094 private static CoordinateTransformationFactory DEFAULT;
095
096 /**
097 * Number for temporary created objects. This number is
098 * incremented each time {@link #getTemporaryName} is
099 * invoked.
100 */
101 private static volatile int temporaryID;
102
103 /**
104 * The underlying math transform factory.
105 */
106 private final MathTransformFactory factory;
107
108 /**
109 * Construct a coordinate transformation factory.
110 *
111 * @param factory The math transform factory to use.
112 */
113 public CoordinateTransformationFactory( final MathTransformFactory factory ) {
114 this.factory = factory;
115 }
116
117 /**
118 * Returns the default coordinate transformation factory.
119 */
120 public static synchronized CoordinateTransformationFactory getDefault() {
121 if ( DEFAULT == null ) {
122 DEFAULT = new CoordinateTransformationFactory( MathTransformFactory.getDefault() );
123 }
124 return DEFAULT;
125 }
126
127 /**
128 * Returns the underlying math transform factory. This factory
129 * is used for constructing {@link MathTransform} objects for
130 * all {@link CoordinateTransformation}.
131 */
132 public final MathTransformFactory getMathTransformFactory() {
133 return factory;
134 }
135
136 /**
137 * Creates a transformation between two coordinate systems.
138 * This method will examine the coordinate systems in order to construct a
139 * transformation between them. This method may fail if no path between the
140 * coordinate systems is found.
141 *
142 * @param sourceCS Input coordinate system.
143 * @param targetCS Output coordinate system.
144 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
145 * @throws CannotCreateTransformException if no transformation path has been found.
146 *
147 */
148 public CoordinateTransformation createFromCoordinateSystems( final CoordinateSystem sourceCS,
149 final CoordinateSystem targetCS )
150 throws CannotCreateTransformException {
151 /////////////////////////////////////////////////////////////////////
152 //// ////
153 //// Geographic --> Geographic, Projected or Geocentric ////
154 //// ////
155 /////////////////////////////////////////////////////////////////////
156 if ( sourceCS instanceof GeographicCoordinateSystem ) {
157 final GeographicCoordinateSystem source = (GeographicCoordinateSystem) sourceCS;
158 if ( targetCS instanceof GeographicCoordinateSystem ) {
159 return createTransformationStep( source, (GeographicCoordinateSystem) targetCS );
160 }
161 if ( targetCS instanceof ProjectedCoordinateSystem ) {
162 return createTransformationStep( source, (ProjectedCoordinateSystem) targetCS );
163 }
164 if ( targetCS instanceof GeocentricCoordinateSystem ) {
165 return createTransformationStep( source, (GeocentricCoordinateSystem) targetCS );
166 }
167 }
168 /////////////////////////////////////////////////////////////////////
169 //// ////
170 //// Projected --> Projected, Geographic or Geocentric ////
171 //// ////
172 /////////////////////////////////////////////////////////////////////
173 if ( sourceCS instanceof ProjectedCoordinateSystem ) {
174 final ProjectedCoordinateSystem source = (ProjectedCoordinateSystem) sourceCS;
175 if ( targetCS instanceof ProjectedCoordinateSystem ) {
176 return createTransformationStep( source, (ProjectedCoordinateSystem) targetCS );
177 }
178 if ( targetCS instanceof GeographicCoordinateSystem ) {
179 return createTransformationStep( source, (GeographicCoordinateSystem) targetCS );
180 }
181 if ( targetCS instanceof GeocentricCoordinateSystem ) {
182 return createTransformationStep( source, (GeocentricCoordinateSystem) targetCS );
183 }
184 }
185 /////////////////////////////////////////////////////////////////////
186 //// ////
187 //// Geocentric --> Geocentric, Horizontal or Compound ////
188 //// ////
189 /////////////////////////////////////////////////////////////////////
190 if ( sourceCS instanceof GeocentricCoordinateSystem ) {
191 final GeocentricCoordinateSystem source = (GeocentricCoordinateSystem) sourceCS;
192 if ( targetCS instanceof GeocentricCoordinateSystem ) {
193 return createTransformationStep( source, (GeocentricCoordinateSystem) targetCS );
194 }
195 try {
196 return createFromCoordinateSystems( targetCS, sourceCS ).inverse();
197 } catch ( TransformException exception ) {
198 final CannotCreateTransformException e = new CannotCreateTransformException(
199 sourceCS,
200 targetCS );
201 e.initCause( exception );
202 throw e;
203 }
204 }
205 /////////////////////////////////////////
206 //// ////
207 //// Vertical --> Vertical ////
208 //// ////
209 /////////////////////////////////////////
210 if ( sourceCS instanceof VerticalCoordinateSystem ) {
211 final VerticalCoordinateSystem source = (VerticalCoordinateSystem) sourceCS;
212 if ( targetCS instanceof VerticalCoordinateSystem ) {
213 return createTransformationStep( source, (VerticalCoordinateSystem) targetCS );
214 }
215 }
216 /////////////////////////////////////////
217 //// ////
218 //// Temporal --> Temporal ////
219 //// ////
220 /////////////////////////////////////////
221 if ( sourceCS instanceof TemporalCoordinateSystem ) {
222 final TemporalCoordinateSystem source = (TemporalCoordinateSystem) sourceCS;
223 if ( targetCS instanceof TemporalCoordinateSystem ) {
224 return createTransformationStep( source, (TemporalCoordinateSystem) targetCS );
225 }
226 }
227 ///////////////////////////////////////////
228 //// ////
229 //// Compound --> various CS ////
230 //// ////
231 ///////////////////////////////////////////
232 if ( sourceCS instanceof CompoundCoordinateSystem ) {
233 final CompoundCoordinateSystem source = (CompoundCoordinateSystem) sourceCS;
234 if ( targetCS instanceof CompoundCoordinateSystem ) {
235 return createTransformationStep( source, (CompoundCoordinateSystem) targetCS );
236 }
237 if ( targetCS instanceof GeocentricCoordinateSystem ) {
238 return createTransformationStep( source, (GeocentricCoordinateSystem) targetCS );
239 }
240 /*
241 * Try a loosely transformation. For example, the source CS may be
242 * a geographic + vertical coordinate systems, will the target CS
243 * may be only the geographic part. The code below will try to
244 * discart one or more dimension.
245 */
246 final CoordinateSystem headSourceCS = source.getHeadCS();
247 final CoordinateSystem tailSourceCS = source.getTailCS();
248 final int dimHeadCS = headSourceCS.getDimension();
249 final int dimSource = source.getDimension();
250 CoordinateTransformation step2;
251 int lower, upper;
252 try {
253 lower = 0;
254 upper = dimHeadCS;
255 step2 = createFromCoordinateSystems( headSourceCS, targetCS );
256 } catch ( CannotCreateTransformException exception ) {
257 /*
258 * If we can't construct a transformation from the head CS,
259 * then try a transformation from the tail CS. If this step
260 * fails also, then the head CS will be taken as the raison
261 * for the failure.
262 */
263 try {
264 lower = dimHeadCS;
265 upper = dimSource;
266 step2 = createFromCoordinateSystems( tailSourceCS, targetCS );
267 } catch ( CannotCreateTransformException ignore ) {
268 CannotCreateTransformException e = new CannotCreateTransformException(
269 sourceCS,
270 targetCS );
271 e.initCause( exception );
272 throw e;
273 }
274 }
275 /*
276 * A coordinate transformation from the head or tail part of 'sourceCS'
277 * has been succesfully contructed. Now, build a matrix transform that
278 * will select only the corresponding ordinates from input arrays, and
279 * pass them to the transform.
280 */
281 final MathTransform step1 = factory.createSubMathTransform(
282 lower,
283 upper,
284 factory.createIdentityTransform( dimSource ) );
285 final MathTransform transform = factory.createConcatenatedTransform(
286 step1,
287 step2.getMathTransform() );
288 return createFromMathTransform( sourceCS, targetCS, step2.getTransformType(), transform );
289 }
290 throw new CannotCreateTransformException( sourceCS, targetCS );
291 }
292
293 /**
294 * Creates a transformation between two temporal coordinate systems. This
295 * method is automatically invoked by {@link #createFromCoordinateSystems
296 * createFromCoordinateSystems(...)}. The default implementation checks if
297 * both coordinate systems use the same datum, and then adjusts for axis
298 * orientation, units and epoch.
299 *
300 * @param sourceCS Input coordinate system.
301 * @param targetCS Output coordinate system.
302 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
303 * @throws CannotCreateTransformException if no transformation path has been found.
304 */
305 protected CoordinateTransformation createTransformationStep(
306 final TemporalCoordinateSystem sourceCS,
307 final TemporalCoordinateSystem targetCS )
308 throws CannotCreateTransformException {
309 if ( Utilities.equals( sourceCS.getTemporalDatum(), targetCS.getTemporalDatum() ) ) {
310 // Compute the epoch shift
311 double epochShift = sourceCS.getEpoch().getTime() - targetCS.getEpoch().getTime();
312 epochShift = targetCS.getUnits( 0 ).convert( epochShift / ( 24 * 60 * 60 * 1000 ),
313 Unit.DAY );
314
315 // Get the affine transform, add the epoch
316 // shift and returns the resulting transform.
317 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
318 final int translationColumn = matrix.getNumCol() - 1;
319 if ( translationColumn >= 0 ) // Paranoiac check: should always be 1.
320 {
321 final double translation = matrix.getElement( 0, translationColumn );
322 matrix.setElement( 0, translationColumn, translation + epochShift );
323 }
324 final MathTransform transform = factory.createAffineTransform( matrix );
325 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION, transform );
326 }
327 throw new CannotCreateTransformException( sourceCS, targetCS );
328 }
329
330 /**
331 * Creates a transformation between two vertical coordinate systems. This
332 * method is automatically invoked by {@link #createFromCoordinateSystems
333 * createFromCoordinateSystems(...)}. The default implementation checks if
334 * both coordinate systems use the same datum, and then adjusts for axis
335 * orientation and units.
336 *
337 * @param sourceCS Input coordinate system.
338 * @param targetCS Output coordinate system.
339 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
340 * @throws CannotCreateTransformException if no transformation path has been found.
341 */
342 protected CoordinateTransformation createTransformationStep(
343 final VerticalCoordinateSystem sourceCS,
344 final VerticalCoordinateSystem targetCS )
345 throws CannotCreateTransformException {
346 if ( Utilities.equals( sourceCS.getVerticalDatum(), targetCS.getVerticalDatum() ) ) {
347 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
348 final MathTransform transform = factory.createAffineTransform( matrix );
349 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION, transform );
350 }
351 throw new CannotCreateTransformException( sourceCS, targetCS );
352 }
353
354 /**
355 * Creates a transformation between two geographic coordinate systems. This
356 * method is automatically invoked by {@link #createFromCoordinateSystems
357 * createFromCoordinateSystems(...)}. The default implementation can adjust
358 * axis order and orientation (e.g. transforming from <code>(NORTH,WEST)</code>
359 * to <code>(EAST,NORTH)</code>), performs units conversion and apply Bursa Wolf
360 * transformation if needed.
361 *
362 * @param sourceCS Input coordinate system.
363 * @param targetCS Output coordinate system.
364 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
365 * @throws CannotCreateTransformException if no transformation path has been found.
366 */
367 protected CoordinateTransformation createTransformationStep(
368 final GeographicCoordinateSystem sourceCS,
369 final GeographicCoordinateSystem targetCS )
370 throws CannotCreateTransformException {
371 final HorizontalDatum sourceDatum = sourceCS.getHorizontalDatum();
372 final HorizontalDatum targetDatum = targetCS.getHorizontalDatum();
373 final Ellipsoid sourceEllipsoid = sourceDatum.getEllipsoid();
374 final Ellipsoid targetEllipsoid = targetDatum.getEllipsoid();
375 if ( !Utilities.equals( sourceEllipsoid, targetEllipsoid ) )
376 try {
377 /*
378 * If the two geographic coordinate systems use differents ellipsoid,
379 * convert from the source to target ellipsoid through the geocentric
380 * coordinate system.
381 */
382 final String name = getTemporaryName();
383 final GeocentricCoordinateSystem gcs1 = new GeocentricCoordinateSystem( name,
384 sourceDatum );
385 final GeocentricCoordinateSystem gcs3 = new GeocentricCoordinateSystem( name,
386 targetDatum );
387 CoordinateTransformation step1 = createTransformationStep( sourceCS, gcs1 );
388 CoordinateTransformation step2 = createTransformationStep( gcs1, gcs3 );
389 CoordinateTransformation step3 = createTransformationStep( targetCS, gcs3 ).inverse();
390 return concatenate( step1, step2, step3 );
391
392 } catch ( TransformException exception ) {
393 CannotCreateTransformException e = new CannotCreateTransformException( sourceCS,
394 targetCS );
395 throw e;
396 }
397 /*
398 * Swap axis order, and rotate the longitude
399 * coordinate if prime meridians are different.
400 */
401 final Matrix matrix = swapAndScaleGeoAxis( sourceCS, targetCS );
402 // TODO: We should ensure that longitude is in range [-180..+180�].
403 final MathTransform transform = factory.createAffineTransform( matrix );
404 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION, transform );
405 }
406
407 /**
408 * Creates a transformation between two projected coordinate systems. This
409 * method is automatically invoked by {@link #createFromCoordinateSystems
410 * createFromCoordinateSystems(...)}. The default implementation can adjust
411 * axis order and orientation. It also performs units conversion if it
412 * is the only extra change needed. Otherwise, it performs three steps:
413 *
414 * <ol>
415 * <li>Unproject <code>sourceCS</code>.</li>
416 * <li>Transform from <code>sourceCS.geographicCS</code> to <code>targetCS.geographicCS</code>.</li>
417 * <li>Project <code>targetCS</code>.</li>
418 * </ol>
419 *
420 * @param sourceCS Input coordinate system.
421 * @param targetCS Output coordinate system.
422 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
423 * @throws CannotCreateTransformException if no transformation path has been found.
424 */
425 protected CoordinateTransformation createTransformationStep(
426 final ProjectedCoordinateSystem sourceCS,
427 final ProjectedCoordinateSystem targetCS )
428 throws CannotCreateTransformException {
429 if ( Utilities.equals( sourceCS.getProjection(), targetCS.getProjection() )
430 && Utilities.equals( sourceCS.getHorizontalDatum(), targetCS.getHorizontalDatum() ) ) {
431 // This special case is necessary for createTransformationStep(GeographicCS,ProjectedCS).
432 // 'swapAndScaleAxis(...) takes care of axis orientation and units. Datum and projection
433 // have just been checked above. Prime meridien is not checked (TODO: should we do???)
434 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
435 final MathTransform transform = factory.createAffineTransform( matrix );
436 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION, transform );
437 }
438 final GeographicCoordinateSystem sourceGeo = sourceCS.getGeographicCoordinateSystem();
439 final GeographicCoordinateSystem targetGeo = targetCS.getGeographicCoordinateSystem();
440 final CoordinateTransformation step1 = createTransformationStep( sourceCS, sourceGeo );
441 final CoordinateTransformation step2 = createTransformationStep( sourceGeo, targetGeo );
442 final CoordinateTransformation step3 = createTransformationStep( targetGeo, targetCS );
443 return concatenate( step1, step2, step3 );
444 }
445
446 /**
447 * Creates a transformation between a geographic and a projected coordinate systems.
448 * This method is automatically invoked by {@link #createFromCoordinateSystems
449 * createFromCoordinateSystems(...)}.
450 *
451 * @param sourceCS Input coordinate system.
452 * @param targetCS Output coordinate system.
453 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
454 * @throws CannotCreateTransformException if no transformation path has been found.
455 */
456 protected CoordinateTransformation createTransformationStep(
457 final GeographicCoordinateSystem sourceCS,
458 final ProjectedCoordinateSystem targetCS )
459 throws CannotCreateTransformException {
460
461 final ProjectedCoordinateSystem stepProjCS = normalize( targetCS );
462 final GeographicCoordinateSystem stepGeoCS = stepProjCS.getGeographicCoordinateSystem();
463
464 final Projection projection = stepProjCS.getProjection();
465 final MathTransform mapProjection = factory.createParameterizedTransform( projection );
466 // should perform a datum shift
467 final CoordinateTransformation step1 = createTransformationStep( sourceCS, stepGeoCS );
468 final CoordinateTransformation step2 = createFromMathTransform( stepGeoCS, stepProjCS,
469 TransformType.CONVERSION,
470 mapProjection );
471
472 // TODO
473 // this is possibly an error because stepProjCS and targetCS
474 // are identical
475 //final CoordinateTransformation step3 = createTransformationStep( stepProjCS, targetCS );
476 //return concatenate( step1, step2, step3 );
477 // it seems it must look like this
478 return concatenate( step1, step2, null );
479 }
480
481 /**
482 * Creates a transformation between a projected and a geographic coordinate systems.
483 * This method is automatically invoked by {@link #createFromCoordinateSystems
484 * createFromCoordinateSystems(...)}. The default implementation returns
485 * <code>{@link #createTransformationStep(GeographicCoordinateSystem, ProjectedCoordinateSystem)
486 * createTransformationStep}(targetCS, sourceCS).{@link MathTransform#inverse() inverse()}</code>.
487 *
488 * @param sourceCS Input coordinate system.
489 * @param targetCS Output coordinate system.
490 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
491 * @throws CannotCreateTransformException if no transformation path has been found.
492 */
493 protected CoordinateTransformation createTransformationStep(
494 final ProjectedCoordinateSystem sourceCS,
495 final GeographicCoordinateSystem targetCS )
496 throws CannotCreateTransformException {
497 try {
498 return createTransformationStep( targetCS, sourceCS ).inverse();
499 } catch ( NoninvertibleTransformException exception ) {
500 final CannotCreateTransformException e = new CannotCreateTransformException( sourceCS,
501 targetCS );
502 throw e;
503 }
504 }
505
506 /**
507 * Creates a transformation between two geocentric coordinate systems. This
508 * method is automatically invoked by {@link #createFromCoordinateSystems
509 * createFromCoordinateSystems(...)}. The default implementation can adjust
510 * for axis order and orientation, adjust for prime meridian, performs units
511 * conversion and apply Bursa Wolf transformation if needed.
512 *
513 * @param sourceCS Input coordinate system.
514 * @param targetCS Output coordinate system.
515 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
516 * @throws CannotCreateTransformException if no transformation path has been found.
517 */
518 protected CoordinateTransformation createTransformationStep(
519 final GeocentricCoordinateSystem sourceCS,
520 final GeocentricCoordinateSystem targetCS )
521 throws CannotCreateTransformException {
522 final HorizontalDatum sourceHD = sourceCS.getHorizontalDatum();
523 final HorizontalDatum targetHD = targetCS.getHorizontalDatum();
524 if ( Utilities.equals( sourceHD, targetHD ) ) {
525 if ( Utilities.equals( sourceCS.getPrimeMeridian(), targetCS.getPrimeMeridian() ) ) {
526 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
527 final MathTransform transform = factory.createAffineTransform( matrix );
528 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION,
529 transform );
530 }
531 // If prime meridians are not the same, performs the full transformation.
532 }
533 if ( !PrimeMeridian.GREENWICH.equals( sourceCS.getPrimeMeridian() )
534 || !PrimeMeridian.GREENWICH.equals( targetCS.getPrimeMeridian() ) ) {
535 throw new CannotCreateTransformException(
536 "Rotation of prime meridian not yet implemented" );
537 }
538 // Transform between differents ellipsoids
539 // using Bursa Wolf parameters.
540 final Matrix step1 = swapAndScaleAxis( sourceCS, GeocentricCoordinateSystem.DEFAULT );
541 final Matrix step2 = getWGS84Parameters( sourceHD );
542 final Matrix step3 = getWGS84Parameters( targetHD );
543 final Matrix step4 = swapAndScaleAxis( GeocentricCoordinateSystem.DEFAULT, targetCS );
544
545 if ( step2 != null && step3 != null )
546 try {
547 // Note: GMatrix.mul(GMatrix) is equivalents to AffineTransform.concatenate(...):
548 // First transform by the supplied transform and then transform the result
549 // by the original transform.
550 step3.invert(); // Invert in place.
551 step4.mul( step3 ); // step4 = step4*step3
552 step4.mul( step2 ); // step4 = step4*step3*step2
553 step4.mul( step1 ); // step4 = step4*step3*step2*step1
554 final MathTransform transform = factory.createAffineTransform( step4 );
555 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION,
556 transform );
557 } catch ( SingularMatrixException exception ) {
558 final CannotCreateTransformException e = new CannotCreateTransformException(
559 sourceCS,
560 targetCS );
561 throw e;
562 }
563 throw new CannotCreateTransformException(
564 Resources.format( ResourceKeys.BURSA_WOLF_PARAMETERS_REQUIRED ) );
565 }
566
567 /**
568 * Creates a transformation between a geographic and a geocentric coordinate systems.
569 * Since the source coordinate systems doesn't have a vertical axis, height above the
570 * ellipsoid is assumed equals to zero everywhere. This method is automatically invoked
571 * by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
572 *
573 * @param sourceCS Input geographic coordinate system.
574 * @param targetCS Output coordinate system.
575 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
576 * @throws CannotCreateTransformException if no transformation path has been found.
577 */
578 protected CoordinateTransformation createTransformationStep(final GeographicCoordinateSystem sourceCS,
579 final GeocentricCoordinateSystem targetCS )
580 throws CannotCreateTransformException {
581 if ( !PrimeMeridian.GREENWICH.equals( targetCS.getPrimeMeridian() ) ) {
582 throw new CannotCreateTransformException( "Rotation of prime meridian not yet implemented" );
583 }
584 final MathTransform step1 =
585 factory.createAffineTransform( swapAndScaleGeoAxis( sourceCS, GeographicCoordinateSystem.WGS84 ) );
586 final MathTransform step2 = getGeocentricTransform( "Ellipsoid_To_Geocentric", 2,
587 sourceCS.getHorizontalDatum().getEllipsoid() );
588 return createFromMathTransform( sourceCS, targetCS, TransformType.CONVERSION,
589 factory.createConcatenatedTransform( step1, step2 ) );
590 }
591
592 /**
593 * Creates a transformation between a projected and a geocentric coordinate systems.
594 * This method is automatically invoked by {@link #createFromCoordinateSystems
595 * createFromCoordinateSystems(...)}. This method doesn't need to be public since
596 * its decomposition in two step should be general enough.
597 *
598 * @param sourceCS Input projected coordinate system.
599 * @param targetCS Output coordinate system.
600 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
601 * @throws CannotCreateTransformException if no transformation path has been found.
602 */
603 private CoordinateTransformation createTransformationStep(
604 final ProjectedCoordinateSystem sourceCS,
605 final GeocentricCoordinateSystem targetCS )
606 throws CannotCreateTransformException {
607 final GeographicCoordinateSystem sourceGCS = sourceCS.getGeographicCoordinateSystem();
608 final CoordinateTransformation step1 = createTransformationStep( sourceCS, sourceGCS );
609 final CoordinateTransformation step2 = createTransformationStep( sourceGCS, targetCS );
610 return concatenate( step1, step2 );
611 }
612
613 /**
614 * Creates a transformation between a compound and a geocentric coordinate systems.
615 * The compound coordinate system <strong>must</strong> be an aggregate of the two
616 * following coordinate systems:
617 *
618 * <ul>
619 * <li>The head must be an instance of {@link HorizontalCoordinateSystem}</li>
620 * <li>The tail must be an instance of {@link VerticalCoordinateSystem}</li>
621 * </ul>
622 *
623 * This method is automatically invoked by {@link #createFromCoordinateSystems
624 * createFromCoordinateSystems(...)}.
625 *
626 * @param sourceCS Input compound coordinate system.
627 * @param targetCS Output coordinate system.
628 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
629 * @throws CannotCreateTransformException if no transformation path has been found.
630 */
631 protected CoordinateTransformation createTransformationStep(
632 final CompoundCoordinateSystem sourceCS,
633 final GeocentricCoordinateSystem targetCS )
634 throws CannotCreateTransformException {
635 /*
636 * Construct a temporary transformation from 'sourceCS' to a standard coordinate
637 * system which is an agregate of a geographic and a vertical coordinate system.
638 * The horizontal datum is preserved, but other properties (vertical datum, axis
639 * order, units, prime meridian...) are "standardized".
640 */
641 final String name = getTemporaryName();
642 final HorizontalDatum datum = OpenGIS.getHorizontalDatum( sourceCS );
643 final CompoundCoordinateSystem stdCS =
644 new CompoundCoordinateSystem( name, new GeographicCoordinateSystem( name, datum ),
645 new VerticalCoordinateSystem( name, VerticalDatum.ELLIPSOIDAL ) );
646 final CoordinateTransformation step1 = createTransformationStep( sourceCS, stdCS );
647 /*
648 * Construct the next steps: conversion to the geocentric coordinate
649 * system, and then conversion to the requested coordinate system.
650 * In summary, the 3 conversions steps are:
651 *
652 * sourceCS --> standardized geographic CS --> standardized geocentric CS --> targetCS
653 */
654 final MathTransform transform = getGeocentricTransform( "Ellipsoid_To_Geocentric", 3,
655 datum.getEllipsoid() );
656 final CoordinateTransformation step2 = createFromMathTransform( stdCS,
657 GeocentricCoordinateSystem.DEFAULT,
658 TransformType.CONVERSION,
659 transform );
660 final CoordinateTransformation step3 = createTransformationStep( GeocentricCoordinateSystem.DEFAULT,
661 targetCS );
662 return concatenate( step1, step2, step3 );
663 }
664
665 /**
666 * Creates a transformation between two compound coordinate systems. This
667 * method is automatically invoked by {@link #createFromCoordinateSystems
668 * createFromCoordinateSystems(...)}.
669 *
670 * @param sourceCS Input coordinate system.
671 * @param targetCS Output coordinate system.
672 * @return A coordinate transformation from <code>sourceCS</code> to <code>targetCS</code>.
673 * @throws CannotCreateTransformException if no transformation path has been found.
674 */
675 protected CoordinateTransformation createTransformationStep(
676 final CompoundCoordinateSystem sourceCS,
677 final CompoundCoordinateSystem targetCS )
678 throws CannotCreateTransformException {
679 final CoordinateSystem headSourceCS = sourceCS.getHeadCS();
680 final CoordinateSystem tailSourceCS = sourceCS.getTailCS();
681 final CoordinateSystem headTargetCS = targetCS.getHeadCS();
682 final CoordinateSystem tailTargetCS = targetCS.getTailCS();
683 if ( tailSourceCS.equivalents( tailTargetCS ) ) {
684 final CoordinateTransformation tr = createFromCoordinateSystems( headSourceCS,
685 headTargetCS );
686 final MathTransform transform = factory.createPassThroughTransform(
687 0,
688 tr.getMathTransform(),
689 tailSourceCS.getDimension() );
690 return createFromMathTransform( sourceCS, targetCS, tr.getTransformType(), transform );
691 }
692 if ( headSourceCS.equivalents( headTargetCS ) ) {
693 final CoordinateTransformation tr = createFromCoordinateSystems( tailSourceCS,
694 tailTargetCS );
695 final MathTransform transform = factory.createPassThroughTransform(
696 headSourceCS.getDimension(),
697 tr.getMathTransform(),
698 0 );
699 return createFromMathTransform( sourceCS, targetCS, tr.getTransformType(), transform );
700 }
701 // TODO: implement others CompoundCoordinateSystem cases.
702 // We could do it in a more general way be creating
703 // and using a 'CompoundTransform' class instead of
704 // of 'PassThroughTransform'. PassThroughTransform
705 // is really a special case of a 'CompoundTransform'
706 // where the head transform is the identity transform.
707 throw new CannotCreateTransformException( sourceCS, targetCS );
708 }
709
710 ////////////////////////////////////////////////////////////////////////
711 //////////// ////////////
712 //////////// HELPER METHODS (private) ////////////
713 //////////// ////////////
714 ////////////////////////////////////////////////////////////////////////
715
716 /**
717 * Returns the WGS84 parameters as an affine transform,
718 * or <code>null</code> if not available.
719 */
720 private static Matrix getWGS84Parameters( final HorizontalDatum datum ) {
721 final WGS84ConversionInfo info = datum.getWGS84Parameters();
722 if ( info != null ) {
723 return info.getAffineTransform();
724 }
725 if ( Ellipsoid.WGS84.equals( datum.getEllipsoid() ) ) {
726 return new Matrix( 4 ); // Identity matrix
727 }
728 return null;
729 }
730
731 /**
732 * Returns a transform between a geocentric
733 * coordinate system and an ellipsoid.
734 *
735 * @param classification either "Ellipsoid_To_Geocentric" or "Geocentric_To_Ellipsoid".
736 * @param dimGeoCS Dimension of the geographic coordinate system (2 or 3).
737 * @param ellipsoid The ellipsoid.
738 * @return The transformation.
739 */
740 private MathTransform getGeocentricTransform( final String classification, final int dimGeoCS,
741 final Ellipsoid ellipsoid ) {
742 final Unit unit = ellipsoid.getAxisUnit();
743 ParameterList param = factory.getMathTransformProvider( classification ).getParameterList();
744 param = param.setParameter( "semi_major", Unit.METRE.convert( ellipsoid.getSemiMajorAxis(),
745 unit ) );
746 param = param.setParameter( "semi_minor", Unit.METRE.convert( ellipsoid.getSemiMinorAxis(),
747 unit ) );
748 try {
749 param = param.setParameter( "dim_geoCS", dimGeoCS );
750 } catch ( IllegalArgumentException exception ) {
751 // The "dim_geoCS" is a custom argument needed for our SEAGIS
752 // implementation. It is not part of OpenGIS's specification.
753 // But if the required dimension is not 3, we can't finish
754 // the operation (TODO: What should we do? Open question...)
755 if ( dimGeoCS != 3 ) {
756 throw exception;
757 }
758 }
759 return factory.createParameterizedTransform( classification, param );
760 }
761
762 /**
763 * Concatenate two transformation steps.
764 *
765 * @param step1 The first step, or <code>null</code> for the identity transform.
766 * @param step2 The second step, or <code>null</code> for the identity transform.
767 * @return A concatenated transform, or <code>null</code> if all arguments was nul.
768 */
769 private CoordinateTransformation concatenate( final CoordinateTransformation step1,
770 final CoordinateTransformation step2 ) {
771 if ( step1 == null )
772 return step2;
773 if ( step2 == null )
774 return step1;
775 if ( !step1.getTargetCS().equivalents( step2.getSourceCS() ) ) {
776 throw new IllegalArgumentException( String.valueOf( step1 ) );
777 }
778 final MathTransform step = factory.createConcatenatedTransform( step1.getMathTransform(),
779 step2.getMathTransform() );
780 final TransformType type = step1.getTransformType().concatenate( step2.getTransformType() );
781
782 return createFromMathTransform( step1.getSourceCS(), step2.getTargetCS(), type, step );
783 }
784
785 /**
786 * Concatenate three transformation steps.
787 *
788 * @param step1 The first step, or <code>null</code> for the identity transform.
789 * @param step2 The second step, or <code>null</code> for the identity transform.
790 * @param step3 The third step, or <code>null</code> for the identity transform.
791 * @return A concatenated transform, or <code>null</code> if all arguments was nul.
792 */
793 private CoordinateTransformation concatenate( final CoordinateTransformation step1,
794 final CoordinateTransformation step2,
795 final CoordinateTransformation step3 ) {
796 if ( step1 == null )
797 return concatenate( step2, step3 );
798 if ( step2 == null )
799 return concatenate( step1, step3 );
800 if ( step3 == null ) {
801 return concatenate( step1, step2 );
802 }
803 if ( !step1.getTargetCS().equivalents( step2.getSourceCS() ) ) {
804 // Current message is for debugging purpose only.
805 throw new IllegalArgumentException( String.valueOf( step1 ) );
806 }
807 if ( !step2.getTargetCS().equivalents( step3.getSourceCS() ) ) {
808 // Current message is for debugging purpose only.
809 throw new IllegalArgumentException( String.valueOf( step3 ) );
810 }
811 final MathTransform step = factory.createConcatenatedTransform(
812 step1.getMathTransform(),
813 factory.createConcatenatedTransform(
814 step2.getMathTransform(),
815 step3.getMathTransform() ) );
816 final TransformType type = step1.getTransformType().concatenate(
817 step2.getTransformType().concatenate(
818 step3.getTransformType() ) );
819 return createFromMathTransform( step1.getSourceCS(), step3.getTargetCS(), type, step );
820 }
821
822 /**
823 * Create a coordinate transform from a math transform.
824 * If the specified math transform is already a coordinate transform, and if source
825 * and target coordinate systems match, then <code>transform</code> is returned with
826 * no change. Otherwise, a new coordinate transform is created.
827 *
828 * @param sourceCS The source coordinate system.
829 * @param targetCS The destination coordinate system.
830 * @param type The transform type.
831 * @param transform The math transform.
832 * @return A coordinate transform using the specified math transform.
833 */
834 private static CoordinateTransformation createFromMathTransform(
835 final CoordinateSystem sourceCS,
836 final CoordinateSystem targetCS,
837 final TransformType type,
838 final MathTransform transform ) {
839 if ( transform instanceof CoordinateTransformation ) {
840 final CoordinateTransformation ct = (CoordinateTransformation) transform;
841 if ( Utilities.equals( ct.getSourceCS(), sourceCS )
842 && Utilities.equals( ct.getTargetCS(), targetCS ) ) {
843 return ct;
844 }
845 }
846 return (CoordinateTransformation) MathTransformFactory.pool.intern( new CoordinateTransformation(
847 null,
848 sourceCS,
849 targetCS,
850 type,
851 transform ) );
852 }
853
854 /**
855 * Returns an affine transform between two coordinate systems. Only units and
856 * axis order (e.g. transforming from (NORTH,WEST) to (EAST,NORTH)) are taken
857 * in account. Other attributes (especially the datum) must be checked before
858 * invoking this method.
859 *
860 * @param sourceCS The source coordinate system. If <code>null</code>, then
861 * (x,y,z,t) axis order is assumed.
862 * @param targetCS The target coordinate system. If <code>null</code>, then
863 * (x,y,z,t) axis order is assumed.
864 */
865 private Matrix swapAndScaleAxis( final CoordinateSystem sourceCS,
866 final CoordinateSystem targetCS )
867 throws CannotCreateTransformException {
868 final AxisOrientation[] sourceAxis = getAxisOrientations( sourceCS, targetCS );
869 final AxisOrientation[] targetAxis = getAxisOrientations( targetCS, sourceCS );
870 final Matrix matrix;
871 try {
872 matrix = Matrix.createAffineTransform( sourceAxis, targetAxis );
873 } catch ( RuntimeException exception ) {
874 final CannotCreateTransformException e = new CannotCreateTransformException( sourceCS,
875 targetCS );
876 throw e;
877 }
878 // Convert units (Optimized case where the conversion
879 // can be applied right into the AffineTransform).
880 final int dimension = matrix.getNumRow() - 1;
881 for ( int i = 0; i < dimension; i++ ) {
882 // TODO: check if units conversion is really linear.
883 final Unit sourceUnit = sourceCS.getUnits( i );
884 final Unit targetUnit = targetCS.getUnits( i );
885 final double offset = targetUnit.convert( 0, sourceUnit );
886 final double scale = targetUnit.convert( 1, sourceUnit ) - offset;
887 matrix.setElement( i, i, scale * matrix.getElement( i, i ) );
888 matrix.setElement( i, dimension, scale * matrix.getElement( i, dimension ) + offset );
889 }
890 return matrix;
891 }
892
893 /**
894 * Returns an affine transform between two geographic coordinate systems. Only
895 * units, axis order (e.g. transforming from (NORTH,WEST) to (EAST,NORTH)) and
896 * prime meridian are taken in account. Other attributes (especially the datum)
897 * must be checked before invoking this method.
898 *
899 * @param sourceCS The source coordinate system.
900 * @param targetCS The target coordinate system.
901 */
902 private Matrix swapAndScaleGeoAxis( final GeographicCoordinateSystem sourceCS,
903 final GeographicCoordinateSystem targetCS )
904 throws CannotCreateTransformException {
905 final Matrix matrix = swapAndScaleAxis( sourceCS, targetCS );
906 for ( int i = targetCS.getDimension(); --i >= 0; ) {
907 // Find longitude ordinate, and apply a rotation if prime meridian are different.
908 final AxisOrientation orientation = targetCS.getAxis( i ).orientation;
909 if ( AxisOrientation.EAST.equals( orientation.absolute() ) ) {
910 final Unit unit = targetCS.getUnits( i );
911 final double sourceLongitude = sourceCS.getPrimeMeridian().getLongitude( unit );
912 final double targetLongitude = targetCS.getPrimeMeridian().getLongitude( unit );
913 final int lastMatrixColumn = matrix.getNumCol() - 1;
914 double rotate = targetLongitude - sourceLongitude;
915 if ( AxisOrientation.WEST.equals( orientation ) )
916 rotate = -rotate;
917 matrix.setElement( i, lastMatrixColumn, matrix.getElement( i, lastMatrixColumn )
918 - rotate );
919 }
920 }
921 return matrix;
922 }
923
924 /**
925 * Returns the axis orientation for the specified coordinate system.
926 * If <code>cs</code> is <code>null</code>, then an array of length
927 * <code>dim.getDimension()</code> is created and filled with
928 * <code>(x,y,z,t)</code> axis orientations.
929 */
930 private static AxisOrientation[] getAxisOrientations( final CoordinateSystem cs,
931 final Dimensioned dim ) {
932 final AxisOrientation[] axis;
933 if ( cs != null ) {
934 axis = new AxisOrientation[cs.getDimension()];
935 for ( int i = 0; i < axis.length; i++ )
936 axis[i] = cs.getAxis( i ).orientation;
937 } else {
938 axis = new AxisOrientation[dim.getDimension()];
939 switch ( axis.length ) {
940 default:
941 for ( int i = axis.length; --i >= 4; )
942 axis[i] = AxisOrientation.OTHER; // fall through
943 case 4:
944 axis[3] = AxisOrientation.FUTURE; // fall through
945 case 3:
946 axis[2] = AxisOrientation.UP; // fall through
947 case 2:
948 axis[1] = AxisOrientation.NORTH; // fall through
949 case 1:
950 axis[0] = AxisOrientation.EAST; // fall through
951 case 0:
952 break;
953 }
954 }
955 return axis;
956 }
957
958 /**
959 * Makes sure that the specified {@link GeographicCoordinateSystem} use standard axis
960 * (longitude and latitude in degrees), Greenwich prime meridian and an ellipsoid
961 * matching projection's parameters. If <code>cs</code> already meets all those
962 * conditions, then it is returned unchanged. Otherwise, a new normalized geographic
963 * coordinate system is created and returned.
964 */
965 private static GeographicCoordinateSystem normalize( GeographicCoordinateSystem cs,
966 final Projection projection ) {
967 HorizontalDatum datum = cs.getHorizontalDatum();
968 Ellipsoid ellipsoid = datum.getEllipsoid();
969 final String name = getTemporaryName();
970 final double semiMajorEll = ellipsoid.getSemiMajorAxis();
971 final double semiMinorEll = ellipsoid.getSemiMinorAxis();
972 final double semiMajorPrj = projection.getValue( "semi_major", semiMajorEll );
973 final double semiMinorPrj = projection.getValue( "semi_minor", semiMinorEll );
974 if ( !hasStandardAxis( cs, Unit.DEGREE )
975 || cs.getPrimeMeridian().getLongitude( Unit.DEGREE ) != 0 ) {
976 cs = null; // Signal that it needs to be reconstructed.
977 }
978 if ( semiMajorEll != semiMajorPrj || semiMinorEll != semiMinorPrj ) {
979 // Those objects are temporary. We assume it is not
980 // a big deal if their name are not very explicit...
981 ellipsoid = new Ellipsoid( name, semiMajorPrj, semiMinorPrj, Unit.METRE );
982 datum = new HorizontalDatum( name, ellipsoid );
983 cs = null; // Signal that it needs to be reconstructed.
984 }
985 if ( cs != null ) {
986 return cs;
987 }
988 return new GeographicCoordinateSystem( name, datum );
989 }
990
991 /**
992 * Makes sure that a {@link ProjectedCoordinateSystem} use standard axis
993 * (x and y in metres) and a normalized {@link GeographicCoordinateSystem}.
994 * If <code>cs</code> already meets all those conditions, then it is
995 * returned unchanged. Otherwise, a new normalized projected coordinate
996 * system is created and returned.
997 */
998 private static ProjectedCoordinateSystem normalize( final ProjectedCoordinateSystem cs ) {
999 final Projection projection = cs.getProjection();
1000 final GeographicCoordinateSystem geoCS = cs.getGeographicCoordinateSystem();
1001 final GeographicCoordinateSystem normalizedGeoCS = normalize( geoCS, projection );
1002
1003 if ( hasStandardAxis( cs, Unit.METRE ) && normalizedGeoCS == geoCS ) {
1004 return cs;
1005 }
1006 return new ProjectedCoordinateSystem( getTemporaryName(), normalizedGeoCS, projection );
1007 }
1008
1009 /**
1010 * Returns <code>true</code> if the specified coordinate
1011 * system use standard axis and standard units.
1012 *
1013 * @param cs The coordinate system to test.
1014 * @paral unit The standard units.
1015 */
1016 private static boolean hasStandardAxis( final HorizontalCoordinateSystem cs, final Unit unit ) {
1017 return cs.getDimension() == 2 /* Just a paranoiac check */
1018 && unit.equals( cs.getUnits( 0 ) ) && unit.equals( cs.getUnits( 1 ) )
1019 && AxisOrientation.EAST.equals( cs.getAxis( 0 ).orientation )
1020 && AxisOrientation.NORTH.equals( cs.getAxis( 1 ).orientation );
1021 }
1022
1023 /**
1024 * Returns a temporary name for generated objects. The first object
1025 * has a name like "Temporary-1", the second is "Temporary-2", etc.
1026 *
1027 */
1028 private static String getTemporaryName() {
1029 return "Temporary-" + ( ++temporaryID );
1030 }
1031
1032 }