/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 2005-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.util.InternationalString;
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.parameter.InvalidParameterValueException;
import org.opengis.referencing.operation.CylindricalProjection;
import org.opengis.referencing.operation.MathTransform;
import org.geotools.measure.Angle;
import org.geotools.measure.Latitude;
import org.geotools.metadata.iso.citation.Citations;
import org.geotools.referencing.NamedIdentifier;
import org.geotools.resources.i18n.Vocabulary;
import org.geotools.resources.i18n.VocabularyKeys;
import org.geotools.resources.i18n.Errors;
import org.geotools.resources.i18n.ErrorKeys;
import static java.lang.Math.*;
/**
* Oblique Mercator Projection. A conformal, oblique, cylindrical projection with the cylinder
* touching the ellipsoid (or sphere) along a great circle path (the central line). The
* {@linkplain Mercator} and {@linkplain TransverseMercator Transverse Mercator} projections can
* be thought of as special cases of the oblique mercator, where the central line is along the
* equator or a meridian, respectively. The Oblique Mercator projection has been used in
* Switzerland, Hungary, Madagascar, Malaysia, Borneo and the panhandle of Alaska.
* <p>
* The Oblique Mercator projection uses a (<var>U</var>,<var>V</var>) coordinate system, with the
* <var>U</var> axis along the central line. During the forward projection, coordinates from the
* ellipsoid are projected conformally to a sphere of constant total curvature, called the
* "aposphere", before being projected onto the plane. The projection coordinates are further
* convented to a (<var>X</var>,<var>Y</var>) coordinate system by rotating the calculated
* (<var>u</var>,<var>v</var>) coordinates to give output (<var>x</var>,<var>y</var>) coordinates.
* The rotation value is usually the same as the projection azimuth (the angle, east of north, of
* the central line), but some cases allow a separate rotation parameter.
* <p>
* There are two forms of the oblique mercator, differing in the origin of their grid coordinates.
* The {@linkplain HotineObliqueMercator Hotine Oblique Mercator} (EPSG code 9812) has grid
* coordinates start at the intersection of the central line and the equator of the aposphere.
* The {@linkplain ObliqueMercator Oblique Mercator} (EPSG code 9815) is the same, except the
* grid coordinates begin at the central point (where the latitude of center and central line
* intersect). ESRI separates these two case by appending {@code "Natural_Origin"} (for the
* {@code "Hotine_Oblique_Mercator"}) and {@code "Center"} (for the {@code "Oblique_Mercator"})
* to the projection names.
* <p>
* Two different methods are used to specify the central line for the oblique mercator:
* 1) a central point and an azimuth, east of north, describing the central line and
* 2) two points on the central line. The EPSG does not use the two point method,
* while ESRI separates the two cases by putting {@code "Azimuth"} and {@code "Two_Point"}
* in their projection names. Both cases use the point where the {@code "latitude_of_center"}
* parameter crosses the central line as the projection's central point.
* The {@linkplain #centralMeridian central meridian} is not a projection parameter,
* and is instead calculated as the intersection between the central line and the
* equator of the aposphere.
* <p>
* For the azimuth method, the central latitude cannot be ±90.0 degrees
* and the central line cannot be at a maximum or minimum latitude at the central point.
* In the two point method, the latitude of the first and second points cannot be
* equal. Also, the latitude of the first point and central point cannot be
* ±90.0 degrees. Furthermore, the latitude of the first point cannot be 0.0 and
* the latitude of the second point cannot be -90.0 degrees. A change of
* 10<sup>-7</sup> radians can allow calculation at these special cases. Snyder's restriction
* of the central latitude being 0.0 has been removed, since the equations appear
* to work correctly in this case.
* <p>
* Azimuth values of 0.0 and ±90.0 degrees are allowed (and used in Hungary
* and Switzerland), though these cases would usually use a Mercator or
* Transverse Mercator projection instead. Azimuth values > 90 degrees cause
* errors in the equations.
* <p>
* The oblique mercator is also called the "Rectified Skew Orthomorphic" (RSO). It appears
* is that the only difference from the oblique mercator is that the RSO allows the rotation
* from the (<var>U</var>,<var>V</var>) to (<var>X</var>,<var>Y</var>) coordinate system to
* be different from the azimuth. This separate parameter is called
* {@code "rectified_grid_angle"} (or {@code "XY_Plane_Rotation"} by ESRI) and is also
* included in the EPSG's parameters for the Oblique Mercator and Hotine Oblique Mercator.
* The rotation parameter is optional in all the non-two point projections and will be
* set to the azimuth if not specified.
* <p>
* Projection cases and aliases implemented by the {@link ObliqueMercator} are:
* <ul>
* <li>{@code Oblique_Mercator} (EPSG code 9815)<br>
* grid coordinates begin at the central point,
* has {@code "rectified_grid_angle"} parameter.</li>
* <li>{@code Hotine_Oblique_Mercator_Azimuth_Center} (ESRI)<br>
* grid coordinates begin at the central point.</li>
* <li>{@code Rectified_Skew_Orthomorphic_Center} (ESRI)<br>
* grid coordinates begin at the central point,
* has {@code "rectified_grid_angle"} parameter.</li>
*
* <li>{@code Hotine_Oblique_Mercator} (EPSG code 9812)<br>
* grid coordinates begin at the interseciton of the central line and aposphere equator,
* has {@code "rectified_grid_angle"} parameter.</li>
* <li>{@code Hotine_Oblique_Mercator_Azimuth_Natural_Origin} (ESRI)<br>
* grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
* <li>{@code Rectified_Skew_Orthomorphic_Natural_Origin} (ESRI)<br>
* grid coordinates begin at the interseciton of the central line and aposphere equator,
* has {@code "rectified_grid_angle"} parameter.</li>
*
* <li>{@code Hotine_Oblique_Mercator_Two_Point_Center} (ESRI)<br>
* grid coordinates begin at the central point.</li>
* <li>{@code Hotine_Oblique_Mercator_Two_Point_Natural_Origin} (ESRI)<br>
* grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
* </ul>
* <p>
* <b>References:</b>
* <ul>
* <li>{@code libproj4} is available at
* <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/">libproj4 Miscellanea</A><br>
* Relevent files are: {@code PJ_omerc.c}, {@code pj_tsfn.c},
* {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}</li>
* <li>John P. Snyder (Map Projections - A Working Manual,
* U.S. Geological Survey Professional Paper 1395, 1987)</li>
* <li>"Coordinate Conversions and Transformations including Formulas",
* EPSG Guidence Note Number 7 part 2, Version 24.</li>
* <li>Gerald Evenden, 2004, <a href="http://members.verizon.net/~vze2hc4d/proj4/omerc.pdf">
* Documentation of revised Oblique Mercator</a></li>
* </ul>
*
* @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Oblique Mercator projection on MathWorld</A>
* @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/hotine_oblique_mercator.html">"hotine_oblique_mercator" on RemoteSensing.org</A>
* @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_mercator.html">"oblique_mercator" on RemoteSensing.org</A>
*
* @since 2.1
*
*
* @source $URL$
* @version $Id$
* @author Gerald I. Evenden (for original code in Proj4)
* @author Rueben Schulz
*/
public class ObliqueMercator extends MapProjection {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = 5382294977124711214L;
/**
* Maximum difference allowed when comparing real numbers.
*/
private static final double EPSILON = 1E-6;
/**
* Maximum difference allowed when comparing latitudes.
*/
private static final double EPSILON_LATITUDE = 1E-10;
//////
////// Map projection parameters. The following are NOT used by the transformation
////// methods, but are stored in order to restitute them in WKT formatting. They
////// are made visible ('protected' access) for documentation purpose and because
////// they are user-supplied parameters, not derived coefficients.
//////
/**
* Latitude of the projection centre. This is similar to the
* {@link #latitudeOfOrigin}, but the latitude of origin is the
* Earth equator on aposphere for the oblique mercator.
*/
protected final double latitudeOfCentre;
/**
* Longitude of the projection centre. This is <strong>NOT</strong> equal
* to the {@link #centralMeridian}, which is the meridian where the
* central line intersects the Earth equator on aposphere.
* <p>
* This parameter applies to the "azimuth" case only.
* It is set to {@link Double#NaN NaN} for the "two points" case.
*/
protected final double longitudeOfCentre;
/**
* The azimuth of the central line passing throught the centre of the projection, in radians.
*/
protected final double azimuth;
/**
* The rectified bearing of the central line, in radians. This is equals to the
* {@linkplain #azimuth} if the {@link Provider#RECTIFIED_GRID_ANGLE "rectified_grid_angle"}
* parameter value is not set.
*/
protected final double rectifiedGridAngle;
// The next parameters still private for now because I'm not
// sure if they should appear in some 'TwoPoint' subclass...
/**
* The latitude of the 1st point used to specify the central line, in radians.
* <p>
* This parameter applies to the "two points" case only.
* It is set to {@link Double#NaN NaN} for the "azimuth" case.
*/
private final double latitudeOf1stPoint;
/**
* The longitude of the 1st point used to specify the central line, in radians.
* <p>
* This parameter applies to the "two points" case only.
* It is set to {@link Double#NaN NaN} for the "azimuth" case.
*/
private final double longitudeOf1stPoint;
/**
* The latitude of the 2nd point used to specify the central line, in radians.
* <p>
* This parameter applies to the "two points" case only.
* It is set to {@link Double#NaN NaN} for the "azimuth" case.
*/
private final double latitudeOf2ndPoint;
/**
* The longitude of the 2nd point used to specify the central line, in radians.
* <p>
* This parameter applies to the "two points" case only.
* It is set to {@link Double#NaN NaN} for the "azimuth" case.
*/
private final double longitudeOf2ndPoint;
//////
////// Map projection coefficients computed from the above parameters.
////// They are the fields used for coordinate transformations.
//////
/**
* Constants used in the transformation.
*/
private final double B, A, E;
/**
* Convenience values equal to {@link #A} / {@link #B},
* {@link #A}×{@link #B}, and {@link #B} / {@link #A}.
*/
private final double ArB, AB, BrA;
/**
* <var>v</var> values when the input latitude is a pole.
*/
private final double v_pole_n, v_pole_s;
/**
* Sine and Cosine values for gamma0 (the angle between the meridian
* and central line at the intersection between the central line and
* the Earth equator on aposphere).
*/
private final double singamma0, cosgamma0;
/**
* Sine and Cosine values for the rotation between (U,V) and
* (X,Y) coordinate systems
*/
private final double sinrot, cosrot;
/**
* <var>u</var> value (in (U,V) coordinate system) of the central point. Used in
* the oblique mercator case. The <var>v</var> value of the central point is 0.0.
*/
private final double u_c;
/**
* {@code true} if using two points on the central line to specify the azimuth.
*/
final boolean twoPoint;
/**
* 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.
*
* @since 2.4
*
* @todo Current implementation assumes a "azimuth" case. We may try to detect the
* "two points" case in a future version if needed.
*/
protected ObliqueMercator(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
this(parameters, Provider.PARAMETERS.descriptors(), false, false);
}
/**
* Constructs a new map projection from the supplied parameters.
*
* @param parameters The parameter values in standard units.
* @param expected The expected parameter descriptors.
* @param twoPoint {@code true} for the "two points" case, or {@code false} for the
* "azimuth" case. The former is used by ESRI but not EPSG.
* @param hotine {@code true} only if invoked by the {@link HotineObliqueMercator}
* constructor.
* @throws ParameterNotFoundException if a mandatory parameter is missing.
*/
ObliqueMercator(final ParameterValueGroup parameters,
final Collection<GeneralParameterDescriptor> expected,
final boolean twoPoint, final boolean hotine)
throws ParameterNotFoundException
{
// Fetch parameters
super(parameters, expected);
this.twoPoint = twoPoint;
/*
* Sets the 'centralMeridian' and 'latitudeOfOrigin' fields to NaN for safety
* (they are not 'ObliqueMercator' parameters, so the super-class initialized
* them to 0.0 by default). The 'centralMeridian' value will be computed later.
*/
centralMeridian = Double.NaN;
latitudeOfOrigin = Double.NaN;
latitudeOfCentre = doubleValue(expected, Provider.LATITUDE_OF_CENTRE, parameters);
/*
* Checks that 'latitudeOfCentre' is not +- 90 degrees.
* Not checking if 'latitudeOfCentere' is 0, since equations behave correctly.
*/
ensureLatitudeInRange(Provider.LATITUDE_OF_CENTRE, latitudeOfCentre, false);
/*
* Computes common constants now. In a previous version of 'ObliqueMercator', those
* constants were computed a little bit later. We compute them now in order to have
* a single 'if (twoPoint) ... else' statement, which help us to keep every fields
* final and catch some potential errors at compile-time (e.g. unitialized fields).
*/
final double com = sqrt(1.0 - excentricitySquared);
final double sinphi0 = sin(latitudeOfCentre);
final double cosphi0 = cos(latitudeOfCentre);
double temp = cosphi0 * cosphi0;
B = sqrt(1.0 + excentricitySquared * (temp * temp) / (1.0 - excentricitySquared));
final double con = 1.0 - excentricitySquared * sinphi0 * sinphi0;
A = B * com / con;
final double D = B * com / (cosphi0 * sqrt(con));
double F = D * D - 1.0;
if (F < 0.0) {
F = 0.0;
} else {
F = sqrt(F);
if (latitudeOfCentre < 0.0) { // Taking sign of 'latitudeOfCentre'
F = -F;
}
}
F = F += D;
E = F * pow(tsfn(latitudeOfCentre, sinphi0), B);
/*
* Computes the constants that depend on the "twoPoint" vs "azimuth" case. In the
* two points case, we compute them from (LAT_OF_1ST_POINT, LONG_OF_1ST_POINT) and
* (LAT_OF_2ND_POINT, LONG_OF_2ND_POINT). For the "azimuth" case, we compute them
* from LONGITUDE_OF_CENTRE and AZIMUTH.
*/
final double gamma0;
if (twoPoint) {
longitudeOfCentre = Double.NaN; // This is for the "azimuth" case only.
latitudeOf1stPoint = doubleValue(expected, Provider_TwoPoint.LAT_OF_1ST_POINT, parameters);
// Checks that latOf1stPoint is not +-90 degrees
ensureLatitudeInRange(Provider_TwoPoint.LAT_OF_1ST_POINT, latitudeOf1stPoint, false);
longitudeOf1stPoint = doubleValue(expected, Provider_TwoPoint.LONG_OF_1ST_POINT, parameters);
ensureLongitudeInRange(Provider_TwoPoint.LONG_OF_1ST_POINT, longitudeOf1stPoint, true);
latitudeOf2ndPoint = doubleValue(expected, Provider_TwoPoint.LAT_OF_2ND_POINT, parameters);
ensureLatitudeInRange(Provider_TwoPoint.LAT_OF_2ND_POINT, latitudeOf2ndPoint, true);
double longitudeOf2ndPoint; // Will be assigned to the field later.
longitudeOf2ndPoint = doubleValue(expected, Provider_TwoPoint.LONG_OF_2ND_POINT, parameters);
ensureLongitudeInRange(Provider_TwoPoint.LONG_OF_2ND_POINT, longitudeOf2ndPoint, true);
/*
* Ensures that (phi1 != phi2), (phi1 != 0°) and (phi2 != -90°),
* as specified in class javadoc.
*/
ParameterDescriptor desc = null;
Object value = null;
if (abs(latitudeOf1stPoint - latitudeOf2ndPoint) < EPSILON_LATITUDE) {
desc = Provider_TwoPoint.LAT_OF_1ST_POINT;
value = Provider_TwoPoint.LAT_OF_2ND_POINT.getName().getCode();
// Exception will be thrown below.
}
if (abs(latitudeOf1stPoint) < EPSILON_LATITUDE) {
desc = Provider_TwoPoint.LAT_OF_1ST_POINT;
value = new Latitude(latitudeOf1stPoint);
// Exception will be thrown below.
}
if (abs(latitudeOf2ndPoint + PI/2.0) < EPSILON_LATITUDE) {
desc = Provider_TwoPoint.LAT_OF_2ND_POINT;
value = new Latitude(latitudeOf2ndPoint);
// Exception will be thrown below.
}
if (desc != null) {
final String name = desc.getName().getCode();
throw new InvalidParameterValueException(Errors.format(
ErrorKeys.ILLEGAL_ARGUMENT_$2, name, value), name, value);
}
/*
* The coefficients for the "two points" case.
*/
final double H = pow(tsfn(latitudeOf1stPoint, sin(latitudeOf1stPoint)), B);
final double L = pow(tsfn(latitudeOf2ndPoint, sin(latitudeOf2ndPoint)), B);
final double Fp = E / H;
final double P = (L - H) / (L + H);
double J = E * E;
J = (J - L * H) / (J + L * H);
double diff = longitudeOf1stPoint - longitudeOf2ndPoint;
if (diff < -PI) {
longitudeOf2ndPoint -= 2.0 * PI;
} else if (diff > PI) {
longitudeOf2ndPoint += 2.0 * PI;
}
this.longitudeOf2ndPoint = longitudeOf2ndPoint;
centralMeridian = rollLongitude(0.5 * (longitudeOf1stPoint + longitudeOf2ndPoint) -
atan(J * tan(0.5 * B * (longitudeOf1stPoint - longitudeOf2ndPoint)) / P) / B);
gamma0 = atan(2.0 * sin(B * rollLongitude(longitudeOf1stPoint - centralMeridian)) /
(Fp - 1.0 / Fp));
azimuth = asin(D * sin(gamma0));
rectifiedGridAngle = azimuth;
} else {
/*
* Computes coefficients for the "azimuth" case. Set the 1st and 2nd points
* to (NaN,NaN) since they are specific to the "two points" case. They are
* involved in WKT formatting only, not in transformation calculation.
*/
latitudeOf1stPoint = Double.NaN;
longitudeOf1stPoint = Double.NaN;
latitudeOf2ndPoint = Double.NaN;
longitudeOf2ndPoint = Double.NaN;
longitudeOfCentre = doubleValue(expected, Provider.LONGITUDE_OF_CENTRE, parameters);
ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE, longitudeOfCentre, true);
azimuth = doubleValue(expected, Provider.AZIMUTH, parameters);
// Already checked for +-360 deg. above.
if ((azimuth > -1.5*PI && azimuth < -0.5*PI) ||
(azimuth > 0.5*PI && azimuth < 1.5*PI))
{
final String name = Provider.AZIMUTH.getName().getCode();
final Angle value = new Angle(toDegrees(azimuth));
throw new InvalidParameterValueException(Errors.format(
ErrorKeys.ILLEGAL_ARGUMENT_$2, name, value), name, value);
}
temp = doubleValue(expected, Provider.RECTIFIED_GRID_ANGLE, parameters);
if (Double.isNaN(temp)) {
temp = azimuth;
}
rectifiedGridAngle = temp;
gamma0 = asin(sin(azimuth) / D);
// Check for asin(+-1.00000001)
temp = 0.5 * (F - 1.0 / F) * tan(gamma0);
if (abs(temp) > 1.0) {
if (abs(abs(temp) - 1.0) > EPSILON) {
throw new IllegalArgumentException(Errors.format(ErrorKeys.TOLERANCE_ERROR));
}
temp = (temp > 0) ? 1.0 : -1.0;
}
centralMeridian = longitudeOfCentre - asin(temp) / B;
}
/*
* More coefficients common to all kind of oblique mercator.
*/
singamma0 = sin(gamma0);
cosgamma0 = cos(gamma0);
sinrot = sin(rectifiedGridAngle);
cosrot = cos(rectifiedGridAngle);
ArB = A / B;
AB = A * B;
BrA = B / A;
v_pole_n = ArB * log(tan(0.5 * (PI/2.0 - gamma0)));
v_pole_s = ArB * log(tan(0.5 * (PI/2.0 + gamma0)));
if (hotine) {
u_c = 0.0;
} else {
if (abs(abs(azimuth) - PI/2.0) < EPSILON_LATITUDE) {
// LongitudeOfCentre = NaN in twoPoint, but azimuth cannot be 90 here (lat1 != lat2)
u_c = A * (longitudeOfCentre - centralMeridian);
} else {
double u_c = abs(ArB * atan2(sqrt(D * D - 1.0), cos(azimuth)));
if (latitudeOfCentre < 0.0) {
u_c = -u_c;
}
this.u_c = u_c;
}
}
}
/**
* {@inheritDoc}
*/
public ParameterDescriptorGroup getParameterDescriptors() {
return (twoPoint) ? Provider_TwoPoint.PARAMETERS : Provider.PARAMETERS;
}
/**
* {@inheritDoc}
*/
@Override
public ParameterValueGroup getParameterValues() {
final ParameterValueGroup values = super.getParameterValues();
final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
// Note: we don't need a "if (twoPoint) ... else" statement since
// the "set" method will actually set the value only if applicable.
set(expected, Provider.LATITUDE_OF_CENTRE, values, latitudeOfCentre);
set(expected, Provider.LONGITUDE_OF_CENTRE, values, longitudeOfCentre);
set(expected, Provider.AZIMUTH, values, azimuth);
set(expected, Provider.RECTIFIED_GRID_ANGLE, values, rectifiedGridAngle);
set(expected, Provider_TwoPoint.LAT_OF_1ST_POINT, values, latitudeOf1stPoint);
set(expected, Provider_TwoPoint.LONG_OF_1ST_POINT, values, longitudeOf1stPoint);
set(expected, Provider_TwoPoint.LAT_OF_2ND_POINT, values, latitudeOf2ndPoint);
set(expected, Provider_TwoPoint.LONG_OF_2ND_POINT, values, longitudeOf2ndPoint);
return values;
}
/**
* {@inheritDoc}
*
* @param x The longitude of the coordinate, in <strong>radians</strong>.
* @param y The latitude of the coordinate, in <strong>radians</strong>.
*/
protected Point2D transformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
double u, v;
if (abs(abs(y) - PI/2.0) > EPSILON) {
double Q = E / pow(tsfn(y, sin(y)), B);
double temp = 1.0 / Q;
double S = 0.5 * (Q - temp);
double V = sin(B * x);
double U = (S * singamma0 - V * cosgamma0) / (0.5 * (Q + temp));
if (abs(abs(U) - 1.0) < EPSILON) {
throw new ProjectionException(ErrorKeys.INFINITE_VALUE_$1, "v");
}
v = 0.5 * ArB * log((1.0 - U) / (1.0 + U));
temp = cos(B * x);
if (abs(temp) < EPSILON_LATITUDE) {
u = AB * x;
} else {
u = ArB * atan2((S * cosgamma0 + V * singamma0), temp);
}
} else {
v = y > 0 ? v_pole_n : v_pole_s;
u = ArB * y;
}
u -= u_c;
x = v * cosrot + u * sinrot;
y = u * cosrot - v * sinrot;
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* {@inheritDoc}
*/
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
throws ProjectionException
{
double v = x * cosrot - y * sinrot;
double u = y * cosrot + x * sinrot + u_c;
double Qp = exp(-BrA * v);
double temp = 1.0 / Qp;
double Sp = 0.5 * (Qp - temp);
double Vp = sin(BrA * u);
double Up = (Vp * cosgamma0 + Sp * singamma0) / (0.5 * (Qp + temp));
if (abs(abs(Up) - 1.0) < EPSILON) {
x = 0.0;
y = Up < 0.0 ? -PI / 2.0 : PI / 2.0;
} else {
y = pow(E / sqrt((1. + Up) / (1. - Up)), 1.0 / B); //calculate t
y = cphi2(y);
x = -atan2((Sp * cosgamma0 - Vp * singamma0), cos(BrA * u)) / B;
}
if (ptDst != null) {
ptDst.setLocation(x,y);
return ptDst;
}
return new Point2D.Double(x,y);
}
/**
* Maximal error (in metres) tolerated for assertion, if enabled.
*
* @param longitude The longitude in decimal degrees.
* @param latitude The latitude in decimal degrees.
* @return The tolerance level for assertions, in meters.
*/
@Override
protected double getToleranceForAssertions(final double longitude, final double latitude) {
if (abs(longitude - centralMeridian)/2 + abs(latitude - latitudeOfCentre) > 10) {
// When far from the valid area, use a larger tolerance.
return 1;
}
return super.getToleranceForAssertions(longitude, latitude);
}
/**
* Returns a hash value for this projection.
*/
@Override
public int hashCode() {
long code = Double.doubleToLongBits(latitudeOfCentre);
code = code*37 + Double.doubleToLongBits(longitudeOfCentre);
code = code*37 + Double.doubleToLongBits(azimuth);
code = code*37 + Double.doubleToLongBits(rectifiedGridAngle);
code = code*37 + Double.doubleToLongBits(latitudeOf1stPoint);
code = code*37 + Double.doubleToLongBits(latitudeOf2ndPoint);
return ((int)code ^ (int)(code >>> 32)) + 37*super.hashCode();
}
/**
* Compares the specified object with this map projection for equality.
*/
@Override
public boolean equals(final Object object) {
if (object == this) {
// Slight optimization
return true;
}
if (super.equals(object)) {
final ObliqueMercator that = (ObliqueMercator) object;
return this.twoPoint == that.twoPoint &&
equals(this.latitudeOfCentre , that.latitudeOfCentre ) &&
equals(this.longitudeOfCentre , that.longitudeOfCentre ) &&
equals(this.azimuth , that.azimuth ) &&
equals(this.rectifiedGridAngle , that.rectifiedGridAngle ) &&
equals(this.latitudeOf1stPoint , that.latitudeOf1stPoint ) &&
equals(this.longitudeOf1stPoint, that.longitudeOf1stPoint) &&
equals(this.latitudeOf2ndPoint , that.latitudeOf2ndPoint ) &&
equals(this.longitudeOf2ndPoint, that.longitudeOf2ndPoint) &&
equals(this.u_c, that.u_c);
// Note: "u_c" is a derived parameter, so in theory we don't need to compare it.
// However we still compare it as a safety, because it takes a different
// value in the "hotine" case.
}
return false;
}
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////// ////////
//////// PROVIDERS ////////
//////// ////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
/**
* The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
* provider} for an {@linkplain ObliqueMercator Oblique Mercator} projection (EPSG code 9815).
*
* @since 2.1
* @version $Id$
* @author Rueben Schulz
*
* @see org.geotools.referencing.operation.DefaultMathTransformFactory
*/
public static class Provider extends AbstractProvider {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = 201776686002266891L;
/**
* The operation parameter descriptor for the {@link #latitudeOfCentre latitudeOfCentre}
* 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 projection centre"),
new NamedIdentifier(Citations.ESRI, "Latitude_Of_Center"),
new NamedIdentifier(Citations.GEOTIFF, "CenterLat")
},
0, -90, 90, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@link #longitudeOfCentre longitudeOfCentre}
* 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 projection centre"),
new NamedIdentifier(Citations.ESRI, "Longitude_Of_Center"),
new NamedIdentifier(Citations.GEOTIFF, "CenterLong")
},
0, -180, 180, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@link #azimuth azimuth}
* parameter value. Valid values range is from -360 to -270, -90 to 90,
* and 270 to 360 degrees. Default value is 0.
*/
public static final ParameterDescriptor AZIMUTH = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "azimuth"),
new NamedIdentifier(Citations.ESRI, "Azimuth"),
new NamedIdentifier(Citations.EPSG, "Azimuth of initial line"),
new NamedIdentifier(Citations.GEOTIFF, "AzimuthAngle")
},
0, -360, 360, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@link #rectifiedGridAngle
* rectifiedGridAngle} parameter value. It is an optional parameter with
* valid values ranging from -360 to 360°. Default value is {@link #azimuth azimuth}.
*/
public static final ParameterDescriptor RECTIFIED_GRID_ANGLE = createOptionalDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "rectified_grid_angle"),
new NamedIdentifier(Citations.EPSG, "Angle from Rectified to Skew Grid"),
new NamedIdentifier(Citations.ESRI, "XY_Plane_Rotation"),
new NamedIdentifier(Citations.GEOTIFF, "RectifiedGridAngle")
},
-360, 360, NonSI.DEGREE_ANGLE);
/**
* The localized name for "Oblique Mercator".
*/
static final InternationalString NAME =
Vocabulary.formatInternational(VocabularyKeys.OBLIQUE_MERCATOR_PROJECTION);
/**
* The parameters group.
*/
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
new NamedIdentifier(Citations.OGC, "Oblique_Mercator"),
new NamedIdentifier(Citations.EPSG, "Oblique Mercator"),
new NamedIdentifier(Citations.EPSG, "Hotine Oblique Mercator (variant B)"),
new NamedIdentifier(Citations.EPSG, "9815"),
new NamedIdentifier(Citations.GEOTIFF, "CT_ObliqueMercator"),
new NamedIdentifier(Citations.ESRI, "Hotine_Oblique_Mercator_Azimuth_Center"),
new NamedIdentifier(Citations.ESRI, "Rectified_Skew_Orthomorphic_Center"),
new NamedIdentifier(Citations.GEOTOOLS, NAME)
}, new ParameterDescriptor[] {
SEMI_MAJOR, SEMI_MINOR,
LONGITUDE_OF_CENTRE, LATITUDE_OF_CENTRE,
AZIMUTH, RECTIFIED_GRID_ANGLE,
SCALE_FACTOR,
FALSE_EASTING, FALSE_NORTHING
});
/**
* Constructs a new provider.
*/
public Provider() {
super(PARAMETERS);
}
/**
* Constructs a new provider.
*
* @param params A description of parameters.
*/
protected Provider(final ParameterDescriptorGroup params) {
super(params);
}
/**
* Returns the operation type for this map projection.
*/
@Override
public Class<CylindricalProjection> getOperationType() {
return CylindricalProjection.class;
}
/**
* 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.
*/
protected MathTransform createMathTransform(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
final Collection<GeneralParameterDescriptor> descriptors = PARAMETERS.descriptors();
return new ObliqueMercator(parameters, descriptors, false, false);
}
}
/**
* The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
* provider} for a {@linkplain ObliqueMercator Oblique Mercator} projection, specified
* with two points on the central line (instead of a central point and azimuth).
*
* @since 2.1
* @version $Id$
* @author Rueben Schulz
*
* @see org.geotools.referencing.operation.DefaultMathTransformFactory
*/
public static class Provider_TwoPoint extends Provider {
/**
* For compatibility with different versions during deserialization.
*/
private static final long serialVersionUID = 7124258885016543889L;
/**
* The operation parameter descriptor for the {@code latitudeOf1stPoint}
* parameter value. Valid values range is from -90 to 90°. Default value is 0.
*/
public static final ParameterDescriptor LAT_OF_1ST_POINT = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.ESRI, "Latitude_Of_1st_Point")
},
0, -90, 90, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@code longitudeOf1stPoint}
* parameter value. Valid values range is from -180 to 180°. Default value is 0.
*/
public static final ParameterDescriptor LONG_OF_1ST_POINT = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.ESRI, "Longitude_Of_1st_Point")
},
0, -180, 180, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@code latitudeOf2ndPoint}
* parameter value. Valid values range is from -90 to 90°. Default value is 0.
*/
public static final ParameterDescriptor LAT_OF_2ND_POINT = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.ESRI, "Latitude_Of_2nd_Point")
},
0, -90, 90, NonSI.DEGREE_ANGLE);
/**
* The operation parameter descriptor for the {@code longitudeOf2ndPoint}
* parameter value. Valid values range is from -180 to 180°. Default value is 0.
*/
public static final ParameterDescriptor LONG_OF_2ND_POINT = createDescriptor(
new NamedIdentifier[] {
new NamedIdentifier(Citations.ESRI, "Longitude_Of_2nd_Point")
},
0, -180, 180, NonSI.DEGREE_ANGLE);
/**
* The parameters group.
*/
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
new NamedIdentifier(Citations.ESRI, "Hotine_Oblique_Mercator_Two_Point_Center"),
new NamedIdentifier(Citations.GEOTOOLS, NAME)
}, new ParameterDescriptor[] {
SEMI_MAJOR, SEMI_MINOR,
LAT_OF_1ST_POINT, LONG_OF_1ST_POINT,
LAT_OF_2ND_POINT, LONG_OF_2ND_POINT,
LATITUDE_OF_CENTRE, SCALE_FACTOR,
FALSE_EASTING, FALSE_NORTHING
});
/**
* Constructs a new provider.
*/
public Provider_TwoPoint() {
super(PARAMETERS);
}
/**
* Constructs a new provider.
*
* @param params A description of parameters.
*/
protected Provider_TwoPoint(final ParameterDescriptorGroup params) {
super(params);
}
/**
* 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.
*/
@Override
protected MathTransform createMathTransform(final ParameterValueGroup parameters)
throws ParameterNotFoundException
{
final Collection<GeneralParameterDescriptor> descriptors = PARAMETERS.descriptors();
return new ObliqueMercator(parameters, descriptors, true, false);
}
}
}