001: /*
002: * Geotools2 - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2002-2005, Geotools Project Managment Committee (PMC)
005: *
006: * This library is free software; you can redistribute it and/or
007: * modify it under the terms of the GNU Lesser General Public
008: * License as published by the Free Software Foundation;
009: * version 2.1 of the License.
010: *
011: * This library is distributed in the hope that it will be useful,
012: * but WITHOUT ANY WARRANTY; without even the implied warranty of
013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
014: * Lesser General Public License for more details.
015: */
016: package org.geotools.referencing.operation.projection;
017:
018: import java.awt.geom.Point2D;
019: import java.util.Collection;
020: import javax.units.NonSI;
021: import javax.units.SI;
022: import javax.units.Unit;
023:
024: import org.opengis.parameter.ParameterDescriptor;
025: import org.opengis.parameter.ParameterDescriptorGroup;
026: import org.opengis.parameter.ParameterNotFoundException;
027: import org.opengis.parameter.ParameterValueGroup;
028: import org.opengis.referencing.operation.ConicProjection;
029: import org.opengis.referencing.operation.MathTransform;
030:
031: import org.geotools.metadata.iso.citation.Citations;
032: import org.geotools.referencing.NamedIdentifier;
033: import org.geotools.resources.XMath;
034: import org.geotools.resources.i18n.ErrorKeys;
035: import org.geotools.resources.i18n.Errors;
036:
037: /**
038: * Krovak Oblique Conformal Conic projection (EPSG code 9819). This projection is
039: * used in the Czech Republic and Slovakia under the name 'Krovak' projection. The
040: * geographic coordinates on the ellipsoid are first reduced to conformal
041: * coordinates on the conformal (Gaussian) sphere. These spherical coordinates
042: * are then projected onto the oblique cone and converted to grid coordinates.
043: * The pseudo standard parallel is defined on the conformal sphere after its
044: * rotation, to obtain the oblique aspect of the projection. It is then the
045: * parallel on this sphere at which the map projection is true to scale; on
046: * the ellipsoid it maps as a complex curve.
047: *
048: * <p>The compulsory parameters are just the ellipsoid characteristics.
049: * All other parameters are optional and have defaults to match the common
050: * usage with Krovak projection.</p>
051: *
052: * <p>In general the axis of Krovak projection are defined as westing and
053: * southing (not easting and northing) and they are also reverted, so if the
054: * value of projected coordinates should (and in <var>y</var>, <var>x</var>
055: * order in Krovak) be positive the 'Axis' parameter for projection should be
056: * defined explicitly like this (in wkt):</p>
057: *
058: * <pre>PROJCS["S-JTSK (Ferro) / Krovak",
059: * .
060: * .
061: * .
062: *
063: * PROJECTION["Krovak"]
064: * PARAMETER["semi_major", 6377397.155],
065: * PARAMETER["semi_minor", 6356078.963],
066: * UNIT["meter",1.0],
067: * AXIS["x", WEST],
068: * AXIS["y", SOUTH]]
069: * </pre>Axis in Krovak:
070: * <pre>
071: * y<------------------+
072: * |
073: * Czech. Rep. |
074: * |
075: * x
076: * </pre>
077: * <p>By default, the axis are 'easting, northing' so the values of projected coordinates
078: * are negative (and in <var>y</var>, <var>x</var> order in Krovak - it is cold
079: * Krovak GIS version).</p>
080: *
081: * <p><strong>References:</strong>
082: * <ul>
083: * <li>Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
084: * Relevant files is: {@code PJ_krovak.c}</li>
085: * <li>"Coordinate Conversions and Transformations including Formulas" available at, <A
086: * HREF="http://www.remotesensing.org/geotiff/proj_list/guid7.html">http://www.remotesensing.org/geotiff/proj_list/guid7.html</A></li>
087: * </ul>
088: * </p>
089: *
090: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/krovak.html">Krovak on
091: * RemoteSensing.org </A>
092: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/guid7.html">Krovak on "Coordinate
093: * Conversions and Transformations including Formulas"</A>
094: * @see <A HREF="http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34e2.html">Krovak on POSC</A>
095: *
096: * @since 2.4
097: * @version $Id: Krovak.java 24563 2007-02-23 00:20:43Z desruisseaux $
098: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/Krovak.java $
099: * @author Jan Jezek
100: * @author Martin Desruisseaux
101: */
102: public class Krovak extends MapProjection {
103: /**
104: * Maximum number of iterations for iterative computations.
105: */
106: private static final int MAXIMUM_ITERATIONS = 15;
107:
108: /**
109: * When to stop the iteration.
110: */
111: private static final double ITERATION_TOLERANCE = 1E-11;
112:
113: /**
114: * Azimuth of the centre line passing through the centre of the projection.
115: * This is equals to the co-latitude of the cone axis at point of intersection
116: * with the ellipsoid.
117: */
118: protected final double azimuth;
119:
120: /**
121: * Latitude of pseudo standard parallel.
122: */
123: protected final double pseudoStandardParallel;
124:
125: /**
126: * Useful variables calculated from parameters defined by user.
127: */
128: private final double sinAzim, cosAzim, n, tanS2, alfa, hae, k1, ka,
129: ro0, rop;
130:
131: /**
132: * Useful constant - 45° in radians.
133: */
134: private static final double s45 = 0.785398163397448;
135:
136: /**
137: * Constructs a new map projection from the supplied parameters.
138: *
139: * @param parameters The parameter values in standard units.
140: * @throws ParameterNotFoundException if a mandatory parameter is missing.
141: */
142: protected Krovak(final ParameterValueGroup parameters)
143: throws ParameterNotFoundException {
144: super (parameters);
145: final Collection expected = getParameterDescriptors()
146: .descriptors();
147: // Fetch parameters from user input.
148: latitudeOfOrigin = doubleValue(expected,
149: Provider.LATITUDE_OF_CENTER, parameters);
150: centralMeridian = doubleValue(expected,
151: Provider.LONGITUDE_OF_CENTER, parameters);
152: azimuth = doubleValue(expected, Provider.AZIMUTH, parameters);
153: pseudoStandardParallel = doubleValue(expected,
154: Provider.PSEUDO_STANDARD_PARALLEL, parameters);
155: scaleFactor = doubleValue(expected, Provider.SCALE_FACTOR,
156: parameters);
157: ensureLatitudeInRange(Provider.LATITUDE_OF_CENTER,
158: latitudeOfOrigin, false);
159: ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTER,
160: centralMeridian, false);
161:
162: // Calculates useful constants.
163: sinAzim = Math.sin(azimuth);
164: cosAzim = Math.cos(azimuth);
165: n = Math.sin(pseudoStandardParallel);
166: tanS2 = Math.tan(pseudoStandardParallel / 2 + s45);
167:
168: final double sinLat, cosLat, cosL2, u0;
169: sinLat = Math.sin(latitudeOfOrigin);
170: cosLat = Math.cos(latitudeOfOrigin);
171: cosL2 = cosLat * cosLat;
172: alfa = Math
173: .sqrt(1 + ((excentricitySquared * (cosL2 * cosL2)) / (1 - excentricitySquared)));
174: hae = alfa * excentricity / 2;
175: u0 = Math.asin(sinLat / alfa);
176:
177: final double g, esl;
178: esl = excentricity * sinLat;
179: g = Math.pow((1 - esl) / (1 + esl), (alfa * excentricity) / 2);
180: k1 = Math.pow(Math.tan(latitudeOfOrigin / 2 + s45), alfa) * g
181: / Math.tan(u0 / 2 + s45);
182: ka = Math.pow(1 / k1, -1 / alfa);
183:
184: final double radius;
185: radius = Math.sqrt(1 - excentricitySquared)
186: / (1 - (excentricitySquared * (sinLat * sinLat)));
187:
188: ro0 = scaleFactor * radius / Math.tan(pseudoStandardParallel);
189: rop = ro0 * Math.pow(tanS2, n);
190: }
191:
192: /**
193: * {@inheritDoc}
194: */
195: public ParameterDescriptorGroup getParameterDescriptors() {
196: return Provider.PARAMETERS;
197: }
198:
199: /**
200: * {@inheritDoc}
201: */
202: public ParameterValueGroup getParameterValues() {
203: final Collection expected = getParameterDescriptors()
204: .descriptors();
205: final ParameterValueGroup values = super .getParameterValues();
206: set(expected, Provider.LATITUDE_OF_CENTER, values,
207: latitudeOfOrigin);
208: set(expected, Provider.LONGITUDE_OF_CENTER, values,
209: centralMeridian);
210: set(expected, Provider.AZIMUTH, values, azimuth);
211: set(expected, Provider.PSEUDO_STANDARD_PARALLEL, values,
212: pseudoStandardParallel);
213: set(expected, Provider.SCALE_FACTOR, values, scaleFactor);
214: return values;
215: }
216:
217: /**
218: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
219: * (units in radians) and stores the result in {@code ptDst} (linear distance
220: * on a unit sphere).
221: */
222: protected Point2D transformNormalized(final double lambda,
223: final double phi, Point2D ptDst) throws ProjectionException {
224: final double esp = excentricity * Math.sin(phi);
225: final double gfi = Math.pow(((1. - esp) / (1. + esp)), hae);
226: final double u = 2 * (Math.atan(Math.pow(Math
227: .tan(phi / 2 + s45), alfa)
228: / k1 * gfi) - s45);
229: final double deltav = -lambda * alfa;
230: final double cosU = Math.cos(u);
231: final double s = Math.asin((cosAzim * Math.sin(u))
232: + (sinAzim * cosU * Math.cos(deltav)));
233: final double d = Math.asin(cosU * Math.sin(deltav)
234: / Math.cos(s));
235: final double eps = n * d;
236: final double ro = rop / Math.pow(Math.tan(s / 2 + s45), n);
237:
238: /* x and y are reverted */
239: final double y = -(ro * Math.cos(eps));
240: final double x = -(ro * Math.sin(eps));
241:
242: if (ptDst != null) {
243: ptDst.setLocation(x, y);
244: return ptDst;
245: }
246: return new Point2D.Double(x, y);
247: }
248:
249: /**
250: * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
251: * and stores the result in {@code ptDst}.
252: */
253: protected Point2D inverseTransformNormalized(final double x,
254: final double y, Point2D ptDst) throws ProjectionException {
255: // x -> southing, y -> westing
256: final double ro = XMath.hypot(x, y);
257: final double eps = Math.atan2(-x, -y);
258: final double d = eps / n;
259: final double s = 2 * (Math.atan(Math.pow(ro0 / ro, 1 / n)
260: * tanS2) - s45);
261: final double cs = Math.cos(s);
262: final double u = Math.asin((cosAzim * Math.sin(s))
263: - (sinAzim * cs * Math.cos(d)));
264: final double kau = ka
265: * Math.pow(Math.tan((u / 2.) + s45), 1 / alfa);
266: final double deltav = Math.asin((cs * Math.sin(d))
267: / Math.cos(u));
268: final double lambda = -deltav / alfa;
269: double phi = 0;
270: double fi1 = u;
271:
272: // iteration calculation
273: for (int i = MAXIMUM_ITERATIONS;;) {
274: fi1 = phi;
275: final double esf = excentricity * Math.sin(fi1);
276: phi = 2. * (Math.atan(kau
277: * Math.pow((1. + esf) / (1. - esf),
278: excentricity / 2.)) - s45);
279: if (Math.abs(fi1 - phi) <= ITERATION_TOLERANCE) {
280: break;
281: }
282: if (--i < 0) {
283: throw new ProjectionException(Errors
284: .format(ErrorKeys.NO_CONVERGENCE));
285: }
286: }
287:
288: if (ptDst != null) {
289: ptDst.setLocation(lambda, phi);
290: return ptDst;
291: }
292: return new Point2D.Double(lambda, phi);
293: }
294:
295: //////////////////////////////////////////////////////////////////////////////////////////
296: //////////////////////////////////////////////////////////////////////////////////////////
297: //////// ////////
298: //////// PROVIDERS ////////
299: //////// ////////
300: //////////////////////////////////////////////////////////////////////////////////////////
301: //////////////////////////////////////////////////////////////////////////////////////////
302:
303: /**
304: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
305: * provider} for an {@linkplain Krovak Krovak} projection (EPSG code 9819).
306: *
307: * @since 2.4
308: * @version $Id: Krovak.java 24563 2007-02-23 00:20:43Z desruisseaux $
309: * @author Jan Jezek
310: *
311: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
312: */
313: public static class Provider extends AbstractProvider {
314: /**
315: * The operation parameter descriptor for the {@linkplain #latitudeOfOrigin
316: * latitude of origin} parameter value. Valid values range is from -90 to 90.
317: * Default value is 49.5.
318: */
319: public static final ParameterDescriptor LATITUDE_OF_CENTER = createDescriptor(
320: new NamedIdentifier[] {
321: new NamedIdentifier(Citations.OGC,
322: "latitude_of_center"),
323: new NamedIdentifier(Citations.EPSG,
324: "Latitude of projection centre"),
325: new NamedIdentifier(Citations.GEOTIFF,
326: "CenterLat") }, 49.5, -90, 90,
327: NonSI.DEGREE_ANGLE);
328:
329: /**
330: * The operation parameter descriptor for the {@linkplain #centralMeridian central
331: * meridian} parameter value. Valid values range is from -180 to 180. Default value
332: * is 24�50' (= 42�50' from Ferro prime meridian).
333: */
334: public static final ParameterDescriptor LONGITUDE_OF_CENTER = createDescriptor(
335: new NamedIdentifier[] {
336: new NamedIdentifier(Citations.OGC,
337: "longitude_of_center"),
338: new NamedIdentifier(Citations.EPSG,
339: "Longitude of projection centre"),
340: new NamedIdentifier(Citations.GEOTIFF,
341: "CenterLong") },
342: 42.5 - 17.66666666666667, -180, 180, NonSI.DEGREE_ANGLE);
343:
344: /**
345: * The operation parameter descriptor for the {@linkplain #azimuth azimuth} parameter
346: * value. Valid values range is from -90 to 90. Default value is 30.28813972222.
347: */
348: public static final ParameterDescriptor AZIMUTH = createDescriptor(
349: new NamedIdentifier[] {
350: new NamedIdentifier(Citations.OGC, "azimuth"),
351: new NamedIdentifier(Citations.EPSG,
352: "Azimuth of initial line"),
353: new NamedIdentifier(Citations.GEOTIFF,
354: "AzimuthAngle") }, 30.28813972222222,
355: 0, 360, NonSI.DEGREE_ANGLE);
356:
357: /**
358: * The operation parameter descriptor for the pseudo {@linkplain #pseudoStandardParallel
359: * pseudi standard parallel} parameter value. Valid values range is from -90 to 90.
360: * Default value is 78.5.
361: */
362: public static final ParameterDescriptor PSEUDO_STANDARD_PARALLEL = createDescriptor(
363: new NamedIdentifier[] {
364: new NamedIdentifier(Citations.OGC,
365: "pseudo_standard_parallel_1"),
366: new NamedIdentifier(Citations.EPSG,
367: "Latitude of Pseudo Standard Parallel") },
368: 78.5, -90, 90, NonSI.DEGREE_ANGLE);
369:
370: /**
371: * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
372: * parameter value. Valid values range is from 0 to infinity. Default value is 1.
373: */
374: public static final ParameterDescriptor SCALE_FACTOR = createDescriptor(
375: new NamedIdentifier[] {
376: new NamedIdentifier(Citations.OGC,
377: "scale_factor"),
378: new NamedIdentifier(Citations.EPSG,
379: "Scale factor on pseudo standard parallel"),
380: new NamedIdentifier(Citations.GEOTIFF,
381: "ScaleAtCenter") }, 0.9999, 0,
382: Double.POSITIVE_INFINITY, Unit.ONE);
383:
384: /**
385: * The parameters group.
386: */
387: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
388: new NamedIdentifier[] {
389: new NamedIdentifier(Citations.OGC, "Krovak"),
390: new NamedIdentifier(Citations.GEOTIFF, "Krovak"),
391: new NamedIdentifier(Citations.EPSG,
392: "Krovak Oblique Conformal Conic"),
393: new NamedIdentifier(Citations.EPSG,
394: "Krovak Oblique Conic Conformal"),
395: new NamedIdentifier(Citations.EPSG, "9819"), },
396: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
397: LATITUDE_OF_CENTER, LONGITUDE_OF_CENTER,
398: AZIMUTH, PSEUDO_STANDARD_PARALLEL,
399: SCALE_FACTOR, FALSE_EASTING, FALSE_NORTHING });
400:
401: /**
402: * Constructs a new provider.
403: */
404: public Provider() {
405: super (PARAMETERS);
406: }
407:
408: /**
409: * Returns the operation type for this map projection.
410: */
411: public Class getOperationType() {
412: return ConicProjection.class;
413: }
414:
415: /**
416: * Creates a transform from the specified group of parameter values.
417: *
418: * @param parameters The group of parameter values.
419: * @return The created math transform.
420: * @throws ParameterNotFoundException if a required parameter was not found.
421: */
422: public MathTransform createMathTransform(
423: final ParameterValueGroup parameters)
424: throws ParameterNotFoundException {
425: return new Krovak(parameters);
426: }
427: }
428:
429: /**
430: * Returns a hash value for this projection.
431: */
432: public int hashCode() {
433: final long code = Double.doubleToLongBits(azimuth)
434: ^ Double.doubleToLongBits(pseudoStandardParallel);
435: return ((int) code ^ (int) (code >>> 32)) + 37
436: * super .hashCode();
437: }
438:
439: /**
440: * Compares the specified object with this map projection for equality.
441: */
442: public boolean equals(final Object object) {
443: if (object == this ) {
444: // Slight optimization
445: return true;
446: }
447: if (super .equals(object)) {
448: final Krovak that = (Krovak) object;
449: return equals(azimuth, that.azimuth)
450: && equals(pseudoStandardParallel,
451: that.pseudoStandardParallel);
452: }
453: return false;
454: }
455: }
|