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.CompoundCRS;
051: import org.deegree.crs.coordinatesystems.CoordinateSystem;
052: import org.deegree.crs.coordinatesystems.GeocentricCRS;
053: import org.deegree.crs.coordinatesystems.GeographicCRS;
054: import org.deegree.crs.coordinatesystems.ProjectedCRS;
055: import org.deegree.crs.exceptions.TransformationException;
056: import org.deegree.crs.projections.ProjectionUtils;
057: import org.deegree.crs.utilities.Matrix;
058: import org.deegree.framework.log.ILogger;
059: import org.deegree.framework.log.LoggerFactory;
060: import org.deegree.i18n.Messages;
061:
062: /**
063: * The <code>GeocentricTransform</code> class is the central access point for all transformations between different
064: * crs's.
065: * <p>
066: * It creates a transformation chain for two given CoordinateSystems by considering their type. For example the
067: * Transformation chain from EPSG:31466 ( a projected crs with underlying geographic crs epsg:4314 using the DHDN datum
068: * and the TransverseMercator Projection) to EPSG:28992 (another projected crs with underlying geographic crs epsg:4289
069: * using the 'new Amersfoort Datum' and the StereographicAzimuthal Projection) would result in following Transformation
070: * Chain:
071: * <ol>
072: * <li>Inverse projection - thus getting the coordinates in lat/lon for geographic crs epsg:4314</li>
073: * <li>Geodetic transformation - thus getting x-y-z coordinates for geographic crs epsg:4314</li>
074: * <li>WGS84 transformation -thus getting the x-y-z coordinates for the WGS84 datum </li>
075: * <li>Inverse WGS84 transformation -thus getting the x-y-z coordinates for the geodetic from epsg:4289</li>
076: * <li>Inverse geodetic - thus getting the lat/lon for epsg:4289</li>
077: * <li>projection - getting the coordinates (in meters) for epsg:28992</li>
078: * </ol>
079: *
080: * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
081: *
082: * @author last edited by: $Author:$
083: *
084: * @version $Revision:$, $Date:$
085: *
086: */
087: public class TransformationFactory {
088: private static ILogger LOG = LoggerFactory
089: .getLogger(TransformationFactory.class);
090:
091: /**
092: * The default coordinate transformation factory. Will be constructed only when first needed.
093: */
094: private static TransformationFactory DEFAULT_INSTANCE = null;
095:
096: private TransformationFactory() {
097: // nottin
098: }
099:
100: /**
101: * @return the default coordinate transformation factory.
102: */
103: public static synchronized TransformationFactory getInstance() {
104: if (DEFAULT_INSTANCE == null) {
105: // DEFAULT_INSTANCE = new CoordinateTransformationFactory(
106: // MathTransformFactory.getDefault() );
107: DEFAULT_INSTANCE = new TransformationFactory();
108: }
109: return DEFAULT_INSTANCE;
110: }
111:
112: /**
113: * Creates a transformation between two coordinate systems. This method will examine the coordinate systems in order
114: * to construct a transformation between them. This method may fail if no path between the coordinate systems is
115: * found.
116: *
117: * @param sourceCRS
118: * Input coordinate system.
119: * @param targetCRS
120: * Output coordinate system.
121: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
122: * @throws TransformationException
123: * @throws TransformationException
124: * if no transformation path has been found.
125: * @throws IllegalArgumentException
126: * if the sourceCRS or targetCRS are <code>null</code>.
127: *
128: */
129: public CRSTransformation createFromCoordinateSystems(
130: final CoordinateSystem sourceCRS,
131: final CoordinateSystem targetCRS)
132: throws TransformationException, IllegalArgumentException {
133: if (sourceCRS == null) {
134: throw new IllegalArgumentException(
135: "The source CRS may not be null");
136: }
137: if (targetCRS == null) {
138: throw new IllegalArgumentException(
139: "The target CRS may not be null");
140: }
141: if ((sourceCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS
142: && sourceCRS.getType() != CoordinateSystem.COMPOUND_CRS
143: && sourceCRS.getType() != CoordinateSystem.PROJECTED_CRS && sourceCRS
144: .getType() != CoordinateSystem.GEOCENTRIC_CRS)
145: || (targetCRS.getType() != CoordinateSystem.GEOGRAPHIC_CRS
146: && targetCRS.getType() != CoordinateSystem.COMPOUND_CRS
147: && targetCRS.getType() != CoordinateSystem.PROJECTED_CRS && targetCRS
148: .getType() != CoordinateSystem.GEOCENTRIC_CRS)) {
149: throw new TransformationException(sourceCRS, targetCRS,
150: "Either the target crs type or the source crs type was unknown");
151: }
152: final boolean sourceIsCompound = (sourceCRS.getType() == CoordinateSystem.COMPOUND_CRS);
153: final boolean targetIsCompound = (targetCRS.getType() == CoordinateSystem.COMPOUND_CRS);
154: CoordinateSystem tmpSourceCRS = sourceCRS;
155: CoordinateSystem tmpTargetCRS = targetCRS;
156: if (sourceIsCompound) {
157: tmpSourceCRS = ((CompoundCRS) sourceCRS).getUnderlyingCRS();
158: }
159: if (targetIsCompound) {
160: tmpTargetCRS = ((CompoundCRS) targetCRS).getUnderlyingCRS();
161: }
162: if (tmpSourceCRS.equals(tmpTargetCRS)) {
163: LOG
164: .logDebug("Source crs and target crs are equal, no transformation needed (returning identity matrix).");
165: final Matrix matrix = new Matrix(
166: sourceCRS.getDimension() + 1);
167: matrix.setIdentity();
168: return createAffineTransform(sourceCRS, targetCRS, matrix);
169: }
170:
171: CRSTransformation result = null;
172: if (tmpSourceCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS) {
173: /**
174: * Geographic --> Geographic, Projected or Geocentric
175: */
176: final GeographicCRS source = (GeographicCRS) tmpSourceCRS;
177: if (tmpTargetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS) {
178: result = createTransformationStep(source,
179: (GeographicCRS) tmpTargetCRS, sourceIsCompound,
180: targetIsCompound);
181: }
182: if (tmpTargetCRS.getType() == CoordinateSystem.PROJECTED_CRS) {
183: result = createTransformationStep(source,
184: (ProjectedCRS) tmpTargetCRS, sourceIsCompound,
185: targetIsCompound);
186: }
187: if (tmpTargetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS) {
188: result = createTransformationStep(source,
189: (GeocentricCRS) tmpTargetCRS, sourceIsCompound);
190: }
191: } else if (tmpSourceCRS.getType() == CoordinateSystem.PROJECTED_CRS) {
192: /**
193: * Projected --> Projected, Geographic or Geocentric
194: */
195: final ProjectedCRS source = (ProjectedCRS) tmpSourceCRS;
196: if (tmpTargetCRS.getType() == CoordinateSystem.PROJECTED_CRS) {
197: result = createTransformationStep(source,
198: (ProjectedCRS) tmpTargetCRS, sourceIsCompound,
199: targetIsCompound);
200: }
201: if (tmpTargetCRS.getType() == CoordinateSystem.GEOGRAPHIC_CRS) {
202: result = createTransformationStep(source,
203: (GeographicCRS) tmpTargetCRS, sourceIsCompound,
204: targetIsCompound);
205: }
206: if (tmpTargetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS) {
207: result = createTransformationStep(source,
208: (GeocentricCRS) tmpTargetCRS, sourceIsCompound);
209: }
210: } else if (tmpSourceCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS) {
211: /**
212: * Geocentric --> Geocentric
213: */
214: final GeocentricCRS source = (GeocentricCRS) sourceCRS;
215: if (tmpTargetCRS.getType() == CoordinateSystem.GEOCENTRIC_CRS) {
216: result = createTransformationStep(source,
217: (GeocentricCRS) tmpTargetCRS);
218: }
219: }
220: if (result == null) {
221: LOG
222: .logDebug("The resulting transformation was null, returning an identity matrix.");
223: final Matrix matrix = new Matrix(
224: sourceCRS.getDimension() + 1);
225: matrix.setIdentity();
226: result = createAffineTransform(sourceCRS, targetCRS, matrix);
227: } else {
228: result = concatenate(createAffineTransform(sourceCRS,
229: sourceCRS, toStandardizedValues(sourceCRS, false)),
230: result, createAffineTransform(targetCRS, targetCRS,
231: toStandardizedValues(targetCRS, true)));
232: }
233: if (LOG.getLevel() == ILogger.LOG_DEBUG) {
234: StringBuilder output = new StringBuilder(
235: "The resulting transformation chain: \n");
236: output = result.getTransformationPath(output);
237: LOG.logDebug(output.toString());
238: }
239: return result;
240:
241: }
242:
243: /**
244: * Creates a matrix, with which incoming values will be transformed to a standardized form. This means, to radians
245: * and meters.
246: *
247: * @param sourceCRS
248: * to create the matrix for.
249: * @param invert
250: * the values. Using the inverted scale, i.e. going from standard to crs specific.
251: * @return the standardized matrix.
252: * @throws TransformationException
253: * if the unit of one of the axis could not be transformed to one of the base units.
254: */
255: private Matrix toStandardizedValues(CoordinateSystem sourceCRS,
256: boolean invert) throws TransformationException {
257: final int dim = sourceCRS.getDimension();
258: Matrix result = null;
259: Axis[] allAxis = sourceCRS.getAxis();
260: for (int i = 0; i < allAxis.length; ++i) {
261: Axis targetAxis = allAxis[i];
262: if (targetAxis != null) {
263: Unit targetUnit = targetAxis.getUnits();
264: if (!Unit.RADIAN.equals(targetUnit)
265: || !Unit.METRE.equals(targetUnit)) {
266: if (!(targetUnit.canConvert(Unit.RADIAN) || targetUnit
267: .canConvert(Unit.METRE))) {
268: throw new TransformationException(
269: Messages
270: .getMessage(
271: "CRS_TRANSFORMATION_NO_APLLICABLE_UNIT",
272: targetUnit));
273: }
274: // lazy instantiation
275: if (result == null) {
276: result = new Matrix(dim + 1);
277: result.setIdentity();
278: }
279: result.setElement(i, i, invert ? 1d / targetUnit
280: .getScale() : targetUnit.getScale());
281: }
282: }
283: }
284: return result;
285: }
286:
287: /**
288: * Creates a transformation between two geographic coordinate systems. This method is automatically invoked by
289: * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust
290: * axis order and orientation (e.g. transforming from <code>(NORTH,WEST)</code> to <code>(EAST,NORTH)</code>),
291: * performs units conversion and apply Bursa Wolf transformation if needed.
292: *
293: * @param sourceCRS
294: * Input coordinate system.
295: * @param targetCRS
296: * Output coordinate system.
297: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
298: * @throws TransformationException
299: * if no transformation path has been found.
300: */
301: private CRSTransformation createTransformationStep(
302: final GeographicCRS sourceCRS,
303: final GeographicCRS targetCRS,
304: final boolean sourceIsCompound,
305: final boolean targetIsCompound)
306: throws TransformationException {
307: final GeodeticDatum sourceDatum = sourceCRS.getGeodeticDatum();
308: final GeodeticDatum targetDatum = targetCRS.getGeodeticDatum();
309: if (sourceDatum.equals(targetDatum)) {
310: LOG.logDebug("The datums of geographic (source): "
311: + sourceCRS.getIdAndName()
312: + " equals geographic (target): "
313: + targetCRS.getIdAndName() + " returning null");
314: return null;
315: }
316: LOG
317: .logDebug("Creating geographic ->geographic transformation: from (source): "
318: + sourceCRS.getIdAndName()
319: + " to(target): "
320: + targetCRS.getIdAndName());
321: final Ellipsoid sourceEllipsoid = sourceDatum.getEllipsoid();
322: final Ellipsoid targetEllipsoid = targetDatum.getEllipsoid();
323: if (sourceEllipsoid != null
324: && !sourceEllipsoid.equals(targetEllipsoid)
325: && (sourceDatum.getWGS84Conversion().hasValues() // check if a conversion needs to take place
326: || targetDatum.getWGS84Conversion().hasValues())) {
327: /*
328: * If the two geographic coordinate systems use different ellipsoid, convert from the source to target
329: * ellipsoid through the geocentric coordinate system.
330: */
331: final GeocentricCRS sourceGCS = new GeocentricCRS(
332: sourceDatum, sourceCRS.getIdentifier(), sourceCRS
333: .getName()
334: + "_Geocentric");
335: final GeocentricCRS targetGCS = new GeocentricCRS(
336: targetDatum, targetCRS.getIdentifier(), targetCRS
337: .getName()
338: + "_Geocentric");
339: CRSTransformation step1 = createTransformationStep(
340: sourceCRS, sourceGCS, sourceIsCompound);
341: LOG.logDebug("Step 1 from geographic to geographic:\n"
342: + step1);
343: CRSTransformation step2 = createTransformationStep(
344: sourceGCS, targetGCS);
345: LOG.logDebug("Step 2 from geographic to geographic:\n"
346: + step2);
347: CRSTransformation step3 = createTransformationStep(
348: targetCRS, targetGCS, targetIsCompound);
349: LOG.logDebug("Step 3 from geographic to geographic:\n"
350: + step3);
351: if (step3 != null) {
352: step3.inverse();// call inverseTransform from step 3.
353: }
354: return concatenate(step1, step2, step3);
355: }
356:
357: /*
358: * Swap axis order, and rotate the longitude coordinate if prime meridians are different.
359: */
360: final Matrix matrix = swapAndScaleGeoAxis(sourceCRS, targetCRS);
361: return createAffineTransform(sourceCRS, targetCRS, matrix);
362: }
363:
364: /**
365: * Creates a transformation between a geographic and a projected coordinate systems. This method is automatically
366: * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
367: *
368: * @param sourceCRS
369: * Input coordinate system.
370: * @param targetCRS
371: * Output coordinate system.
372: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
373: * @throws TransformationException
374: * if no transformation path has been found.
375: */
376: private CRSTransformation createTransformationStep(
377: final GeographicCRS sourceCRS,
378: final ProjectedCRS targetCRS,
379: final boolean sourceIsCompound,
380: final boolean targetIsCompound)
381: throws TransformationException {
382:
383: LOG
384: .logDebug("Creating geographic->projected transformation: from (source): "
385: + sourceCRS.getIdAndName()
386: + " to(target): "
387: + targetCRS.getIdAndName());
388: final GeographicCRS stepGeoCS = targetCRS.getGeographicCRS();
389:
390: // should perform a datum shift
391: final CRSTransformation step1 = createTransformationStep(
392: sourceCRS, stepGeoCS, sourceIsCompound,
393: targetIsCompound);
394: final CRSTransformation step2 = new ProjectionTransform(
395: targetCRS);
396: return concatenate(step1, step2);
397: }
398:
399: /**
400: * Creates a transformation between a geographic and a geocentric coordinate systems. Since the source coordinate
401: * systems doesn't have a vertical axis, height above the ellipsoid is assumed equals to zero everywhere. This
402: * method is automatically invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}.
403: *
404: * @param sourceCRS
405: * Input geographic coordinate system.
406: * @param targetCRS
407: * Output coordinate system.
408: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
409: * @throws TransformationException
410: * if no transformation path has been found.
411: */
412: private CRSTransformation createTransformationStep(
413: final GeographicCRS sourceCRS,
414: final GeocentricCRS targetCRS, boolean hasHeight)
415: throws TransformationException {
416: LOG
417: .logDebug("Creating geographic -> geocentric (helmert) transformation: from (source): "
418: + sourceCRS.getIdAndName()
419: + " to(target): "
420: + targetCRS.getIdAndName());
421: if (!PrimeMeridian.GREENWICH.equals(targetCRS
422: .getGeodeticDatum().getPrimeMeridian())) {
423: throw new TransformationException(
424: "The rotation from "
425: + targetCRS.getGeodeticDatum()
426: .getPrimeMeridian().getIdAndName()
427: + " to Greenwich prime meridian is not yet implemented");
428: }
429: final CRSTransformation step1 = createAffineTransform(
430: sourceCRS, targetCRS, swapAndScaleGeoAxis(sourceCRS,
431: GeographicCRS.WGS84));
432:
433: // TODO maybe something for a pool of transformations?
434: final CRSTransformation step2 = new GeocentricTransform(
435: sourceCRS, targetCRS, hasHeight);
436: return createConcatenatedTransform(step1, step2);
437: }
438:
439: /**
440: * Creates a transformation between two projected coordinate systems. This method is automatically invoked by
441: * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust
442: * axis order and orientation. It also performs units conversion if it is the only extra change needed. Otherwise,
443: * it performs three steps:
444: *
445: * <ol>
446: * <li>Unproject <code>sourceCRS</code>.</li>
447: * <li>Transform from <code>sourceCRS.geographicCS</code> to <code>targetCRS.geographicCS</code>.</li>
448: * <li>Project <code>targetCRS</code>.</li>
449: * </ol>
450: *
451: * @param sourceCRS
452: * Input coordinate system.
453: * @param targetCRS
454: * Output coordinate system.
455: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
456: * @throws TransformationException
457: * if no transformation path has been found.
458: */
459: private CRSTransformation createTransformationStep(
460: final ProjectedCRS sourceCRS, final ProjectedCRS targetCRS,
461: final boolean sourceIsCompound,
462: final boolean targetIsCompound)
463: throws TransformationException {
464: LOG
465: .logDebug("Creating projected -> projected transformation: from (source): "
466: + sourceCRS.getIdAndName()
467: + " to(target): "
468: + targetCRS.getIdAndName());
469: if (sourceCRS.getProjection().equals(targetCRS.getProjection())) {
470: return null;
471: }
472: final GeographicCRS sourceGeo = sourceCRS.getGeographicCRS();
473: final GeographicCRS targetGeo = targetCRS.getGeographicCRS();
474: final CRSTransformation step1 = createTransformationStep(
475: sourceCRS, sourceGeo, sourceIsCompound,
476: sourceIsCompound);
477: LOG.logDebug("Step 1: " + step1);
478: final CRSTransformation step2 = createTransformationStep(
479: sourceGeo, targetGeo, sourceIsCompound,
480: targetIsCompound);
481: LOG.logDebug("Step 2: " + step2);
482: final CRSTransformation step3 = createTransformationStep(
483: targetGeo, targetCRS, targetIsCompound,
484: targetIsCompound);
485: LOG.logDebug("Step 3: " + step3);
486: return concatenate(step1, step2, step3);
487: }
488:
489: /**
490: * Creates a transformation between a projected and a geocentric coordinate systems. This method is automatically
491: * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. This method doesn't need to be
492: * public since its decomposition in two step should be general enough.
493: *
494: * @param sourceCRS
495: * Input projected coordinate system.
496: * @param targetCRS
497: * Output coordinate system.
498: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
499: * @throws TransformationException
500: * if no transformation path has been found.
501: */
502: private CRSTransformation createTransformationStep(
503: final ProjectedCRS sourceCRS,
504: final GeocentricCRS targetCRS,
505: final boolean sourceIsCompound)
506: throws TransformationException {
507: LOG
508: .logDebug("Creating projected -> geocentric transformation: from (source): "
509: + sourceCRS.getIdAndName()
510: + " to(target): "
511: + targetCRS.getIdAndName());
512: final GeographicCRS sourceGCS = sourceCRS.getGeographicCRS();
513:
514: final CRSTransformation step1 = createTransformationStep(
515: sourceCRS, sourceGCS, sourceIsCompound,
516: sourceIsCompound);
517: final CRSTransformation step2 = createTransformationStep(
518: sourceGCS, targetCRS, sourceIsCompound);
519: return concatenate(step1, step2);
520: }
521:
522: /**
523: * Creates a transformation between a projected and a geographic coordinate systems. This method is automatically
524: * invoked by {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation
525: * returns <code>{@link #createTransformationStep(GeographicCRS, ProjectedCRS, boolean, boolean)}
526: * createTransformationStep}(targetCRS, sourceCRS) inverse)</code>.
527: *
528: * @param sourceCRS
529: * Input coordinate system.
530: * @param targetCRS
531: * Output coordinate system.
532: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code> or
533: * <code>null</code> if {@link ProjectedCRS#getGeographicCRS()}.equals targetCRS.
534: * @throws TransformationException
535: * if no transformation path has been found.
536: */
537: private CRSTransformation createTransformationStep(
538: final ProjectedCRS sourceCRS,
539: final GeographicCRS targetCRS,
540: final boolean sourceIsCompound,
541: final boolean targetIsCompound)
542: throws TransformationException {
543: LOG
544: .logDebug("Creating projected->geographic transformation: from (source): "
545: + sourceCRS.getIdAndName()
546: + " to(target): "
547: + targetCRS.getIdAndName());
548: CRSTransformation result = createTransformationStep(targetCRS,
549: sourceCRS, targetIsCompound, sourceIsCompound);
550: result.inverse();
551: return result;
552:
553: }
554:
555: /**
556: * Creates a transformation between two geocentric coordinate systems. This method is automatically invoked by
557: * {@link #createFromCoordinateSystems createFromCoordinateSystems(...)}. The default implementation can adjust for
558: * axis order and orientation, adjust for prime meridian, performs units conversion and apply Bursa Wolf
559: * transformation if needed.
560: *
561: * @param sourceCRS
562: * Input coordinate system.
563: * @param targetCRS
564: * Output coordinate system.
565: * @return A coordinate transformation from <code>sourceCRS</code> to <code>targetCRS</code>.
566: * @throws TransformationException
567: * if no transformation path has been found.
568: */
569: private CRSTransformation createTransformationStep(
570: final GeocentricCRS sourceCRS, final GeocentricCRS targetCRS)
571: throws TransformationException {
572: LOG
573: .logDebug("Creating geocentric->geocetric transformation: from (source): "
574: + sourceCRS.getIdAndName()
575: + " to(target): "
576: + targetCRS.getIdAndName());
577: final GeodeticDatum sourceHD = sourceCRS.getGeodeticDatum();
578: final GeodeticDatum targetHD = targetCRS.getGeodeticDatum();
579: if (sourceHD.equals(targetHD)) {
580: return null;
581: // final Matrix matrix = swapAndScaleAxis( sourceCRS, targetCRS );
582: // return createAffineTransform( sourceCRS, targetCRS, matrix );
583: }
584: if (!PrimeMeridian.GREENWICH
585: .equals(sourceHD.getPrimeMeridian())
586: || !PrimeMeridian.GREENWICH.equals(targetHD
587: .getPrimeMeridian())) {
588: throw new TransformationException(
589: "Rotation of prime meridian not yet implemented");
590: }
591: // Transform between different ellipsoids using Bursa Wolf parameters.
592: Matrix tmp = swapAndScaleAxis(sourceCRS, GeocentricCRS.WGS84);
593: Matrix4d step1 = null;
594: if (tmp != null) {
595: step1 = new Matrix4d();
596: tmp.get(step1);
597: }
598: final Matrix4d step2 = getWGS84Parameters(sourceHD);
599: final Matrix4d step3 = getWGS84Parameters(targetHD);
600: tmp = swapAndScaleAxis(GeocentricCRS.WGS84, targetCRS);
601: Matrix4d step4 = null;
602: if (tmp != null) {
603: step4 = new Matrix4d();
604: tmp.get(step4);
605: }
606: if (step1 == null && step2 == null && step3 == null
607: && step4 == null) {
608: LOG
609: .logDebug("The given geocentric crs's do not need a helmert transformation (but they are not equal), returning identity");
610: step4 = new Matrix4d();
611: step4.setIdentity();
612: } else {
613: LOG.logDebug("step1 matrix: \n " + step1);
614: LOG.logDebug("step2 matrix: \n " + step2);
615: LOG.logDebug("step3 matrix: \n " + step3);
616: LOG.logDebug("step4 matrix: \n " + step4);
617: if (step3 != null) {
618: step3.invert(); // Invert in place.
619: LOG.logDebug("step3 inverted matrix: \n " + step3);
620: }
621: if (step4 != null) {
622: if (step3 != null) {
623: step4.mul(step3); // step4 = step4*step3
624: LOG.logDebug("step4 (after mul3): \n " + step4);
625: }
626: if (step2 != null) {
627: step4.mul(step2); // step4 = step4*step3*step2
628: LOG.logDebug("step4 (after mul2): \n " + step4);
629: }
630: if (step1 != null) {
631: step4.mul(step1); // step4 = step4*step3*step2*step1
632: }
633: } else if (step3 != null) {
634: step4 = step3;
635: if (step2 != null) {
636: step4.mul(step2); // step4 = step3*step2*step1
637: LOG.logDebug("step4 (after mul2): \n " + step4);
638: }
639: if (step1 != null) {
640: step4.mul(step1); // step4 = step3*step2*step1
641: }
642: } else if (step2 != null) {
643: step4 = step2;
644: if (step1 != null) {
645: step4.mul(step1); // step4 = step2*step1
646: }
647: } else {
648: step4 = step1;
649: }
650: }
651:
652: LOG.logDebug("The resulting geo-geo matrix: from( "
653: + sourceCRS.getIdentifier() + ") to("
654: + targetCRS.getIdentifier() + ")\n " + step4);
655: return new MatrixTransform(sourceCRS, targetCRS, step4,
656: "Helmert-Transformation");
657: }
658:
659: /**
660: * Concatenates two existing transforms.
661: *
662: * @param tr1
663: * The first transform to apply to points.
664: * @param tr2
665: * The second transform to apply to points.
666: * @return The concatenated transform.
667: *
668: */
669: private CRSTransformation createConcatenatedTransform(
670: CRSTransformation tr1, CRSTransformation tr2) {
671: if (tr1 == null) {
672: return tr2;
673: }
674: if (tr2 == null) {
675: return tr1;
676: }
677: // if one of the two is an identity transformation, just return the other.
678: if (tr1.isIdentity()) {
679: return tr2;
680: }
681: if (tr2.isIdentity()) {
682: return tr1;
683: }
684:
685: /*
686: * If one transform is the inverse of the other, just return an identitiy transform.
687: */
688: if (tr1.areInverse(tr2)) {
689: LOG
690: .logDebug("Transformation1 and Transformation2 are inverse operations, returning null");
691: return null;
692: }
693:
694: /*
695: * If both transforms use matrix, then we can create a single transform using the concatened matrix.
696: */
697: if (tr1 instanceof MatrixTransform
698: && tr2 instanceof MatrixTransform) {
699: GMatrix m1 = ((MatrixTransform) tr1).getMatrix();
700: GMatrix m2 = ((MatrixTransform) tr2).getMatrix();
701: if (m1 == null) {
702: if (m2 == null) {
703: // both matrices are null, just return the identiy matrix.
704: return new MatrixTransform(tr1.getSourceCRS(), tr1
705: .getTargetCRS(), new GMatrix(tr2
706: .getTargetDimension() + 1, tr1
707: .getSourceDimension() + 1));
708: }
709: return tr2;
710: } else if (m2 == null) {
711: return tr1;
712: }
713: m2.mul(m1);
714: return new MatrixTransform(tr1.getSourceCRS(), tr2
715: .getTargetCRS(), m2);
716: }
717:
718: /*
719: * If one or both math transform are instance of {@link ConcatenatedTransform}, then concatenate <code>tr1</code>
720: * or <code>tr2</code> with one of step transforms.
721: */
722: if (tr1 instanceof ConcatenatedTransform) {
723: final ConcatenatedTransform ctr = (ConcatenatedTransform) tr1;
724: tr1 = ctr.getFirstTransform();
725: tr2 = createConcatenatedTransform(ctr.getSecondTransform(),
726: tr2);
727: }
728: if (tr2 instanceof ConcatenatedTransform) {
729: final ConcatenatedTransform ctr = (ConcatenatedTransform) tr2;
730: tr1 = createConcatenatedTransform(tr1, ctr
731: .getFirstTransform());
732: tr2 = ctr.getSecondTransform();
733: }
734:
735: return new ConcatenatedTransform(tr1, tr2);
736: }
737:
738: /**
739: * Creates an affine transform from a matrix.
740: *
741: * @param matrix
742: * The matrix used to define the affine transform.
743: * @return The affine transform.
744: * @throws TransformationException
745: * if the matrix is not affine.
746: *
747: */
748: private MatrixTransform createAffineTransform(
749: CoordinateSystem source, CoordinateSystem target,
750: final Matrix matrix) throws TransformationException {
751: if (matrix == null) {
752: return null;
753: }
754: if (matrix.isAffine()) {// Affine transform are square.
755: if (matrix.getNumRow() == 3 && matrix.getNumCol() == 3
756: && !matrix.isIdentity()) {
757: LOG
758: .logDebug("Creating affineTransform of incoming matrix:\n"
759: + matrix);
760: Matrix3d tmp = matrix.toAffineTransform();
761: LOG.logDebug("Resulting AffineTransform is:\n" + tmp);
762: return new MatrixTransform(source, target, tmp);
763: }
764: return new MatrixTransform(source, target, matrix);
765: }
766: throw new TransformationException(
767: "Given matrix is not affine, cannot continue");
768: }
769:
770: /**
771: * Returns the WGS84 parameters as an affine transform, or <code>null</code> if not available.
772: */
773: private Matrix4d getWGS84Parameters(final GeodeticDatum datum) {
774: final WGS84ConversionInfo info = datum.getWGS84Conversion();
775: if (info != null && info.hasValues()) {
776: return info.getAsAffineTransform();
777: }
778: return null;
779: }
780:
781: /**
782: * Concatenate two transformation steps.
783: *
784: * @param step1
785: * The first step, or <code>null</code> for the identity transform.
786: * @param step2
787: * The second step, or <code>null</code> for the identity transform.
788: * @return A concatenated transform, or <code>null</code> if all arguments was nul.
789: */
790: private CRSTransformation concatenate(
791: final CRSTransformation step1, final CRSTransformation step2) {
792: if (step1 == null)
793: return step2;
794: if (step2 == null)
795: return step1;
796: return createConcatenatedTransform(step1, step2);
797: }
798:
799: /**
800: * Concatenate three transformation steps.
801: *
802: * @param step1
803: * The first step, or <code>null</code> for the identity transform.
804: * @param step2
805: * The second step, or <code>null</code> for the identity transform.
806: * @param step3
807: * The third step, or <code>null</code> for the identity transform.
808: * @return A concatenated transform, or <code>null</code> if all arguments were <code>null</code>.
809: */
810: private CRSTransformation concatenate(
811: final CRSTransformation step1,
812: final CRSTransformation step2, final CRSTransformation step3) {
813: if (step1 == null) {
814: return concatenate(step2, step3);
815: }
816: if (step2 == null) {
817: return concatenate(step1, step3);
818: }
819: if (step3 == null) {
820: return concatenate(step1, step2);
821: }
822: return createConcatenatedTransform(step1,
823: createConcatenatedTransform(step2, step3));
824: }
825:
826: /**
827: * Returns an affine transform between two coordinate systems. Only units and axis order (e.g. transforming from
828: * (NORTH,WEST) to (EAST,NORTH)) are taken in account. Other attributes (especially the datum) must be checked
829: * before invoking this method.
830: *
831: * @param sourceCRS
832: * The source coordinate system.
833: * @param targetCRS
834: * The target coordinate system.
835: */
836: private Matrix swapAndScaleAxis(final CoordinateSystem sourceCRS,
837: final CoordinateSystem targetCRS)
838: throws TransformationException {
839: LOG.logDebug("Creating swap matrix from: "
840: + sourceCRS.getIdAndName() + " to: "
841: + targetCRS.getIdAndName());
842: final Matrix matrix;
843: try {
844: matrix = new Matrix(sourceCRS.getAxis(), targetCRS
845: .getAxis());
846: } catch (RuntimeException e) {
847: throw new TransformationException(sourceCRS, targetCRS, e
848: .getMessage());
849: }
850: return matrix.isIdentity() ? null : matrix;
851: }
852:
853: /**
854: * Returns an affine transform between two geographic coordinate systems. Only units, axis order (e.g. transforming
855: * from (NORTH,WEST) to (EAST,NORTH)) and prime meridian are taken in account. Other attributes (especially the
856: * datum) must be checked before invoking this method.
857: *
858: * @param sourceCRS
859: * The source coordinate system.
860: * @param targetCRS
861: * The target coordinate system.
862: */
863: private Matrix swapAndScaleGeoAxis(final GeographicCRS sourceCRS,
864: final GeographicCRS targetCRS)
865: throws TransformationException {
866: LOG.logDebug("Creating geo swap matrix from: "
867: + sourceCRS.getIdAndName() + " to: "
868: + targetCRS.getIdAndName());
869: final Matrix matrix = swapAndScaleAxis(sourceCRS, targetCRS);
870: if (matrix != null) {
871: Axis[] targetAxis = targetCRS.getAxis();
872: final int lastMatrixColumn = matrix.getNumCol() - 1;
873: for (int i = 1; --i >= 0;) {
874: // Find longitude, and apply a translation if prime meridians are different.
875: final int orientation = targetAxis[i].getOrientation();
876: if (Axis.AO_EAST == Math.abs(orientation)) {
877: final Unit unit = targetAxis[i].getUnits();
878: final double sourceLongitude = sourceCRS
879: .getGeodeticDatum().getPrimeMeridian()
880: .getLongitude(unit);
881: final double targetLongitude = targetCRS
882: .getGeodeticDatum().getPrimeMeridian()
883: .getLongitude(unit);
884: if (Math.abs(sourceLongitude - targetLongitude) > ProjectionUtils.EPS11) {
885: double translation = targetLongitude
886: - sourceLongitude;
887: if (Axis.AO_WEST == orientation) {
888: translation = -translation;
889: }
890: // add the translation to the matrix translate element of this axis
891: matrix.setElement(i, lastMatrixColumn, matrix
892: .getElement(i, lastMatrixColumn)
893: - translation);
894: }
895: }
896: }
897: }
898: return matrix;
899: }
900: }
|