/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 2006-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 java.util.Collection;
import javax.measure.unit.NonSI;
import org.opengis.parameter.GeneralParameterDescriptor;
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.metadata.iso.citation.Citations;
import org.geotools.referencing.NamedIdentifier;
import org.geotools.resources.i18n.ErrorKeys;
import static java.lang.Math.*;
/**
* Lambert Azimuthal Equal Area (EPSG code 9820).
* <p>
* <b>References:</b>
* <ul>
* <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
* <li> John P. Snyder (Map Projections - A Working Manual,
* U.S. Geological Survey Professional Paper 1395)</li>
* </ul>
*
* @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection on MathWorld</A>
* @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area" on RemoteSensing.org</A>
*
* @since 2.4
* @version $Id$
* @source $URL$
* @author Gerald Evenden (for original code in Proj4)
* @author Beate Stollberg
* @author Martin Desruisseaux
*/
public class LambertAzimuthalEqualArea extends MapProjection {
/** For cross-version compatibility. */
private static final long serialVersionUID = 1639914708790574760L;
/** Maximum difference allowed when comparing real numbers. */
private static final double EPSILON = 1E-7;
/** Epsilon for the comparaison of small quantities. */
private static final double FINE_EPSILON = 1E-10;
/** Epsilon for the comparaison of latitudes. */
private static final double EPSILON_LATITUDE = 1E-10;
/** Constants for authalic latitude. */
private static final double P00 = 0.33333333333333333333,
P01 = 0.17222222222222222222,
P02 = 0.10257936507936507936,
P10 = 0.06388888888888888888,
P11 = 0.06640211640211640211,
P20 = 0.01641501294219154443;
/** The projection mode. */
static final int OBLIQUE=0, EQUATORIAL=1, NORTH_POLE=2, SOUTH_POLE=3;
/** The projection mode for this particular instance. */
final int mode;
/** Constant parameters. */
final double sinb1, cosb1, xmf, ymf, mmf, qp, dd, rq;
/** Coefficients for authalic latitude. */
private final double APA0, APA1, APA2;
/**
* Constructs a new map projection from the supplied parameters.
*
* @param parameters The parameter values in standard units.
* @throws ParameterNotFoundException if a mandatory parameter is missing.
*/
protected LambertAzimuthalEqualArea(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
// Fetch parameters
super(parameters);
final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
latitudeOfOrigin = doubleValue(expected, Provider.LATITUDE_OF_CENTRE, parameters);
centralMeridian = doubleValue(expected, Provider.LONGITUDE_OF_CENTRE, parameters);
ensureLatitudeInRange (Provider.LATITUDE_OF_CENTRE, latitudeOfOrigin, true);
ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE, centralMeridian, true);
/*
* Detects the mode (oblique, etc.).
*/
final double t = abs(latitudeOfOrigin);
if (abs(t - PI/2) < EPSILON_LATITUDE) {
mode = latitudeOfOrigin < 0.0 ? SOUTH_POLE : NORTH_POLE;
} else if (abs(t) < EPSILON_LATITUDE) {
mode = EQUATORIAL;
} else {
mode = OBLIQUE;
}
/*
* Computes the constants for authalic latitude.
*/
final double es2 = excentricitySquared * excentricitySquared;
final double es3 = excentricitySquared * es2;
APA0 = P02 * es3 + P01 * es2 + P00 * excentricitySquared;
APA1 = P11 * es3 + P10 * es2;
APA2 = P20 * es3;
final double sinphi;
qp = qsfn(1);
rq = sqrt(0.5 * qp);
mmf = 0.5 / (1 - excentricitySquared);
sinphi = sin(latitudeOfOrigin);
if (isSpherical) {
sinb1 = sin(latitudeOfOrigin);
cosb1 = cos(latitudeOfOrigin);
} else {
sinb1 = qsfn(sinphi) / qp;
cosb1 = sqrt(1.0 - sinb1 * sinb1);
}
switch (mode) {
case NORTH_POLE: // Fall through
case SOUTH_POLE: {
dd = 1.0;
xmf = ymf = rq;
break;
}
case EQUATORIAL: {
dd = 1.0 / rq;
xmf = 1.0;
ymf = 0.5 * qp;
break;
}
case OBLIQUE: {
dd = cos(latitudeOfOrigin) /
(sqrt(1.0 - excentricitySquared * sinphi * sinphi) * rq * cosb1);
xmf = rq * dd;
ymf = rq / dd;
break;
}
default: {
throw new AssertionError(mode);
}
}
}
/**
* {@inheritDoc}
*/
public ParameterDescriptorGroup getParameterDescriptors() {
return Provider.PARAMETERS;
}
/**
* {@inheritDoc}
*/
@Override
public ParameterValueGroup getParameterValues() {
final ParameterValueGroup values = super.getParameterValues();
final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
set(expected, Provider.LATITUDE_OF_CENTRE, values, latitudeOfOrigin);
set(expected, Provider.LONGITUDE_OF_CENTRE, values, centralMeridian);
return values;
}
/**
* 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(final double lambda, final double phi, Point2D ptDst)
throws ProjectionException
{
final double coslam = cos(lambda);
final double sinlam = sin(lambda);
final double sinphi = sin(phi);
double q = qsfn(sinphi);
final double sinb, cosb, b, c, x, y;
switch (mode) {
case OBLIQUE: {
sinb = q / qp;
cosb = sqrt(1.0 - sinb * sinb);
c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
b = sqrt(2.0 / c);
y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
x = xmf * b * cosb * sinlam;
break;
}
case EQUATORIAL: {
sinb = q / qp;
cosb = sqrt(1.0 - sinb * sinb);
c = 1.0 + cosb * coslam;
b = sqrt(2.0 / c);
y = ymf * b * sinb;
x = xmf * b * cosb * sinlam;
break;
}
case NORTH_POLE: {
c = (PI / 2) + phi;
q = qp - q;
if (q >= 0.0) {
b = sqrt(q);
x = b * sinlam;
y = coslam * -b;
} else {
x = y = 0.;
}
break;
}
case SOUTH_POLE: {
c = phi - (PI / 2);
q = qp + q;
if (q >= 0.0) {
b = sqrt(q);
x = b * sinlam;
y = coslam * +b;
} else {
x = y = 0.;
}
break;
}
default: {
throw new AssertionError(mode);
}
}
if (abs(c) < EPSILON_LATITUDE) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Transforms the specified (<var>x</var>,<var>y</var>) coordinate
* and stores the result in {@code ptDst}.
*/
@Override
@SuppressWarnings("fallthrough")
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
final double lambda, phi;
switch (mode) {
case EQUATORIAL: // Fall through
case OBLIQUE: {
x /= dd;
y *= dd;
final double rho = hypot(x, y);
if (rho < FINE_EPSILON) {
lambda = 0.0;
phi = latitudeOfOrigin;
} else {
double sCe, cCe, q, ab;
sCe = 2.0 * asin(0.5 * rho / rq);
cCe = cos(sCe);
sCe = sin(sCe);
x *= sCe;
if (mode == OBLIQUE) {
ab = cCe * sinb1 + y * sCe * cosb1 / rho;
q = qp * ab;
y = rho * cosb1 * cCe - y * sinb1 * sCe;
} else {
ab = y * sCe / rho;
q = qp * ab;
y = rho * cCe;
}
lambda = atan2(x, y);
phi = authlat(asin(ab));
}
break;
}
case NORTH_POLE: {
y = -y;
// Fall through
}
case SOUTH_POLE: {
final double q = x*x + y*y;
if (q == 0) {
lambda = 0.;
phi = latitudeOfOrigin;
} else {
double ab = 1.0 - q / qp;
if (mode == SOUTH_POLE) {
ab = -ab;
}
lambda = atan2(x, y);
phi = authlat(asin(ab));
}
break;
}
default: {
throw new AssertionError(mode);
}
}
if (ptDst != null) {
ptDst.setLocation(lambda, phi);
return ptDst;
}
return new Point2D.Double(lambda, phi);
}
/**
* Provides the transform equations for the spherical case.
*
* @version $Id$
* @author Martin Desruisseaux
*/
private static final class Spherical extends LambertAzimuthalEqualArea {
/**
* For cross-version compatibility.
*/
private static final long serialVersionUID = 2091431369806844342L;
/**
* Constructs a new map projection from the suplied parameters.
*
* @param parameters The parameter values in standard units.
* @throws ParameterNotFoundException if a mandatory parameter is missing.
*/
protected Spherical(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
super(parameters);
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(final double lambda, final double phi, Point2D ptDst)
throws ProjectionException
{
// Compute using ellipsoidal formulas, for comparaison later.
assert (ptDst = super.transformNormalized(lambda, phi, ptDst)) != null;
final double sinphi = sin(phi);
final double cosphi = cos(phi);
final double coslam = cos(lambda);
double x,y;
switch (mode) {
case EQUATORIAL: {
y = 1.0 + cosphi * coslam;
if (y <= FINE_EPSILON) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
y = sqrt(2.0 / y);
x = y * cosphi * sin(lambda);
y *= sinphi;
break;
}
case OBLIQUE: {
y = 1.0 + sinb1 * sinphi + cosb1 * cosphi * coslam;
if (y <= FINE_EPSILON) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
y = sqrt(2.0 / y);
x = y * cosphi * sin(lambda);
y *= cosb1 * sinphi - sinb1 * cosphi * coslam;
break;
}
case NORTH_POLE: {
if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
y = (PI/4) - phi * 0.5;
y = 2.0 * sin(y);
x = y * sin(lambda);
y *= -coslam;
break;
}
case SOUTH_POLE: {
if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
y = (PI/4) - phi * 0.5;
y = 2.0 * cos(y);
x = y * sin(lambda);
y *= +coslam;
break;
}
default: {
throw new AssertionError(mode);
}
}
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>) coordinate
* and stores the result in {@code ptDst} using equations for a sphere.
*/
@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;
double lambda, phi;
final double rh = hypot(x, y);
phi = rh * 0.5;
if (phi > 1.0) {
throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
}
phi = 2.0 * asin(phi);
switch (mode) {
case EQUATORIAL: {
final double sinz = sin(phi);
final double cosz = cos(phi);
phi = abs(rh) <= FINE_EPSILON ? 0.0 : asin(y * sinz / rh);
x *= sinz;
y = cosz * rh;
lambda = (y == 0) ? 0.0 : atan2(x, y);
break;
}
case OBLIQUE: {
final double sinz = sin(phi);
final double cosz = cos(phi);
phi = abs(rh) <= FINE_EPSILON ? latitudeOfOrigin :
asin(cosz * sinb1 + y * sinz * cosb1 / rh);
x *= sinz * cosb1;
y = (cosz - sin(phi) * sinb1) * rh;
lambda = (y == 0) ? 0.0 : atan2(x, y);
break;
}
case NORTH_POLE: {
phi = (PI / 2) - phi;
lambda = atan2(x, -y);
break;
}
case SOUTH_POLE: {
phi -= (PI / 2);
lambda = atan2(x, y);
break;
}
default: {
throw new AssertionError(mode);
}
}
assert checkInverseTransform(lambda, phi, ptDst);
if (ptDst != null) {
ptDst.setLocation(lambda, phi);
return ptDst;
}
return new Point2D.Double(lambda, phi);
}
}
/**
* Calculates <var>q</var>, Snyder equation (3-12)
*
* @param sinphi sin of the latitude <var>q</var> is calculated for.
* @return <var>q</var> from Snyder equation (3-12).
*/
private double qsfn(final double sinphi) {
if (excentricity >= EPSILON) {
final double con = excentricity * sinphi;
return ((1.0 - excentricitySquared) * (sinphi / (1.0 - con*con) -
(0.5 / excentricity) * log((1.0 - con) / (1.0 + con))));
} else {
return sinphi + sinphi;
}
}
/**
* Determines latitude from authalic latitude.
*/
private double authlat(final double beta) {
final double t = beta + beta;
return beta + APA0 * sin(t) + APA1 * sin(t+t) + APA2 * sin(t+t+t);
}
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////// ////////
//////// PROVIDERS ////////
//////// ////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
/**
* The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
* provider} for an {@linkplain LambertAzimuthalEqualArea Lambert Equal Area} projection
* (EPSG code 9820).
*
* @since 2.4
* @version $Id$
* @author Beate Stollberg
*
* @see org.geotools.referencing.operation.DefaultMathTransformFactory
*/
public static class Provider extends AbstractProvider {
/**
* For cross-version compatibility.
*/
private static final long serialVersionUID = 3877793025552244132L;
/**
* The operation parameter descriptor for the {@link #latitudeOfOrigin}
* parameter value. Valid values range is from -90 to 90°. Default value is 0.
*/
public static final ParameterDescriptor LATITUDE_OF_CENTRE = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "latitude_of_center"),
new NamedIdentifier(Citations.EPSG, "Latitude of natural origin"),
new NamedIdentifier(Citations.EPSG, "Spherical latitude of origin"),
new NamedIdentifier(Citations.ESRI, "latitude_of_origin"),
new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLat")
},
0, -90, 90, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@link #centralMeridian}
* parameter value. Valid values range is from -180 to 180°. Default value is 0.
*/
public static final ParameterDescriptor LONGITUDE_OF_CENTRE = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "longitude_of_center"),
new NamedIdentifier(Citations.EPSG, "Longitude of natural origin"),
new NamedIdentifier(Citations.EPSG, "Spherical longitude of origin"),
new NamedIdentifier(Citations.ESRI, "central_meridian"),
new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLong")
},
0, -180, 180, NonSI.DEGREE_ANGLE);
/**
* The parameters group.
*/
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "Lambert_Azimuthal_Equal_Area"),
new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area"),
new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area (Spherical)"),
new NamedIdentifier(Citations.GEOTIFF, "CT_LambertAzimEqualArea"),
new NamedIdentifier(Citations.EPSG, "9820"),
}, new ParameterDescriptor[] {
SEMI_MAJOR, SEMI_MINOR,
LATITUDE_OF_CENTRE, LONGITUDE_OF_CENTRE,
FALSE_EASTING, FALSE_NORTHING
});
/**
* Constructs a new provider.
*/
public Provider() {
super(PARAMETERS);
}
/**
* Creates a transform from the specified group of parameter values.
*
* @param parameters The group of parameter values.
* @return The created math transform.
* @throws ParameterNotFoundException if a required parameter was not found.
*/
public MathTransform createMathTransform(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
return isSpherical(parameters) ? new Spherical(parameters) :
new LambertAzimuthalEqualArea(parameters);
}
}
}