Statistics
| Revision:

root / branches / libProjection_v2_0_prep / libraries / libJCRS / src / org / geotools / referencing / operation / projection / IdrKrovak.java @ 27137

History | View | Annotate | Download (19.6 KB)

1
/*
2
 *    Geotools2 - OpenSource mapping toolkit
3
 *    http://geotools.org
4
 *    (C) 2002-2005, Geotools Project Managment Committee (PMC)
5
 *
6
 *    This library is free software; you can redistribute it and/or
7
 *    modify it under the terms of the GNU Lesser General Public
8
 *    License as published by the Free Software Foundation;
9
 *    version 2.1 of the License.
10
 *
11
 *    This library is distributed in the hope that it will be useful,
12
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
13
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14
 *    Lesser General Public License for more details.
15
 */
16
package org.geotools.referencing.operation.projection;
17

    
18
// J2SE dependencies and extensions
19
import java.awt.geom.Point2D;
20
import java.util.Collection;
21

    
22
import javax.units.SI;
23
import javax.units.NonSI;
24
import javax.units.Unit;
25

    
26
import org.geotools.metadata.iso.citation.CitationImpl;
27
import org.geotools.referencing.NamedIdentifier;
28
import org.geotools.resources.XMath;
29
import org.geotools.resources.cts.ResourceKeys;
30
import org.geotools.resources.cts.Resources;
31
import org.opengis.parameter.ParameterDescriptor;
32
import org.opengis.parameter.ParameterDescriptorGroup;
33
import org.opengis.parameter.ParameterNotFoundException;
34
import org.opengis.parameter.ParameterValueGroup;
35
import org.opengis.referencing.operation.ConicProjection;
36
import org.opengis.referencing.operation.MathTransform;
37

    
38

    
39
/**
40
 * Krovak Oblique Conformal Conic projection (EPSG code 9819). This projection is
41
 * used in the Czech Republic and Slovakia under the name 'Krovak' projection. The
42
 * geographic coordinates on the ellipsoid are first reduced to conformal
43
 * coordinates on the conformal (Gaussian) sphere. These spherical coordinates
44
 * are then projected onto the oblique cone and converted to grid coordinates.
45
 * The pseudo standard parallel is defined on the conformal sphere after its
46
 * rotation, to obtain the oblique aspect of the projection. It is then the
47
 * parallel on this sphere at which the map projection is true to scale; on
48
 * the ellipsoid it maps as a complex curve.
49
 *
50
 * <p>The compulsory parameters are just the ellipsoid characteristics.
51
 * All other parameters are optional and have defaults to match the common
52
 * usage with Krovak projection.</p>
53
 *
54
 * <p>In general the axis of Krovak projection are defined as westing and
55
 * southing (not easting and northing) and they are also reverted, so if the
56
 * value of projected coordinates should (and in <var>y</var>, <var>x</var>
57
 * order in Krovak) be positive the 'Axis' parameter for projection should be
58
 * defined explicitly like this (in wkt):</p>
59
 * 
60
 * <pre>PROJCS["S-JTSK (Ferro) / Krovak",  
61
 *         .                                                              
62
 *         .                                                              
63
 *         .
64
 *                                                                       
65
 *     PROJECTION["Krovak"]                                         
66
 *     PARAMETER["semi_major", 6377397.155],  
67
 *     PARAMETER["semi_minor", 6356078.963],                   
68
 *     UNIT["meter",1.0],                                  
69
 *     AXIS["x", WEST],                     
70
 *     AXIS["y", SOUTH]]                                              
71
 *     </pre>Axis in Krovak:
72
 * <pre>                                                              
73
 *   y<------------------+                                                                                                  
74
 *                       |                                             
75
 *    Czech. Rep.        | 
76
 *                       |                                                                   
77
 *                       x                              
78
 * </pre>
79
 * <p>By default, the axis are 'easting, northing' so the values of projected coordinates
80
 * are negative and in (and in <var>y</var>, <var>x</var> order in Krovak - it is cold
81
 * Krovak GIS version).</p>
82
 *
83
 * <p><strong>References:</strong>
84
 *  <ul>
85
 *      <li>Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
86
 *      Relevant files is: {@code PJ_krovak.c}</li>
87
 *      <li>"Coordinate Conversions and Transformations including Formulas" available at, <A
88
 *      HREF="http://www.remotesensing.org/geotiff/proj_list/guid7.html">http://www.remotesensing.org/geotiff/proj_list/guid7.html</A></li>
89
 *  </ul>
90
 * </p>
91
 *
92
 * @since 2.4
93
 * @version $Id: IdrKrovak.java 12357 2007-06-27 09:05:41Z dguerrero $
94
 * @source $URL: http://svn.geotools.org/geotools/tags/2.3.0/module/referencing/src/org/geotools/referencing/operation/projection/Krovak.java $
95
 * @author Jan Jezek
96
 * @author Martin Desruisseaux
97
 *
98
 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/krovak.html">Krovak on
99
 *      RemoteSensing.org </A>
100
 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/guid7.html">Krovak on "Coordinate
101
 *      Conversions and Transformations including Formulas"</A>
102
 */
