/*
* 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.ParameterDescriptor;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.parameter.ParameterNotFoundException;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.operation.MathTransform;
import org.geotools.resources.i18n.ErrorKeys;
import org.geotools.referencing.NamedIdentifier;
import org.geotools.metadata.iso.citation.Citations;
import static java.lang.Math.*;
/**
* Provides the transform equations for the Oblique Stereographic (EPSG code 9809).
* The formulas used below are not from the EPSG, but rather those of the
* "Oblique Stereographic Alternative" in the {@code libproj4} package
* written by Gerald Evenden. His work is acknowledged here and greatly appreciated.
*
* The forward equations used in {@code libproj4} are the same as those given in the
* UNB reports for the Double Stereographic. The inverse equations are similar,
* but use different methods to iterate for the latitude.
*
* References:
*
* - {@code libproj4} is available at
* libproj4 Miscellanea
* Relevent files are: {@code PJ_sterea.c}, {@code pj_gauss.c},
* {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}
* - Gerald Evenden.
* "Supplementary PROJ.4 Notes - Oblique Stereographic Alternative"
* - "Coordinate Conversions and Transformations including Formulas",
* EPSG Guidence Note Number 7, Version 19.
* - Krakiwsky, E.J., D.B. Thomson, and R.R. Steeves. 1977. A Manual
* For Geodetic Coordinate Transformations in the Maritimes.
* Geodesy and Geomatics Engineering, UNB. Technical Report No. 48.
* - Thomson, D.B., M.P. Mepham and R.R. Steeves. 1977.
* The Stereographic Double Projection.
* Surveying Engineering, University of New Brunswick. Technical Report No. 46.
*
*
* @since 2.4
*
*
* @source $URL$
* @version $Id$
* @author Gerald I. Evenden (for original code in Proj4)
* @author Rueben Schulz
*/
public class ObliqueStereographic extends StereographicUSGS {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = -1454098847621943639L;
/*
* The tolerance used for the inverse iteration. This is smaller
* than the tolerance in the {@code StereographicOblique} superclass.
*/
private static final double ITERATION_TOLERANCE = 1E-14;
/**
* Maximum number of iterations for iterative computations.
*/
private static final int MAXIMUM_ITERATIONS = 15;
/**
* Maximum difference allowed when comparing real numbers.
*/
private static final double EPSILON = 1E-6;
/*
* Contstants used in the forward and inverse gauss methods.
*/
private final double C, K, ratexp;
/*
* Constants for the EPSG stereographic transform.
*/
private final double phic0, cosc0, sinc0, R2;
/**
* Constructs an oblique stereographic projection (EPSG equations).
*
* @param parameters The group of parameter values.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
protected ObliqueStereographic(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
this(parameters, Provider.PARAMETERS);
}
/**
* Constructs an oblique stereographic projection (EPSG equations).
*
* @param parameters The group of parameter values.
* @param descriptor The expected parameter descriptor.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
ObliqueStereographic(final ParameterValueGroup parameters,
final ParameterDescriptorGroup descriptor)
throws ParameterNotFoundException
{
super(parameters, descriptor);
// Compute constants
final double sphi = sin(latitudeOfOrigin);
double cphi = cos(latitudeOfOrigin);
cphi *= cphi;
R2 = 2.0 * sqrt(1-excentricitySquared) / (1-excentricitySquared * sphi * sphi);
C = sqrt(1. + excentricitySquared * cphi * cphi / (1. - excentricitySquared));
phic0 = asin(sphi / C);
sinc0 = sin(phic0);
cosc0 = cos(phic0);
ratexp = 0.5 * C * excentricity;
K = tan(0.5 * phic0 + PI/4) /
(pow(tan(0.5 * latitudeOfOrigin + PI/4), C) * srat(excentricity * sphi, ratexp));
}
/**
* Transforms the specified (λ,φ) 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 USGS formulas, for comparaison later.
assert (ptDst = super.transformNormalized(x, y, ptDst)) != null;
y = 2.0 * atan(K * pow(tan(0.5 * y + PI/4), C) * srat(excentricity * sin(y), ratexp)) - PI/2;
x *= C;
final double sinc = sin(y);
final double cosc = cos(y);
final double cosl = cos(x);
final double k = R2 / (1.0 + sinc0 * sinc + cosc0 * cosc * cosl);
x = k * cosc * sin(x);
y = k * (cosc0 * sinc - sinc0 * cosc * cosl);
assert checkTransform(x, y, ptDst, 0.1);
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Transforms the specified (x,y) coordinates
* and stores the result in {@code ptDst}.
*/
@Override
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
// Compute using USGS formulas, for comparaison later.
assert (ptDst = super.inverseTransformNormalized(x, y, ptDst)) != null;
final double rho = hypot(x, y);
if (abs(rho) < EPSILON) {
x = 0.0;
y = phic0;
} else {
final double ce = 2.0 * atan2(rho, R2);
final double sinc = sin(ce);
final double cosc = cos(ce);
x = atan2(x * sinc, rho * cosc0 * cosc - y * sinc0 * sinc);
y = (cosc * sinc0) + (y * sinc * cosc0 / rho);
if (abs(y) >= 1.0) {
y = (y < 0.0) ? -PI/2.0 : PI/2.0;
} else {
y = asin(y);
}
}
// Begin pj_inv_gauss(...) method inlined
x /= C;
double num = pow(tan(0.5 * y + PI/4)/K, 1.0/C);
for (int i=MAXIMUM_ITERATIONS;;) {
double phi = 2.0 * atan(num * srat(excentricity * sin(y), -0.5 * excentricity)) - PI/2;
if (abs(phi - y) < ITERATION_TOLERANCE) {
break;
}
y = phi;
if (--i < 0) {
throw new ProjectionException(ErrorKeys.NO_CONVERGENCE);
}
}
// End pj_inv_gauss(...) method inlined
// TODO: the tolerance in the following assertion is quite large for
// an angle in radians. We should check if this is normal.
assert checkInverseTransform(x, y, ptDst, 0.01);
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* A simple function used by the transforms.
*/
private static double srat(double esinp, double exp) {
return pow((1.0 - esinp) / (1.0 + esinp), exp);
}
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////// ////////
//////// PROVIDERS ////////
//////// ////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
/**
* The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
* provider} for a stereographic projection of any kind. The equations used are the one from
* EPSG.
*
* @since 2.4
* @source $URL$
* @version $Id$
* @author Rueben Schulz
*
* @see org.geotools.referencing.operation.DefaultMathTransformFactory
*/
public static final class Provider extends Stereographic.Provider {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = 6505988910141381354L;
/**
* The parameters group.
*/
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "Oblique_Stereographic"),
new NamedIdentifier(Citations.EPSG, "Oblique Stereographic"),
new NamedIdentifier(Citations.EPSG, "Roussilhe"),
new NamedIdentifier(Citations.EPSG, "9809"),
new NamedIdentifier(Citations.GEOTIFF, "CT_ObliqueStereographic"),
new NamedIdentifier(Citations.ESRI, "Double_Stereographic"),
new NamedIdentifier(Citations.GEOTOOLS, NAME)
}, new ParameterDescriptor[] {
SEMI_MAJOR, SEMI_MINOR,
CENTRAL_MERIDIAN, LATITUDE_OF_ORIGIN,
SCALE_FACTOR,
FALSE_EASTING, FALSE_NORTHING
});
/**
* Constructs a new provider.
*/
public Provider() {
super(PARAMETERS);
}
/**
* Creates the general case.
*/
@Override
MathTransform createMathTransform(final ParameterValueGroup parameters,
final ParameterDescriptorGroup descriptor)
throws ParameterNotFoundException
{
return new ObliqueStereographic(parameters, descriptor);
}
}
}