001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2007, GeoTools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation;
010: * version 2.1 of the License.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: *
017: * This package contains documentation from OpenGIS specifications.
018: * OpenGIS consortium's work is fully acknowledged here.
019: */
020: package org.geotools.referencing.operation;
021:
022: // J2SE dependencies and extensions
023: import java.util.List;
024: import java.util.Iterator;
025: import java.util.Collection;
026: import java.util.LinkedList;
027: import java.util.logging.Level;
028: import java.util.logging.LogRecord;
029: import javax.units.Unit;
030: import javax.units.SI;
031:
032: // OpenGIS dependencies
033: import org.opengis.parameter.ParameterValue;
034: import org.opengis.parameter.ParameterValueGroup;
035: import org.opengis.parameter.ParameterDescriptor;
036: import org.opengis.parameter.ParameterDescriptorGroup;
037: import org.opengis.parameter.GeneralParameterValue;
038: import org.opengis.parameter.GeneralParameterDescriptor;
039: import org.opengis.referencing.IdentifiedObject;
040: import org.opengis.referencing.crs.ProjectedCRS;
041: import org.opengis.referencing.cs.AxisDirection;
042: import org.opengis.referencing.cs.CoordinateSystem;
043: import org.opengis.referencing.operation.Matrix;
044: import org.opengis.referencing.operation.Conversion;
045: import org.opengis.referencing.operation.MathTransform;
046:
047: // Geotools dependencies
048: import org.geotools.resources.Utilities;
049: import org.geotools.referencing.CRS;
050: import org.geotools.referencing.cs.AbstractCS;
051: import org.geotools.referencing.AbstractIdentifiedObject;
052: import org.geotools.referencing.factory.ReferencingFactory;
053: import org.geotools.referencing.operation.matrix.XMatrix;
054: import org.geotools.referencing.operation.matrix.MatrixFactory;
055: import org.geotools.referencing.operation.transform.AbstractMathTransform;
056: import org.geotools.referencing.operation.transform.ConcatenatedTransform;
057: import org.geotools.referencing.operation.projection.MapProjection; // For javadoc
058: import org.geotools.resources.i18n.LoggingKeys;
059: import org.geotools.resources.i18n.Logging;
060:
061: /**
062: * Returns a conversion from a source to target projected CRS, if this conversion
063: * is representable as an affine transform. More specifically, if all projection
064: * parameters are identical except the following ones:
065: * <P>
066: * <UL>
067: * <LI>{@link MapProjection.AbstractProvider#SCALE_FACTOR scale_factor}</LI>
068: * <LI>{@link MapProjection.AbstractProvider#FALSE_EASTING false_easting}</LI>
069: * <LI>{@link MapProjection.AbstractProvider#FALSE_NORTHING false_northing}</LI>
070: * </UL>
071: * <P>
072: * Then the conversion between two projected CRS can sometime be represented as a linear
073: * conversion. For example if only false easting/northing differ, then the coordinate conversion
074: * is simply a translation.
075: *
076: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/ProjectionAnalyzer.java $
077: * @version $Id: ProjectionAnalyzer.java 24581 2007-02-26 01:36:46Z desruisseaux $
078: * @author Martin Desruisseaux
079: */
080: final class ProjectionAnalyzer {
081: /**
082: * The map projection.
083: */
084: private final Conversion projection;
085:
086: /**
087: * The affine transform applied on geographic coordinates before the projection.
088: * In Geotools {@link MapProjection} implementation, this is the axis swapping and
089: * scaling needed in order to get standard (<var>longitude</var>,<var>latitude</var>)
090: * axis in degrees. Can be {@code null} if none.
091: * <p>
092: * This is not needed for {@code ProjectionAnalyzer} working, but is stored anyway
093: * for debugging purpose.
094: */
095: private final Matrix geographicScale;
096:
097: /**
098: * The affine transform applied on projected coordinates after the projection.
099: * In Geotools {@link MapProjection} implementation, this is the axis swapping
100: * and scaling needed in order to get standard (<var>x</var>,<var>y</var>) axis
101: * in metres. Can be {@code null} if none.
102: */
103: private final Matrix projectedScale;
104:
105: /**
106: * The transform for the map projection alone, without the {@link #geographicScale}
107: * and {@link #projectedScale} parts. In Geotools implementation, it should be an
108: * instance of {@link MapProjection}. May be {@code null} if we can't handle the
109: * {@linkplain #projection}.
110: */
111: private final MathTransform transform;
112:
113: /**
114: * The map projection parameters values, or a copy of them.
115: */
116: private List/*<GeneralParameterValue>*/parameters;
117:
118: /**
119: * Constructs a {@code ProjectionAnalyzer} for the specified projected CRS. This constructor
120: * inspects the {@linkplain ProjectedCRS#getConversionFromBase conversion from base} and splits
121: * {@link ConcatenatedTransform} in their {@link #geographicScale}, {@link #projectedScale}
122: * and {@link #transform} components.
123: */
124: private ProjectionAnalyzer(final ProjectedCRS crs) {
125: Matrix geographicScale = null;
126: Matrix projectedScale = null;
127: projection = crs.getConversionFromBase();
128: MathTransform candidate = projection.getMathTransform();
129: while (candidate instanceof ConcatenatedTransform) {
130: final ConcatenatedTransform ctr = (ConcatenatedTransform) candidate;
131: if (ctr.transform1 instanceof LinearTransform) {
132: if (geographicScale != null) {
133: // Should never happen with ConcatenatedTransform.create(...) implementation.
134: throw new IllegalStateException(String
135: .valueOf(candidate));
136: }
137: geographicScale = ((LinearTransform) ctr.transform1)
138: .getMatrix();
139: candidate = ctr.transform2;
140: continue;
141: }
142: if (ctr.transform2 instanceof LinearTransform) {
143: if (projectedScale != null) {
144: // Should never happen with ConcatenatedTransform.create(...) implementation.
145: throw new IllegalStateException(String
146: .valueOf(candidate));
147: }
148: projectedScale = ((LinearTransform) ctr.transform2)
149: .getMatrix();
150: candidate = ctr.transform1;
151: continue;
152: }
153: // Both transforms are non-linear. We can not handle that.
154: candidate = null;
155: break;
156: }
157: //
158: // TODO: We need to handle PassthroughTransform here in some future version
159: // (when we will want better handling of 3D coordinates).
160: //
161: /*
162: * We should really fetch the parameters from the MathTransform as much as we can, since
163: * this is the most robust source of information (the one which is the most likely to be
164: * an accurate description of the map projection without the above geographic and projected
165: * scale components). However if we are not able to query the math transform, we will query
166: * the Conversion object as a fallback and hope that it describes only the map projection
167: * part, as in Geotools implementation.
168: */
169: ParameterValueGroup group = null;
170: if (candidate instanceof AbstractMathTransform) {
171: group = ((AbstractMathTransform) candidate)
172: .getParameterValues();
173: }
174: if (group == null) {
175: /*
176: * Fallback path only if we don't have a Geotools MapProjection implementation.
177: * NOTE: it is uncertain that we should call 'swapAndScaleAxis'. If the CS has
178: * standard axis, it will not hurt since we should get the identity transform.
179: * If the CS doesn't have standard axis, then 'projectedScale' should be non-
180: * null and 'swapAndScaleAxis' is not needed. But if none of the above hold,
181: * then some axis swapping is probably done straight into the unknown 'transform'
182: * implementation and we need to "guess" what it is. Those rules are somewhat
183: * heuristic; the previous "if" branch for Geotools MapProjection implementations
184: * should be more determinist.
185: */
186: group = projection.getParameterValues();
187: if (projectedScale == null) {
188: final CoordinateSystem cs = crs.getCoordinateSystem();
189: projectedScale = AbstractCS.swapAndScaleAxis(AbstractCS
190: .standard(cs), cs);
191: }
192: }
193: if (group != null) {
194: parameters = group.values();
195: }
196: this .geographicScale = geographicScale;
197: this .projectedScale = projectedScale;
198: this .transform = candidate;
199: }
200:
201: /**
202: * Returns the {@linkplain #transform} parameter descriptor, or {@code null} if none.
203: */
204: private ParameterDescriptorGroup getTransformDescriptor() {
205: return (transform instanceof AbstractMathTransform) ? ((AbstractMathTransform) transform)
206: .getParameterDescriptors()
207: : null;
208: }
209:
210: /**
211: * Returns the affine transform applied after the <em>normalized</em> projection in order to
212: * get the same projection than {@link #transform}. The normalized projection is a imaginary
213: * transform (we don't have a {@link MathTransform} instance for it, but we don't need) with
214: * {@code "scale factor"} == 1, {@code "false easting"} == 0 and {@code "false northing"} == 0.
215: * In other words, this method extracts the above-cited parameters in an affine transform.
216: * <p>
217: * As a side effect, this method removes from the {@linkplain #parameters} list
218: * all the above-cited ones parameters.
219: *
220: * @return The affine transform.
221: */
222: private XMatrix normalizedToProjection() {
223: parameters = new LinkedList(parameters); // Keep the original list unchanged.
224: /*
225: * Creates a matrix which will conceptually stands between the normalized transform and
226: * the 'projectedScale' transform. The matrix dimensions are selected accordingly using
227: * a robust code when possible, but the result should be a 3x3 matrix most of the time.
228: */
229: final int sourceDim = (transform != null) ? transform
230: .getTargetDimensions() : 2;
231: final int targetDim = (projectedScale != null) ? projectedScale
232: .getNumCol() - 1 : sourceDim;
233: final XMatrix matrix = MatrixFactory.create(targetDim + 1,
234: sourceDim + 1);
235: /*
236: * Search for "scale factor", "false easting" and "false northing" parameters.
237: * All parameters found are removed from the list. Note: we assume that linear
238: * units in the "normalized projection" are metres, as specified in the legacy
239: * OGC 01-009 specification, and that unit conversions (if needed) are applied
240: * by 'projectedScale'. However there is no way I can see to ensure that since
241: * it is really a matter of how the map projection is implemented (for example
242: * the unit conversion factor could be merged with the "scale_factor" -- not a
243: * clean approach and I do not recommand, but some could do in order to save a
244: * few multiplications).
245: *
246: * We need "false easting" and "false northing" offsets in either user's unit or
247: * in metres, depending if the unit conversions are applied in 'transform' or in
248: * 'projectedScale' respectively. We assume the later, which stands for Geotools
249: * implementation and is closer to OGC 01-009 spirit. But we will log a warning
250: * in case of doubt.
251: */
252: Unit unit = null;
253: String warning = null;
254: for (final Iterator it = parameters.iterator(); it.hasNext();) {
255: final GeneralParameterValue parameter = (GeneralParameterValue) it
256: .next();
257: if (parameter instanceof ParameterValue) {
258: final ParameterValue value = (ParameterValue) parameter;
259: // TODO: remove the cast below when we will be allowed to compile for J2SE 1.5.
260: final ParameterDescriptor descriptor = (ParameterDescriptor) value
261: .getDescriptor();
262: if (Number.class.isAssignableFrom(descriptor
263: .getValueClass())) {
264: if (nameMatches(descriptor, "scale_factor")) {
265: final double scale = value.doubleValue();
266: for (int i = Math.min(sourceDim, targetDim); --i >= 0;) {
267: matrix.setElement(i, i, matrix.getElement(
268: i, i)
269: * scale);
270: }
271: } else {
272: final int d;
273: if (nameMatches(descriptor, "false_easting")) {
274: d = 0;
275: } else if (nameMatches(descriptor,
276: "false_northing")) {
277: d = 1;
278: } else {
279: continue;
280: }
281: final double offset = value
282: .doubleValue(SI.METER);
283: if (!Double.isNaN(offset)
284: && offset != value.doubleValue()) {
285: // See the above comment about units. The above check could have been
286: // replaced by "if (!SI.METER.equals(unit))", but the above avoid the
287: // warning in the very common case where 'offset == 0'.
288: unit = value.getUnit();
289: warning = descriptor.getName().getCode();
290: }
291: matrix.setElement(d, sourceDim, matrix
292: .getElement(d, sourceDim)
293: + offset);
294: }
295: it.remove();
296: }
297: }
298: }
299: if (warning != null) {
300: final LogRecord record = Logging.format(Level.WARNING,
301: LoggingKeys.APPLIED_UNIT_CONVERSION_$3, warning,
302: unit, SI.METER);
303: record.setSourceClassName(getClass().getName());
304: record.setSourceMethodName("createLinearConversion"); // This is the public method.
305: ReferencingFactory.LOGGER.log(record);
306: }
307: return matrix;
308: }
309:
310: /**
311: * Checks if the parameter in the two specified list contains the same values.
312: * The order parameter order is irrelevant. The common parameters are removed
313: * from both lists.
314: */
315: private static boolean parameterValuesEqual(
316: final List/*<GeneralParameterValue>*/source,
317: final List/*<GeneralParameterValue>*/target,
318: final double errorTolerance) {
319: search: for (final Iterator targetIter = target.iterator(); targetIter
320: .hasNext();) {
321: final GeneralParameterValue targetPrm = (GeneralParameterValue) targetIter
322: .next();
323: for (final Iterator sourceIter = source.iterator(); sourceIter
324: .hasNext();) {
325: final GeneralParameterValue sourcePrm = (GeneralParameterValue) sourceIter
326: .next();
327: if (!nameMatches(sourcePrm.getDescriptor(), targetPrm
328: .getDescriptor())) {
329: continue;
330: }
331: if (sourcePrm instanceof ParameterValue
332: && targetPrm instanceof ParameterValue) {
333: final ParameterValue sourceValue = (ParameterValue) sourcePrm;
334: final ParameterValue targetValue = (ParameterValue) targetPrm;
335: // TODO: Remove the cast when we will be allowed to compile for J2SE 1.5.
336: if (Number.class
337: .isAssignableFrom(((ParameterDescriptor) targetPrm
338: .getDescriptor()).getValueClass())) {
339: final double sourceNum, targetNum;
340: final Unit unit = targetValue.getUnit();
341: if (unit != null) {
342: sourceNum = sourceValue.doubleValue(unit);
343: targetNum = targetValue.doubleValue(unit);
344: } else {
345: sourceNum = sourceValue.doubleValue();
346: targetNum = targetValue.doubleValue();
347: }
348: double error = (targetNum - sourceNum);
349: if (targetNum != 0)
350: error /= targetNum;
351: if (!(Math.abs(error) <= errorTolerance)) { // '!' for trapping NaN
352: return false;
353: }
354: } else {
355: // The parameter do not hold a numerical value. It may be a
356: // String, etc. Use the generic Object.equals(Object) method.
357: if (!Utilities.equals(sourceValue.getValue(),
358: targetValue.getValue())) {
359: return false;
360: }
361: }
362: } else {
363: // The GeneralParameter is not a ParameterValue instance. It is probably a
364: // ParameterValueGroup. Compare all child elements without processing them.
365: if (!Utilities.equals(targetPrm, sourcePrm)) {
366: return false;
367: }
368: }
369: // End of processing a pair of matching parameters. The values are equal or
370: // were one of the special cases processed above. Continue with a new pair.
371: sourceIter.remove();
372: targetIter.remove();
373: continue search;
374: }
375: // End of iteration in the 'source' parameters. If we reach this point, then we
376: // have found a target parameter without matching source parameter. We consider
377: // the two projections as different kind.
378: return false;
379: }
380: // End of iteration in the 'target' parameters, which should now be empty.
381: // Check if there is any unmatched parameter left in the supplied list.
382: assert target.isEmpty();
383: return source.isEmpty();
384: }
385:
386: /**
387: * Applies {@code normalizedToProjection} first, then {@link #projectedScale}.
388: */
389: private XMatrix applyProjectedScale(
390: final XMatrix normalizedToProjection) {
391: if (projectedScale == null) {
392: return normalizedToProjection;
393: }
394: final XMatrix scale = MatrixFactory.create(projectedScale);
395: scale.multiply(normalizedToProjection);
396: return scale;
397: }
398:
399: /**
400: * Returns a conversion from a source to target projected CRS, if this conversion
401: * is representable as an affine transform. If no linear conversion has been found
402: * between the two CRS, then this method returns {@code null}.
403: *
404: * @param sourceCRS The source coordinate reference system.
405: * @param targetCRS The target coordinate reference system.
406: * @param errorTolerance Relative error tolerance for considering two parameter values as
407: * equal. This is usually a small number like {@code 1E-10}.
408: * @return The conversion from {@code sourceCRS} to {@code targetCRS} as an
409: * affine transform, or {@code null} if no linear transform has been found.
410: */
411: public static Matrix createLinearConversion(
412: final ProjectedCRS sourceCRS, final ProjectedCRS targetCRS,
413: final double errorTolerance) {
414: /*
415: * Checks if the datum are the same. To be stricter, we could compare the 'baseCRS'
416: * instead. But this is not always needed. For example we don't really care if the
417: * underlying geographic CRS use different axis order or units. What matter are the
418: * axis order and units of the projected CRS.
419: *
420: * Actually, checking for 'baseCRS' causes an infinite loop (until StackOverflowError)
421: * in CoordinateOperationFactory, because it prevents this method to recognize that the
422: * transform between two projected CRS is the identity transform even if their underlying
423: * geographic CRS use different axis order.
424: */
425: if (!CRS.equalsIgnoreMetadata(sourceCRS.getDatum(), targetCRS
426: .getDatum())) {
427: return null;
428: }
429: final ProjectionAnalyzer source = new ProjectionAnalyzer(
430: sourceCRS);
431: final ProjectionAnalyzer target = new ProjectionAnalyzer(
432: targetCRS);
433: if (!nameMatches(source.projection.getMethod(),
434: target.projection.getMethod())) {
435: /*
436: * In theory, we can not find a linear conversion if the operation method is
437: * not the same. In practice, it still hapen in some occasions. For example
438: * "Transverse Mercator" and "Transverse Mercator (South Oriented)" are two
439: * distinct operation methods in EPSG point of view, but in Geotools the South
440: * Oriented case is implemented as a "Transverse Mercator" concatenated with
441: * an affine transform performing the axis flip, which is a linear conversion.
442: *
443: * We may be tempted to compare the 'source.transform' and 'target.transform'
444: * implementation classes, but this is not robust enough. For example it is
445: * possible to implement the "Oblique Mercator" and "Hotine Oblique Mercator"
446: * projections with a single class. But both cases can have identical user-
447: * supplied parameters and still be different projections (they differ in the
448: * origin of their grid coordinates).
449: *
450: * As a compromise, we compare the method name declared by the math transform,
451: * in addition of the method name declared by the conversion (the check above).
452: */
453: final ParameterDescriptorGroup sourceDsc = source
454: .getTransformDescriptor();
455: final ParameterDescriptorGroup targetDsc = source
456: .getTransformDescriptor();
457: if (sourceDsc == null || targetDsc == null
458: || !nameMatches(sourceDsc, targetDsc)) {
459: return null;
460: }
461: }
462: /*
463: * Extracts the "scale_factor", "false_easting" and "false_northing" parameters
464: * as affine transforms. All remaining parameters must be identical.
465: */
466: if (source.parameters == null || target.parameters == null) {
467: return null;
468: }
469: XMatrix sourceScale = source.normalizedToProjection();
470: XMatrix targetScale = target.normalizedToProjection();
471: if (!parameterValuesEqual(source.parameters, target.parameters,
472: errorTolerance)) {
473: return null;
474: }
475: /*
476: * Creates the matrix (including axis order changes and unit conversions),
477: * and apply the scale and translation inferred from the "false_easting"
478: * parameter and its friends. We perform the conversion in three conceptual
479: * steps (in the end, everything is bundle in a single matrix):
480: *
481: * 1) remove the old false northing/easting
482: * 2) apply the scales
483: * 3) add the new false northing/easting
484: */
485: targetScale = target.applyProjectedScale(targetScale);
486: sourceScale = source.applyProjectedScale(sourceScale);
487: sourceScale.invert();
488: targetScale.multiply(sourceScale);
489: if (targetScale.isIdentity(errorTolerance)) {
490: targetScale.setIdentity();
491: }
492: return targetScale;
493: }
494:
495: /**
496: * Temporary convenience method to be deleted when we will be allow
497: * to use static imports in J2SE 1.5.
498: */
499: private static boolean nameMatches(final IdentifiedObject object,
500: final String name) {
501: return AbstractIdentifiedObject.nameMatches(object, name);
502: }
503:
504: /**
505: * Temporary convenience method to be deleted when we will be allow
506: * to use static imports in J2SE 1.5.
507: */
508: private static boolean nameMatches(final IdentifiedObject o1,
509: final IdentifiedObject o2) {
510: return AbstractIdentifiedObject.nameMatches(o1, o2);
511: }
512: }
|