001 //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_testing/src/org/deegree/crs/projections/azimuthal/StereographicAlternative.java $
002 /*----------------------------------------------------------------------------
003 This file is part of deegree, http://deegree.org/
004 Copyright (C) 2001-2009 by:
005 Department of Geography, University of Bonn
006 and
007 lat/lon GmbH
008
009 This library is free software; you can redistribute it and/or modify it under
010 the terms of the GNU Lesser General Public License as published by the Free
011 Software Foundation; either version 2.1 of the License, or (at your option)
012 any later version.
013 This library is distributed in the hope that it will be useful, but WITHOUT
014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
015 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
016 details.
017 You should have received a copy of the GNU Lesser General Public License
018 along with this library; if not, write to the Free Software Foundation, Inc.,
019 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
020
021 Contact information:
022
023 lat/lon GmbH
024 Aennchenstr. 19, 53177 Bonn
025 Germany
026 http://lat-lon.de/
027
028 Department of Geography, University of Bonn
029 Prof. Dr. Klaus Greve
030 Postfach 1147, 53001 Bonn
031 Germany
032 http://www.geographie.uni-bonn.de/deegree/
033
034 e-mail: info@deegree.org
035 ----------------------------------------------------------------------------*/
036
037 package org.deegree.crs.projections.azimuthal;
038
039 import static org.deegree.crs.projections.ProjectionUtils.EPS11;
040 import static org.deegree.crs.projections.ProjectionUtils.HALFPI;
041 import static org.deegree.crs.projections.ProjectionUtils.QUARTERPI;
042
043 import javax.vecmath.Point2d;
044
045 import org.deegree.crs.Identifiable;
046 import org.deegree.crs.components.Unit;
047 import org.deegree.crs.coordinatesystems.GeographicCRS;
048 import org.deegree.crs.exceptions.ProjectionException;
049 import org.deegree.framework.log.ILogger;
050 import org.deegree.framework.log.LoggerFactory;
051
052 /**
053 * <code>StereographicAlternative</code> projection may be imagined to be a projection of the earth's surface onto a
054 * plane in contact with the earth at a single tangent point from the opposite end of the diameter through that tangent
055 * point.
056 * <p>
057 * An alternative approach is given by Snyder {@link StereographicAzimuthal}, where, instead of defining a single
058 * conformal sphere at the origin point, the conformal latitude at each point on the ellipsoid is computed. The
059 * conformal longitude is then always equivalent to the geodetic longitude. This approach is a valid alternative to the
060 * above, but gives slightly different results away from the origin point. It is therefore considered by EPSG to be a
061 * different projection method. Hence this implementation.
062 * </p>
063 * <p>
064 * This projection is best known in its polar form and is frequently used for mapping polar areas where it complements
065 * the Universal Transverse Mercator used for lower latitudes. Its spherical form has also been widely used by the US
066 * Geological Survey for planetary mapping and the mapping at small scale of continental hydrocarbon provinces. In its
067 * transverse or oblique ellipsoidal forms it is useful for mapping limited areas centered on the point where the plane
068 * of the projection is regarded as tangential to the ellipsoid., e.g. the Netherlands. The tangent point is the origin
069 * of the projected coordinate system and the meridian through it is regarded as the central meridian. In order to
070 * reduce the scale error at the extremities of the projection area it is usual to introduce a scale factor of less than
071 * unity at the origin such that a unit scale factor applies on a near circle centered at the origin and some distance
072 * from it.
073 * </p>
074 *
075 * <p>
076 * The coordinate transformation from geographical to projected coordinates is executed via the distance and azimuth of
077 * the point from the center point or origin. For a sphere the formulas are relatively simple. For the ellipsoid the
078 * same formulas are used but with auxiliary latitudes, known as conformal latitudes, substituted for the geodetic
079 * latitudes of the spherical formulas for the origin and the point .
080 * </p>
081 * <quote>from <a
082 * href="http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34j.html">http://www.posc.org/</a></quote>
083 *
084 * <p>
085 * Determinations of oblique projections on an ellipsoid can be difficult to solve and result in long, complex
086 * computations. Because conformal transformations can be made multiple time without loss of the conformal property a
087 * method of determining oblique projections involves conformal transformation of the elliptical coordinates to
088 * coordinates on a conformal sphere. The transformed coordinates can now be translated/rotated on the sphere and then
089 * converted to planar coordinates with a conformal spherical projection. is Ce/2 1 − e sin φ 2 arctan K tanC (π/4 +
090 * φ/2) − π/2 (3.6) χ = 1 + e sin φ λc = Cλ (3.7) √ 1 − e2 Rc = (3.8) 1 − e2 sin2 φ0
091 * </p>
092 * From the <a href="http://members.verizon.net/~gerald.evenden/proj4/manual.pdf">libproj4-manual</a> by Gerald I.
093 * Evenden
094 *
095 *
096 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
097 *
098 * @author last edited by: $Author: rbezema $
099 *
100 * @version $Revision: 19653 $, $Date: 2009-09-15 14:56:30 +0200 (Di, 15 Sep 2009) $
101 *
102 */
103 public class StereographicAlternative extends AzimuthalProjection {
104
105 private final static ILogger LOG = LoggerFactory.getLogger( StereographicAlternative.class );
106
107 private static final long serialVersionUID = -444957611911311171L;
108
109 private double sinc0 = 0;
110
111 private double cosc0 = 0;
112
113 /**
114 * Two times the radius of the conformal sphere, which is denoted by Gerald I. Evenden as R_c
115 */
116 private final double R2;
117
118 /**
119 * THe latitude of on the conformal sphere which is denoted by Gerald I. Evenden as Chi_0
120 */
121 private final double latitudeOnCS;
122
123 /**
124 * The exponent of the calculation of conformal latitude Chi in Gerald I. Evenden (3.6)
125 */
126 private final double clExponent;
127
128 /**
129 * The value K in Gerald I. Evenden (3.11)
130 */
131 private final double K;
132
133 /**
134 * The central geographic latitude of the projection on the conformal sphere, denoted by Gerald I. Evenden (3.9)
135 * with C
136 */
137 private final double centralGeographicLatitude;
138
139 // private double chi = 0;
140
141 // private double radiusOfCS = 0;
142
143 /**
144 * @param geographicCRS
145 * @param falseNorthing
146 * @param falseEasting
147 * @param naturalOrigin
148 * @param units
149 * @param scale
150 * @param id
151 * an identifiable instance containing information about this projection
152 */
153 public StereographicAlternative( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
154 Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
155 super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, true/* conformal */, false, id );
156
157 // if (!(P->en = pj_gauss_ini(P->e, P->phi0, &(P->phic0), &R))) E_ERROR_0;
158 double es = getSquaredEccentricity();
159 double sinPhi0 = getSinphi0();
160 double cosPhiSquare = getCosphi0() * getCosphi0();// Math.cos( getProjectionLatitude() );
161
162 // R_c
163 double radiusOfCS = Math.sqrt( 1. - es ) / ( 1. - es * sinPhi0 * sinPhi0 );
164 centralGeographicLatitude = Math.sqrt( 1. + es * cosPhiSquare * cosPhiSquare / ( 1. - es ) );
165 // xsi_0
166 latitudeOnCS = Math.asin( sinPhi0 / centralGeographicLatitude );
167 clExponent = 0.5 * centralGeographicLatitude * getEccentricity();
168 K = Math.tan( .5 * latitudeOnCS + QUARTERPI )
169 / ( Math.pow( Math.tan( .5 * getProjectionLatitude() + QUARTERPI ), centralGeographicLatitude ) * srat(
170 getEccentricity()
171 * sinPhi0,
172 clExponent ) );
173 sinc0 = Math.sin( latitudeOnCS );
174 cosc0 = Math.cos( latitudeOnCS );
175 R2 = 2 * radiusOfCS;
176 }
177
178 /**
179 * Sets the id of this projection to epsg::9809 (Oblique Stereographic)
180 *
181 * @param geographicCRS
182 * @param falseNorthing
183 * @param falseEasting
184 * @param naturalOrigin
185 * @param units
186 * @param scale
187 */
188 public StereographicAlternative( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
189 Point2d naturalOrigin, Unit units, double scale ) {
190 this( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, new Identifiable( "EPSG::9809" ) );
191 }
192
193 /*
194 * (non-Javadoc)
195 *
196 * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
197 */
198 @Override
199 public Point2d doInverseProjection( double x, double y )
200 throws ProjectionException {
201 if ( LOG.isDebug() ) {
202 LOG.logDebug( "(sterea incoming inv proj) x:" + x + ", y:" + y );
203 }
204 Point2d result = new Point2d();
205 x -= getFalseEasting();
206 y -= getFalseNorthing();
207
208 x /= getScaleFactor();
209 y /= getScaleFactor();
210 double rho = Math.hypot( x, y );
211
212 if ( rho > EPS11 ) {
213 double c = 2 * Math.atan2( rho, R2 );
214 double sinc = Math.sin( c );
215 double cosc = Math.cos( c );
216 result.y = Math.asin( cosc * sinc0 + y * sinc * cosc0 / rho );
217 result.x = Math.atan2( x * sinc, rho * cosc0 * cosc - y * sinc0 * sinc );
218 } else {
219 result.y = latitudeOnCS;
220 result.x = 0;
221 }
222 result = pj_inv_gauss( result );
223 result.x += getProjectionLongitude();
224 if ( LOG.isDebug() ) {
225 LOG.logDebug( "(outgoing) lam:" + Math.toDegrees( result.x ) + ", phi:" + Math.toDegrees( result.y ) );
226 }
227 return result;
228 }
229
230 /*
231 * (non-Javadoc)
232 *
233 * @see org.deegree.crs.projections.Projection#doProjection(double, double)
234 */
235 @Override
236 public Point2d doProjection( double lambda, double phi )
237 throws ProjectionException {
238 if ( LOG.isDebug() ) {
239 LOG.logDebug( "(sterea incoming proj) lam:" + Math.toDegrees( lambda ) + ", phi:" + Math.toDegrees( phi ) );
240 }
241 Point2d result = new Point2d();
242 lambda -= getProjectionLongitude();
243 Point2d lp = pj_gauss( lambda, phi );
244 double sinc = Math.sin( lp.y );
245 double cosc = Math.cos( lp.y );
246 double cosl = Math.cos( lp.x );
247
248 double k = getScaleFactor() * ( R2 / ( 1. + sinc0 * sinc + cosc0 * cosc * cosl ) );
249 result.x = k * cosc * Math.sin( lp.x );
250 result.y = k * ( cosc0 * sinc - sinc0 * cosc * cosl );
251
252 result.x += getFalseEasting();
253 result.y += getFalseNorthing();
254 if ( LOG.isDebug() ) {
255 LOG.logDebug( "(outgoing) x:" + result.x + ", y:" + result.y );
256 }
257 return result;
258 }
259
260 /*
261 * (non-Javadoc)
262 *
263 * @see org.deegree.crs.projections.Projection#getDeegreeSpecificName()
264 */
265 @Override
266 public String getImplementationName() {
267 return "stereographicAlternative";
268 }
269
270 private static double srat( double esinp, double exp ) {
271 return Math.pow( ( 1. - esinp ) / ( 1. + esinp ), exp );
272 }
273
274 /**
275 * To determine the inverse solution, geographic coordinates from Gaussian sphere coordinates, execute with the
276 * initial value of φi−1 = χ and φi−1 iteratively replaced by φ until |φ − φi−1 | is less than an acceptable error
277 * value. (taken from the proj-lib manual.
278 */
279 private Point2d pj_inv_gauss( Point2d slp )
280 throws ProjectionException {
281 Point2d elp = new Point2d();
282
283 elp.x = slp.x / centralGeographicLatitude;
284 double num = Math.pow( Math.tan( .5 * slp.y + QUARTERPI ) / K, 1. / centralGeographicLatitude );
285 int MAX_ITER = 20;
286 int i = MAX_ITER;
287 for ( ; i > 0; --i ) {
288 elp.y = 2. * Math.atan( num * srat( getEccentricity() * Math.sin( slp.y ), -.5 * getEccentricity() ) )
289 - HALFPI;
290 if ( Math.abs( ( elp.y - slp.y ) ) < EPS11 ) {
291 break;
292 }
293 slp.y = elp.y;
294 }
295 /* convergence failed */
296 if ( i == 0 ) {
297 throw new ProjectionException( "No convertgence while calculation the inverse gaus approximation" );
298 }
299 return elp;
300 }
301
302 /**
303 * The conformal transformation of ellipsoid coordinates (φ, λ) to conformal sphere coordinates (χ, λ_c ), where λ
304 * is relative to the longitude of projection origin, R_c is radius of the conformal sphere. χ_0 is the latitude on
305 * the conformal sphere at the central geographic latitude of the projection.
306 */
307 private Point2d pj_gauss( double lambda, double phi ) {
308 Point2d slp = new Point2d();
309 slp.y = 2.
310 * Math.atan( K * Math.pow( Math.tan( .5 * phi + QUARTERPI ), centralGeographicLatitude )
311 * srat( getEccentricity() * Math.sin( phi ), clExponent ) ) - HALFPI;
312 slp.x = centralGeographicLatitude * ( lambda );
313 return slp;
314 }
315
316 }