103
public class IdrKrovak extends MapProjection {
104
    /**
105
     * Maximum number of iterations for iterative computations.
106
     */
107
    private static final int MAXIMUM_ITERATIONS = 15;
108
    
109
    /**
110
     * When to stop the iteration.
111
     */
112
    private static final double ITERATION_TOLERANCE = 1E-11;
113

    
114
    /**
115
     * Azimuth of the centre line passing through the centre of the projection.
116
     * This is equals to the co-latitude of the cone axis at point of intersection
117
     * with the ellipsoid.
118
     */
119
    protected final double azimuth;
120

    
121
    /**
122
     * Latitude of pseudo standard parallel.
123
     */
124
    protected final double pseudoStandardParallel;
125

    
126
    /**
127
     * Useful variables calculated from parameters defined by user.
128
     */
129
    private final double sinAzim, cosAzim, n, tanS2, alfa, hae, k1, ka, 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 IdrKrovak(final ParameterValueGroup parameters) throws ParameterNotFoundException {
143
        super(parameters);
144
        final Collection expected = getParameterDescriptors().descriptors();
145
        // Fetch parameters from user input.
146
        latitudeOfOrigin       = doubleValue(expected, Provider.LATITUDE_OF_CENTER,       parameters);
147
        centralMeridian        = doubleValue(expected, Provider.LONGITUDE_OF_CENTER,      parameters);
148
        azimuth                = doubleValue(expected, Provider.AZIMUTH,                  parameters);
149
        pseudoStandardParallel = doubleValue(expected, Provider.PSEUDO_STANDARD_PARALLEL, parameters);
150
        scaleFactor            = doubleValue(expected, Provider.SCALE_FACTOR,             parameters);
151
        ensureLatitudeInRange (Provider.LATITUDE_OF_CENTER,  latitudeOfOrigin, false);
152
        ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTER, centralMeridian,  false);
153

    
154
        // Calculates useful constants.
155
        sinAzim = Math.sin(azimuth);
156
        cosAzim = Math.cos(azimuth);
157
        n       = Math.sin(pseudoStandardParallel);
158
        tanS2   = Math.tan(pseudoStandardParallel / 2 + s45);
159

    
160
        final double sinLat, cosLat, cosL2, u0;
161
        sinLat = Math.sin(latitudeOfOrigin);
162
        cosLat = Math.cos(latitudeOfOrigin);
163
        cosL2  = cosLat * cosLat;
164
        alfa   = Math.sqrt(1 + ((excentricitySquared * (cosL2*cosL2)) / (1 - excentricitySquared)));
165
        hae    = alfa * excentricity / 2;
166
        u0     = Math.asin(sinLat / alfa);
167

    
168
        final double g, esl;
169
        esl = excentricity * sinLat;
170
        g   = Math.pow((1 - esl) / (1 + esl), (alfa * excentricity) / 2);
171
        k1  = Math.pow(Math.tan(latitudeOfOrigin/2 + s45), alfa) * g / Math.tan(u0/2 + s45);
172
        ka  = Math.pow(1 / k1, -1 / alfa);
173

    
174
        final double radius;
175
        radius = Math.sqrt(1 - excentricitySquared) / (1 - (excentricitySquared * (sinLat * sinLat)));
176

    
177
        ro0 = scaleFactor * radius / Math.tan(pseudoStandardParallel);
178
        rop = ro0 * Math.pow(tanS2, n);
179
    }
