/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 1999-2008, Open Source Geospatial Foundation (OSGeo)
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation;
* version 2.1 of the License.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* This package contains formulas from the PROJ package of USGS.
* USGS's work is fully acknowledged here. This derived work has
* been relicensed under LGPL with Frank Warmerdam's permission.
*/
package org.geotools.referencing.operation.projection;
import java.awt.geom.Point2D;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.parameter.ParameterNotFoundException;
import org.geotools.resources.i18n.ErrorKeys;
import static java.lang.Math.*;
/**
* The USGS oblique/equatorial case of the Stereographic projection. This is similar but
* <strong>NOT</strong> equal to EPSG code 9809 ({@code "Oblique_Stereographic"} EPSG name).
* The later is rather implemented by {@link ObliqueStereographic}.
* <p>
* This class is not public in order to keep names that closely match the ones in common usage
* (i.e. this projection is called just "Stereographic" in ESRI). Furthermore, the "USGS" name
* is not really accurate for a class to be extended by {@link ObliqueStereographic}.
*
* @since 2.4
* @source $URL$
* @version $Id$
* @author Gerald I. Evenden (for original code in Proj4)
* @author André Gosselin
* @author Martin Desruisseaux (PMO, IRD)
* @author Rueben Schulz
*/
class StereographicUSGS extends Stereographic {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = 948619442800459871L;
/**
* Maximum number of iterations for iterative computations.
*/
private static final int MAXIMUM_ITERATIONS = 15;
/**
* Difference allowed in iterative computations.
*/
private static final double ITERATION_TOLERANCE = 1E-10;
/**
* Maximum difference allowed when comparing real numbers.
*/
private static final double EPSILON = 1E-6;
/**
* Constants used for the oblique projections. All those constants are completly determined by
* {@link #latitudeOfOrigin}. Concequently, there is no need to test them in {@link #hashCode}
* or {@link #equals} methods.
*/
final double k0, sinphi0, cosphi0, chi1, sinChi1, cosChi1;
/**
* Constructs an oblique stereographic projection (USGS equations).
*
* @param parameters The group of parameter values.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
protected StereographicUSGS(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
this(parameters, Provider.PARAMETERS);
}
/**
* Constructs an oblique stereographic projection (USGS equations).
*
* @param parameters The group of parameter values.
* @param descriptor The expected parameter descriptor.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
StereographicUSGS(final ParameterValueGroup parameters,
final ParameterDescriptorGroup descriptor)
throws ParameterNotFoundException
{
super(parameters, descriptor);
if (abs(latitudeOfOrigin) < EPSILON) { // Equatorial
latitudeOfOrigin = 0;
cosphi0 = 1.0;
sinphi0 = 0.0;
chi1 = 0.0;
cosChi1 = 1.0;
sinChi1 = 0.0;
} else { // Oblique
cosphi0 = cos(latitudeOfOrigin);
sinphi0 = sin(latitudeOfOrigin);
chi1 = 2.0 * atan(ssfn(latitudeOfOrigin, sinphi0)) - PI/2;
cosChi1 = cos(chi1);
sinChi1 = sin(chi1);
}
// part of (14 - 15)
k0 = 2.0 * msfn(sinphi0, cosphi0);
}
/**
* Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
* (units in radians) and stores the result in {@code ptDst} (linear distance
* on a unit sphere).
*/
protected Point2D transformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
final double chi = 2.0 * atan(ssfn(y, sin(y))) - PI/2;
final double sinChi = sin(chi);
final double cosChi = cos(chi);
final double cosChi_cosLon = cosChi * cos(x);
final double A = k0 / cosChi1 / (1 + sinChi1*sinChi + cosChi1*cosChi_cosLon);
x = A * cosChi * sin(x);
y = A * (cosChi1 * sinChi - sinChi1 * cosChi_cosLon);
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Transforms the specified (<var>x</var>,<var>y</var>) coordinates
* and stores the result in {@code ptDst}.
*/
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
final double rho = hypot(x, y);
final double ce = 2.0 * atan2(rho * cosChi1, k0);
final double cosce = cos(ce);
final double since = sin(ce);
final boolean rhoIs0 = abs(rho) < EPSILON;
final double chi = rhoIs0 ? chi1 : asin(cosce*sinChi1 + (y*since*cosChi1 / rho));
final double tp = tan(PI/4.0 + chi/2.0);
// parts of (21-36) used to calculate longitude
final double t = x*since;
final double ct = rho*cosChi1*cosce - y*sinChi1*since;
// Compute latitude using iterative technique (3-4)
final double halfe = excentricity / 2.0;
double phi0 = chi;
for (int i=MAXIMUM_ITERATIONS;;) {
final double esinphi = excentricity * sin(phi0);
final double phi = 2*atan(tp*pow((1+esinphi)/(1-esinphi), halfe)) - PI/2;
if (abs(phi - phi0) < ITERATION_TOLERANCE) {
// TODO: checking rho may be redundant
x = rhoIs0 || (abs(t)<EPSILON && abs(ct)<EPSILON) ? 0.0 : atan2(t, ct);
y = phi;
break;
}
phi0 = phi;
if (--i < 0) {
throw new ProjectionException(ErrorKeys.NO_CONVERGENCE);
}
}
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Maximal error (in metres) tolerated for assertions, if enabled.
*/
@Override
protected double getToleranceForAssertions(final double longitude, final double latitude) {
final double delta = abs(longitude - centralMeridian)/2 +
abs(latitude - latitudeOfOrigin);
if (delta > 40) {
return 0.5;
}
if (delta > 15) {
return 0.1;
}
return super.getToleranceForAssertions(longitude, latitude);
}
/**
* Computes part of function (3-1) from Snyder.
*/
final double ssfn(double phi, double sinphi) {
sinphi *= excentricity;
return tan(PI/4 + phi/2.0) * pow((1-sinphi) / (1+sinphi), excentricity/2.0);
}
/**
* Provides the transform equations for the spherical case of the
* Stereographic projection.
*
* @version $Id$
* @author Martin Desruisseaux (PMO, IRD)
* @author Rueben Schulz
*/
static final class Spherical extends StereographicUSGS {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = -8558594307755820783L;
/**
* A constant used in the transformations. This constant hides the {@code k0}
* constant from the ellipsoidal case. The spherical and ellipsoidal {@code k0}
* are not computed in the same way, and we preserve the ellipsoidal {@code k0}
* in {@link Stereographic} in order to allow assertions to work.
*/
private static final double k0 = 2;
/**
* Constructs a spherical oblique stereographic projection.
*
* @param parameters The group of parameter values.
* @param descriptor The expected parameter descriptor.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
Spherical(final ParameterValueGroup parameters, final ParameterDescriptorGroup descriptor)
throws ParameterNotFoundException
{
super(parameters, descriptor);
ensureSpherical();
}
/**
* Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
* (units in radians) and stores the result in {@code ptDst} (linear distance
* on a unit sphere).
*/
@Override
protected Point2D transformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
// Compute using ellipsoidal formulas, for comparaison later.
assert (ptDst = super.transformNormalized(x, y, ptDst)) != null;
final double coslat = cos(y);
final double sinlat = sin(y);
final double coslon = cos(x);
double f = 1.0 + sinphi0*sinlat + cosphi0*coslat*coslon; // (21-4)
if (f < EPSILON) {
throw new ProjectionException(ErrorKeys.VALUE_TEND_TOWARD_INFINITY);
}
f = k0 / f;
x = f * coslat * sin(x); // (21-2)
y = f * (cosphi0 * sinlat - sinphi0 * coslat * coslon); // (21-3)
assert checkTransform(x, y, ptDst);
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Transforms the specified (<var>x</var>,<var>y</var>) coordinates
* and stores the result in {@code ptDst}.
*/
@Override
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
// Compute using ellipsoidal formulas, for comparaison later.
assert (ptDst = super.inverseTransformNormalized(x, y, ptDst)) != null;
final double rho = hypot(x, y);
if (abs(rho) < EPSILON) {
y = latitudeOfOrigin;
x = 0.0;
} else {
final double c = 2.0 * atan(rho/k0);
final double cosc = cos(c);
final double sinc = sin(c);
final double ct = rho*cosphi0*cosc - y*sinphi0*sinc; // (20-15)
final double t = x*sinc; // (20-15)
y = asin(cosc*sinphi0 + y*sinc*cosphi0/rho); // (20-14)
x = (abs(ct)<EPSILON && abs(t)<EPSILON) ? 0.0 : atan2(t, ct);
}
assert checkInverseTransform(x, y, ptDst);
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
}
}