001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2005-2006, Geotools Project Managment Committee (PMC)
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: package org.geotools.referencing.operation.projection;
018:
019: // J2SE dependencies
020: import java.util.Collection;
021: import java.awt.geom.Point2D;
022:
023: // OpenGIS dependencies
024: import org.opengis.parameter.ParameterValueGroup;
025: import org.opengis.parameter.ParameterDescriptor;
026: import org.opengis.parameter.ParameterDescriptorGroup;
027: import org.opengis.parameter.ParameterNotFoundException;
028: import org.opengis.referencing.operation.MathTransform;
029: import org.opengis.referencing.ReferenceIdentifier;
030:
031: // Geotools dependencies
032: import org.geotools.math.Complex;
033: import org.geotools.referencing.NamedIdentifier;
034: import org.geotools.metadata.iso.citation.Citations;
035:
036: /**
037: * The NZMG (New Zealand Map Grid) projection.
038: * <p>
039: * This is an implementation of algorithm published by
040: * <a href="http://www.govt.nz/record?recordid=28">Land Information New Zealand</a>.
041: * The algorithm is documented <a href="http://www.linz.govt.nz/rcs/linz/6137/">here</a>.
042: *
043: * @since 2.2
044: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/NewZealandMapGrid.java $
045: * @version $Id: NewZealandMapGrid.java 25262 2007-04-23 21:11:16Z desruisseaux $
046: * @author Justin Deoliveira
047: * @author Martin Desruisseaux
048: *
049: * @todo The algorithm uses complex numbers, which is not very well supported in Java. This
050: * implementation uses {@linkplain Complex} as a support class. Various instances of
051: * {@link Complex} are created once for ever at {@code NewZealandMapGrid} construction
052: * time, in order to avoid creating up to 6 objects for every point to be projected.
053: * The downside is that transformation methods must be synchronized. The cost should
054: * be small for simple applications, but may become important for multi-thread applications.
055: * Furthermore, those fields raise a slight serialization issue.
056: * <p>
057: * The most efficient fix in Java would be to expand inline all {@link Complex} operations
058: * like {@link Complex#add add} (easy), {@link Complex#multiply multiply} (more tedious),
059: * <cite>etc.</cite>, until we get a code using only {@code double} primitives on the stack
060: * and no {@link Complex} objects on the heap (except the {@code A} and {@code B} constants).
061: * But it would make the code significantly more difficult to read.
062: * <p>
063: * An elegant fix would have been "lightweight objects" allocated on the stack (something
064: * similar to {@code struct} in C#), if such thing existed in the Java language.
065: */
066: public class NewZealandMapGrid extends MapProjection {
067: /**
068: * For compatibility with different versions during deserialization.
069: */
070: private static final long serialVersionUID = 8394817836243729133L;
071:
072: /**
073: * Coefficients for forward and inverse projection.
074: */
075: private static final Complex[] A = {
076: new Complex(0.7557853228, 0.0),
077: new Complex(0.249204646, 0.003371507),
078: new Complex(-0.001541739, 0.041058560),
079: new Complex(-0.10162907, 0.01727609),
080: new Complex(-0.26623489, -0.36249218),
081: new Complex(-0.6870983, -1.1651967) };
082:
083: /**
084: * Coefficients for inverse projection.
085: */
086: private static final Complex[] B = {
087: new Complex(1.3231270439, 0.0),
088: new Complex(-0.577245789, -0.007809598),
089: new Complex(0.508307513, -0.112208952),
090: new Complex(-0.15094762, 0.18200602),
091: new Complex(1.01418179, 1.64497696),
092: new Complex(1.9660549, 2.5127645) };
093:
094: /**
095: * Coefficients for inverse projection.
096: */
097: private static final double[] TPHI = new double[] { 1.5627014243,
098: 0.5185406398, -0.03333098, -0.1052906, -0.0368594,
099: 0.007317, 0.01220, 0.00394, -0.0013 };
100:
101: /**
102: * Coefficients for forward projection.
103: */
104: private static final double[] TPSI = new double[] { 0.6399175073,
105: -0.1358797613, 0.063294409, -0.02526853, 0.0117879,
106: -0.0055161, 0.0026906, -0.001333, 0.00067, -0.00034 };
107:
108: /**
109: * A temporary complex number used during transform calculation. Created once for
110: * ever in order to avoid new object creation for every point to be transformed.
111: */
112: private transient final Complex theta = new Complex();
113:
114: /**
115: * An other temporary complex number created once for ever for the same reason than
116: * {@link #theta}. This number is usually equals to some other complex number raised
117: * to some power.
118: */
119: private transient final Complex power = new Complex();
120:
121: /**
122: * An other temporary complex number created once for ever for the same reason than
123: * {@link #theta}.
124: *
125: * @todo Need to reassign those fields on deserialization.
126: */
127: private transient final Complex z = new Complex(),
128: t = new Complex(), num = new Complex(),
129: denom = new Complex();
130:
131: /**
132: * Constructs a new map projection with default parameter values.
133: */
134: protected NewZealandMapGrid() {
135: this ((ParameterValueGroup) Provider.PARAMETERS.createValue());
136: }
137:
138: /**
139: * Constructs a new map projection from the supplied parameters.
140: *
141: * @param parameters The parameter values in standard units.
142: * @throws ParameterNotFoundException if a mandatory parameter is missing.
143: */
144: protected NewZealandMapGrid(final ParameterValueGroup parameters)
145: throws ParameterNotFoundException {
146: super (parameters);
147: }
148:
149: /**
150: * {@inheritDoc}
151: */
152: public ParameterDescriptorGroup getParameterDescriptors() {
153: return Provider.PARAMETERS;
154: }
155:
156: /**
157: * Must be overridden because {@link Provider} uses instances of
158: * {@link ModifiedParameterDescriptor}. This hack was needed because the New Zeland map
159: * projection uses particular default values for parameters like "False Easting", etc.
160: */
161: final boolean isExpectedParameter(final Collection expected,
162: final ParameterDescriptor param) {
163: return ModifiedParameterDescriptor.contains(expected, param);
164: }
165:
166: /**
167: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
168: * (units in radians) and stores the result in {@code ptDst} (linear distance
169: * on a unit sphere).
170: */
171: protected synchronized Point2D transformNormalized(final double x,
172: final double y, final Point2D ptDst)
173: throws ProjectionException {
174: final double dphi = (y - latitudeOfOrigin)
175: * (180 / Math.PI * 3600E-5);
176: double dphi_pow_i = dphi;
177: double dpsi = 0;
178: for (int i = 0; i < TPSI.length; i++) {
179: dpsi += (TPSI[i] * dphi_pow_i);
180: dphi_pow_i *= dphi;
181: }
182: power.real = theta.real = dpsi;
183: power.imag = theta.imag = x;
184: z.multiply(A[0], power);
185: for (int i = 1; i < A.length; i++) {
186: power.multiply(power, theta);
187: z.addMultiply(z, A[i], power);
188: }
189: if (ptDst != null) {
190: ptDst.setLocation(z.imag, z.real);
191: return ptDst;
192: }
193: return new Point2D.Double(z.imag, z.real);
194: }
195:
196: /**
197: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
198: * and stores the result in {@code ptDst}.
199: */
200: protected synchronized Point2D inverseTransformNormalized(
201: final double x, final double y, final Point2D ptDst)
202: throws ProjectionException {
203: power.real = z.real = y;
204: power.imag = z.imag = x;
205: theta.multiply(B[0], z);
206: for (int j = 1; j < B.length; j++) {
207: power.multiply(power, z);
208: theta.addMultiply(theta, B[j], power);
209: }
210: // increasing the number of iterations through this loop decreases
211: // the error in the calculation, but 3 iterations gives 10-3 accuracy
212: for (int j = 0; j < 3; j++) {
213: power.power(theta, 2);
214: num.addMultiply(z, A[1], power);
215: for (int k = 2; k < A.length; k++) {
216: power.multiply(power, theta);
217: t.multiply(A[k], power);
218: t.multiply(t, k);
219: num.add(num, t);
220: }
221:
222: power.real = 1;
223: power.imag = 0;
224: denom.copy(A[0]);
225: for (int k = 1; k < A.length; k++) {
226: power.multiply(power, theta);
227: t.multiply(A[k], power);
228: t.multiply(t, k + 1);
229: denom.add(denom, t);
230: }
231: theta.divide(num, denom);
232: }
233:
234: final double dpsi = theta.real;
235: double dpsi_pow_i = dpsi;
236: double dphi = TPHI[0] * dpsi;
237: for (int i = 1; i < TPHI.length; i++) {
238: dpsi_pow_i *= dpsi;
239: dphi += (TPHI[i] * dpsi_pow_i);
240: }
241:
242: dphi = dphi / (180 / Math.PI * 3600E-5) + latitudeOfOrigin;
243: if (ptDst != null) {
244: ptDst.setLocation(theta.imag, dphi);
245: return ptDst;
246: }
247: return new Point2D.Double(theta.imag, dphi);
248: }
249:
250: //////////////////////////////////////////////////////////////////////////////////////////
251: //////////////////////////////////////////////////////////////////////////////////////////
252: //////// ////////
253: //////// PROVIDERS ////////
254: //////// ////////
255: //////////////////////////////////////////////////////////////////////////////////////////
256: //////////////////////////////////////////////////////////////////////////////////////////
257:
258: /**
259: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
260: * provider} for {@linkplain NewZealandMapGrid New Zealand Map Grid} (EPSG code 27200).
261: *
262: * @since 2.2
263: * @version $Id: NewZealandMapGrid.java 25262 2007-04-23 21:11:16Z desruisseaux $
264: * @author Justin Deoliveira
265: *
266: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
267: */
268: public static class Provider extends AbstractProvider {
269: /**
270: * For compatibility with different versions during deserialization.
271: */
272: private static final long serialVersionUID = -7716733400419275656L;
273:
274: /**
275: * The parameters group.
276: */
277: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
278: new ReferenceIdentifier[] {
279: new NamedIdentifier(Citations.OGC,
280: "New_Zealand_Map_Grid"),
281: new NamedIdentifier(Citations.EPSG,
282: "New Zealand Map Grid"),
283: new NamedIdentifier(Citations.EPSG, "27200") },
284: new ParameterDescriptor[] {
285: new ModifiedParameterDescriptor(SEMI_MAJOR,
286: 6378388.0),
287: new ModifiedParameterDescriptor(SEMI_MINOR,
288: 6378388.0 * (1 - 1 / 297.0)),
289: new ModifiedParameterDescriptor(
290: LATITUDE_OF_ORIGIN, -41.0),
291: new ModifiedParameterDescriptor(
292: CENTRAL_MERIDIAN, 173.0),
293: new ModifiedParameterDescriptor(FALSE_EASTING,
294: 2510000.0),
295: new ModifiedParameterDescriptor(FALSE_NORTHING,
296: 6023150.0) });
297:
298: /**
299: * Constructs a new provider.
300: */
301: public Provider() {
302: super (PARAMETERS);
303: }
304:
305: /**
306: * Creates a transform from the specified group of parameter values. This method doesn't
307: * check for the spherical case, since the New Zealand Map Grid projection is used with
308: * the International 1924 ellipsoid.
309: *
310: * @param parameters The group of parameter values.
311: * @return The created math transform.
312: * @throws ParameterNotFoundException if a required parameter was not found.
313: */
314: public MathTransform createMathTransform(
315: final ParameterValueGroup parameters)
316: throws ParameterNotFoundException {
317: return new NewZealandMapGrid(parameters);
318: }
319: }
320: }
|