180

    
181
    /**
182
     * {@inheritDoc}
183
     */
184
    public ParameterDescriptorGroup getParameterDescriptors() {
185
        return Provider.PARAMETERS;
186
    }
187

    
188
    /**
189
     * {@inheritDoc}
190
     */
191
    public ParameterValueGroup getParameterValues() {
192
        final Collection expected = getParameterDescriptors().descriptors();
193
        final ParameterValueGroup values = super.getParameterValues();
194
        set(expected, Provider.LATITUDE_OF_CENTER,       values, latitudeOfOrigin      );
195
        set(expected, Provider.LONGITUDE_OF_CENTER,      values, centralMeridian       );
196
        set(expected, Provider.AZIMUTH,                  values, azimuth               );
197
        set(expected, Provider.PSEUDO_STANDARD_PARALLEL, values, pseudoStandardParallel);
198
        set(expected, Provider.SCALE_FACTOR,             values, scaleFactor           );
199
        return values;
200
    }
201

    
202
    /**
203
     * Transforms the specified (<var>&lambda;</var>,<var>&phi;</var>) coordinates
204
     * (units in radians) and stores the result in {@code ptDst} (linear distance
205
     * on a unit sphere).
206
     */
207
    protected Point2D transformNormalized(final double lambda, final double phi, Point2D ptDst)
208
            throws ProjectionException
209
    {
210
        final double esp = excentricity * Math.sin(phi);
211
        final double gfi = Math.pow(((1. - esp) / (1. + esp)), hae);
212
        final double u   = 2 * (Math.atan(Math.pow(Math.tan(phi/2 + s45), alfa) / k1 * gfi) - s45);
213
        final double deltav = -lambda * alfa;
214
        final double cosU = Math.cos(u);
215
        final double s = Math.asin((cosAzim * Math.sin(u)) + (sinAzim * cosU * Math.cos(deltav)));
216
        final double d = Math.asin(cosU * Math.sin(deltav) / Math.cos(s));
217
        final double eps = n * d;
218
        final double ro = rop / Math.pow(Math.tan(s/2 + s45), n);
219

    
220
        /* x and y are reverted  */
221
        final double y = -(ro * Math.cos(eps));
222
        final double x = -(ro * Math.sin(eps));
223

    
224
        if (ptDst != null) {
225
            ptDst.setLocation(x, y);
226
            return ptDst;
227
        }
228
        return new Point2D.Double(x, y);
229
    }
230

    
231
    /**
232
     * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
233
     * and stores the result in {@code ptDst}.
234
     */
235
    protected Point2D inverseTransformNormalized(final double x, final double y, Point2D ptDst)
236
            throws ProjectionException
237
    {
238
        // x -> southing, y -> westing
239
        final double ro  = XMath.hypot(x, y);
240
        final double eps = Math.atan2(-x, -y);
241
        final double d   = eps / n;
242
        final double s   = 2 * (Math.atan(Math.pow(ro0/ro, 1/n) * tanS2) - s45);
243
        final double cs  = Math.cos(s);
244
        final double u   = Math.asin((cosAzim * Math.sin(s)) - (sinAzim * cs * Math.cos(d)));
245
        final double kau = ka * Math.pow(Math.tan((u / 2.) + s45), 1 / alfa);
246
        final double deltav = Math.asin((cs * Math.sin(d)) / Math.cos(u));
247
        final double lambda = -deltav / alfa;
248
        double phi = 0;
249
        double fi1 = u;
250

    
251
        // iteration calculation
252
        for (int i=MAXIMUM_ITERATIONS;;) {
253
            fi1 = phi;
254
            final double esf = excentricity * Math.sin(fi1);
255
            phi = 2. * (Math.atan(kau * Math.pow((1. + esf) / (1. - esf), excentricity/2.)) - s45);
256
            if (Math.abs(fi1 - phi) <= ITERATION_TOLERANCE) {
257
                break;
258
            }
259
            if (--i < 0) {
260
                    throw new ProjectionException(Resources.format(ResourceKeys.ERROR_NO_CONVERGENCE));
261
            }
262
        }
263

    
264
        if (ptDst != null) {
265
            ptDst.setLocation(lambda, phi);
266
            return ptDst;
267
        }
268
        return new Point2D.Double(lambda, phi);
269
    }
270

    
271

    
272

    
273

    
274
    //////////////////////////////////////////////////////////////////////////////////////////
275
    //////////////////////////////////////////////////////////////////////////////////////////
276
    ////////                                                                          ////////
277
    ////////                                 PROVIDER                                 ////////
278
    ////////                                                                          ////////
279
    //////////////////////////////////////////////////////////////////////////////////////////
280
    //////////////////////////////////////////////////////////////////////////////////////////
281

    
282
    /**
283
     * The {@link org.geotools.referencing.operation.MathTransformProvider}
284
     * for an {@link IdrKrovak krovak} projection.
285
     *
286
     * @author jezekjan
287
     *
288
     * @see org.geotools.referencing.operation.DefaultMathTransformFactory
289
     */
290
    public static class Provider extends AbstractProvider {
291
        /**
292
         * The operation parameter descriptor for the {@linkPlain #latitudeOfOrigin
293
         * latitude of origin} parameter value. Valid values range is from -90 to 90.
294
         * Default value is 49.5.
295
         */
296
        public static final ParameterDescriptor LATITUDE_OF_CENTER = createDescriptor(
297
                new NamedIdentifier[] {
298
                    new NamedIdentifier(CitationImpl.OGC,     "latitude_of_center"),
299
                    new NamedIdentifier(CitationImpl.EPSG,    "Latitude of projection center"),
300
                    new NamedIdentifier(CitationImpl.EPSG,    "Latitude of projection centre"),
301
                    new NamedIdentifier(CitationImpl.GEOTIFF, "CenterLat")
302
                }, 49.5, -90, 90, NonSI.DEGREE_ANGLE);
303

    
304
        /**
305
         * The operation parameter descriptor for the {@linkPlain #centralMeridian central
306
         * meridian} parameter value. Valid values range is from -180 to 180. Default value
307
         * is 2450' (= 4250' from Ferro prime meridian).
308
         */
309
        public static final ParameterDescriptor LONGITUDE_OF_CENTER = createDescriptor(
310
                new NamedIdentifier[] {
311
                    new NamedIdentifier(CitationImpl.OGC,     "longitude_of_center"),
312
                    new NamedIdentifier(CitationImpl.EPSG,    "Longitude of projection center"),
313
                    new NamedIdentifier(CitationImpl.EPSG,    "Longitude of projection centre"),
314
                    new NamedIdentifier(CitationImpl.GEOTIFF, "CenterLong")
315
                }, 42.5-17.66666666666667, -180, 180, NonSI.DEGREE_ANGLE);
316

    
317
        /**
318
         * The operation parameter descriptor for the {@linkPlain #azimuth azimuth} parameter
319
         * value. Valid values range is from -90 to 90. Default value is 30.28813972222.
320
         */
321
        public static final ParameterDescriptor AZIMUTH = createDescriptor(
322
                new NamedIdentifier[] {
323
                    new NamedIdentifier(CitationImpl.OGC,     "azimuth"),
324
                    new NamedIdentifier(CitationImpl.EPSG,    "Azimuth of the center line"),
325
                    new NamedIdentifier(CitationImpl.EPSG,    "Azimuth of initial line"),
326
                    new NamedIdentifier(CitationImpl.GEOTIFF, "AzimuthAngle")
327
                }, 30.28813972222222, 0, 360, NonSI.DEGREE_ANGLE);
328

    
329
        /**
330
         * The operation parameter descriptor for the pseudo {@linkplain #pseudoStandardParallel
331
         * pseudi standard parallel} parameter value. Valid values range is from -90 to 90.
332
         * Default value is 78.5.
333
         */
334
        public static final ParameterDescriptor PSEUDO_STANDARD_PARALLEL =
335
                createDescriptor(new NamedIdentifier[] {
336
                    new NamedIdentifier(CitationImpl.OGC,  "pseudo_standard_parallel_1"),
337
                    new NamedIdentifier(CitationImpl.EPSG, "Latitude of Pseudo Standard Parallel")
338
                }, 78.5, -90, 90, NonSI.DEGREE_ANGLE);
339

    
340
        /**
341
         * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
342
         * parameter value. Valid values range is from 0 to infinity. Default value is 1.
343
         */
344
        public static final ParameterDescriptor SCALE_FACTOR = createDescriptor(
345
                new NamedIdentifier[] {
346
                    new NamedIdentifier(CitationImpl.OGC,  "scale_factor"),
347
                    new NamedIdentifier(CitationImpl.EPSG, "Scale factor on the pseudo standard line"),
348
                    new NamedIdentifier(CitationImpl.EPSG, "Scale factor on pseudo standard parallel"),
349
                    new NamedIdentifier(CitationImpl.GEOTIFF, "ScaleAtCenter")
350
                }, 0.9999, 0, Double.POSITIVE_INFINITY, Unit.ONE);
351

    
352
        /**
353
         * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
354
         * parameter value. Valid values range is from 0 to infinity. Default value is 1.
355
         */
356
        public static final ParameterDescriptor FALSEEASTING = createDescriptor(
357
                new NamedIdentifier[] {
358
                    new NamedIdentifier(CitationImpl.EPSG, "Easting at projection centre")
359
                }, 0.0, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, SI.METER);
360

    
361
        /**
362
         * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
363
         * parameter value. Valid values range is from 0 to infinity. Default value is 1.
364
         */
365
        public static final ParameterDescriptor FALSENORTHING = createDescriptor(
366
                new NamedIdentifier[] {
367
                    new NamedIdentifier(CitationImpl.EPSG, "Northing at projection centre")
368
                }, 0.0, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, SI.METER);
369

    
370
        /**
371
         * The parameters group.
372
         */
373
        static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
374
            new NamedIdentifier(CitationImpl.OGC,     "Krovak"),
375
            new NamedIdentifier(CitationImpl.GEOTIFF, "Krovak"),
376
            new NamedIdentifier(CitationImpl.EPSG,    "Krovak Oblique Conformal Conic"),
377
            new NamedIdentifier(CitationImpl.EPSG,    "Krovak Oblique Conic Conformal"),
378
            new NamedIdentifier(CitationImpl.EPSG,    "9819"),
379
            new NamedIdentifier(new CitationImpl("IDR"), "IDR")
380
                },
381
                new ParameterDescriptor[] {
382
                    SEMI_MAJOR, SEMI_MINOR, LATITUDE_OF_CENTER, LONGITUDE_OF_CENTER,
383
                    AZIMUTH, PSEUDO_STANDARD_PARALLEL, SCALE_FACTOR,
384
                    FALSEEASTING, FALSENORTHING
385
                });
386

    
387
        /**
388
         * Constructs a new provider. 
389
         */
390
        public Provider() {
391
            super(PARAMETERS);
392
        }
393

    
394
        /**
395
         * Returns the operation type for this map projection.
396
         */
397
        protected Class getOperationType() {
398
            return ConicProjection.class;
399
        }
400

    
401
        /**
402
         * Creates a transform from the specified group of parameter values.
403
         *
404
         * @param parameters The group of parameter values.
405
         * @return The created math transform.
406
         * @throws ParameterNotFoundException if a required parameter was not found.
407
         */
408
        public MathTransform createMathTransform(final ParameterValueGroup parameters)
409
                throws ParameterNotFoundException
410
        {
411
            return new IdrKrovak(parameters);
412
        }
413
    }
414

    
415
    /**
416
     * Returns a hash value for this projection.
417
     */
418
    public int hashCode() {
419
        final long code = Double.doubleToLongBits(azimuth) ^
420
                          Double.doubleToLongBits(pseudoStandardParallel);
421
        return ((int)code ^ (int)(code >>> 32)) + 37*super.hashCode();
422
    }
423
    
424
    /**
425
     * Compares the specified object with this map projection for equality.
426
     */
427
    public boolean equals(final Object object) {
428
        if (object == this) {
429
            // Slight optimization
430
            return true;
431
        }
432
        if (super.equals(object)) {
433
            final IdrKrovak that = (IdrKrovak) object;
434
            return equals(azimuth, that.azimuth) &&
435
                   equals(pseudoStandardParallel, that.pseudoStandardParallel);
436
        }
437
        return false;
438
    }
439
}