Return-Path: X-Original-To: apmail-sis-commits-archive@www.apache.org Delivered-To: apmail-sis-commits-archive@www.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 3474A182A0 for ; Mon, 19 Oct 2015 18:18:08 +0000 (UTC) Received: (qmail 53243 invoked by uid 500); 19 Oct 2015 18:17:58 -0000 Delivered-To: apmail-sis-commits-archive@sis.apache.org Received: (qmail 53151 invoked by uid 500); 19 Oct 2015 18:17:56 -0000 Mailing-List: contact commits-help@sis.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: sis-dev@sis.apache.org Delivered-To: mailing list commits@sis.apache.org Received: (qmail 53142 invoked by uid 99); 19 Oct 2015 18:17:56 -0000 Received: from Unknown (HELO spamd3-us-west.apache.org) (209.188.14.142) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 19 Oct 2015 18:17:56 +0000 Received: from localhost (localhost [127.0.0.1]) by spamd3-us-west.apache.org (ASF Mail Server at spamd3-us-west.apache.org) with ESMTP id 67481180E26 for ; Mon, 19 Oct 2015 18:17:56 +0000 (UTC) X-Virus-Scanned: Debian amavisd-new at spamd3-us-west.apache.org X-Spam-Flag: NO X-Spam-Score: 1.79 X-Spam-Level: * X-Spam-Status: No, score=1.79 tagged_above=-999 required=6.31 tests=[KAM_ASCII_DIVIDERS=0.8, KAM_LAZY_DOMAIN_SECURITY=1, T_RP_MATCHES_RCVD=-0.01] autolearn=disabled Received: from mx1-eu-west.apache.org ([10.40.0.8]) by localhost (spamd3-us-west.apache.org [10.40.0.10]) (amavisd-new, port 10024) with ESMTP id T9RZ04mqzVqJ for ; Mon, 19 Oct 2015 18:17:40 +0000 (UTC) Received: from mailrelay1-us-west.apache.org (mailrelay1-us-west.apache.org [209.188.14.139]) by mx1-eu-west.apache.org (ASF Mail Server at mx1-eu-west.apache.org) with ESMTP id B7F7120751 for ; Mon, 19 Oct 2015 18:17:39 +0000 (UTC) Received: from svn01-us-west.apache.org (svn.apache.org [10.41.0.6]) by mailrelay1-us-west.apache.org (ASF Mail Server at mailrelay1-us-west.apache.org) with ESMTP id C51E7E0718 for ; Mon, 19 Oct 2015 18:17:38 +0000 (UTC) Received: from svn01-us-west.apache.org (localhost [127.0.0.1]) by svn01-us-west.apache.org (ASF Mail Server at svn01-us-west.apache.org) with ESMTP id 2ECFB3A081D for ; Mon, 19 Oct 2015 18:17:38 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 8bit Subject: svn commit: r1709462 - in /sis/branches/JDK7: ./ core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ Date: Mon, 19 Oct 2015 18:17:37 -0000 To: commits@sis.apache.org From: desruisseaux@apache.org X-Mailer: svnmailer-1.0.9 Message-Id: <20151019181738.2ECFB3A081D@svn01-us-west.apache.org> Author: desruisseaux Date: Mon Oct 19 18:17:37 2015 New Revision: 1709462 URL: http://svn.apache.org/viewvc?rev=1709462&view=rev Log: Merge projection work from JDK8 branch. Modified: sis/branches/JDK7/ (props changed) sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NormalizedProjectionTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java Propchange: sis/branches/JDK7/ ------------------------------------------------------------------------------ --- svn:mergeinfo (original) +++ svn:mergeinfo Mon Oct 19 18:17:37 2015 @@ -1,4 +1,4 @@ /sis/branches/Android:1430670-1480699 /sis/branches/JDK6:1394913-1508480 -/sis/branches/JDK8:1584960-1709036 +/sis/branches/JDK8:1584960-1709461 /sis/trunk:1394364-1508466,1519089-1519674 Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -305,10 +305,9 @@ final class Initializer { } /** - * Returns the radius of the conformal sphere at a given latitude. - * The radius of conformal sphere is computed from ρ, which is the radius of curvature in - * the meridian at latitude φ, and ν which is the radius of curvature in the prime vertical, - * as below: + * Returns the radius of the conformal sphere (assuming a semi-major axis length of 1) at a given latitude. + * The radius of conformal sphere is computed from ρ, which is the radius of curvature in the meridian at + * latitude φ, and ν which is the radius of curvature in the prime vertical, as below: * *
Rc = √(ρ⋅ν) = √(1 – ℯ²) / (1 – ℯ²sin²φ)
* @@ -321,9 +320,9 @@ final class Initializer { */ final double radiusOfConformalSphere(final double sinφ) { final DoubleDouble Rc = verbatim(1); - Rc.subtract(excentricitySquared); //-- 1 - ℯ² - Rc.sqrt(); //-- √(1 - ℯ²) - Rc.divide(rν2(sinφ)); //-- √(1 - ℯ²) / (1 - ℯ²sin²φ) + Rc.subtract(excentricitySquared); // 1 - ℯ² + Rc.sqrt(); // √(1 - ℯ²) + Rc.divide(rν2(sinφ)); // √(1 - ℯ²) / (1 - ℯ²sin²φ) return Rc.value; } Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -353,7 +353,7 @@ public class LambertConicConformal exten * comparing two {@code LambertConicConformal} projections or formatting them in debug mode. */ @Override - String[] getInternalParameterNames() { + final String[] getInternalParameterNames() { return new String[] {"n"}; } @@ -362,7 +362,7 @@ public class LambertConicConformal exten * comparing two {@code LambertConicConformal} projections or formatting them in debug mode. */ @Override - double[] getInternalParameterValues() { + final double[] getInternalParameterValues() { return new double[] {n}; } @@ -405,7 +405,7 @@ public class LambertConicConformal exten * the first non-linear one moved to the "normalize" affine transform, and the linear operations * applied after the last non-linear one moved to the "denormalize" affine transform. */ - final double θ = srcPts[srcOff ]; // θ = λ⋅n + final double θ = srcPts[srcOff ]; // θ = λ⋅n (ignoring longitude of origin) final double φ = srcPts[srcOff+1]; // Sign may be reversed final double absφ = abs(φ); final double sinθ = sin(θ); @@ -430,9 +430,9 @@ public class LambertConicConformal exten if (!derivate) { return null; } - // - // End of map projection. Now compute the derivative. - // + /* + * End of map projection. Now compute the derivative. + */ final double dρ; if (sinφ != 1) { dρ = n * dy_dφ(sinφ, cos(φ)) * ρ; @@ -444,7 +444,7 @@ public class LambertConicConformal exten } /** - * Transforms the specified (x,y) coordinates and stores the (θ,φ) result in {@code dstPts}. + * Converts the specified (x,y) coordinates and stores the (θ,φ) result in {@code dstPts}. * * @throws ProjectionException if the point can not be converted. */ Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -404,7 +404,7 @@ public class Mercator extends ConformalP } /** - * Transforms the specified (x,y) coordinates + * Converts the specified (x,y) coordinates * and stores the result in {@code dstPts} (angles in radians). * * @throws ProjectionException if the point can not be converted. Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -23,7 +23,6 @@ import org.opengis.referencing.operation import org.opengis.referencing.operation.MathTransform; import org.opengis.referencing.operation.MathTransformFactory; import org.opengis.util.FactoryException; - import org.apache.sis.parameter.Parameters; import org.apache.sis.util.Workaround; import org.apache.sis.referencing.operation.matrix.Matrix2; @@ -36,19 +35,27 @@ import static org.apache.sis.internal.re /** * Oblique Stereographic projection (EPSG code 9809). - * The formulas used below are from the EPSG guide. + * See the Stereographic projection + * on MathWorld for an overview. + * + *
Description
+ * The directions starting from the central point are true, but the areas and the lengths become + * increasingly deformed as one moves away from the center. This projection is frequently used + * for mapping polar areas, but can also be used for other limited areas centered on a point. * - *
References
- *
    - *
  • {@code libproj4} is available at - * EPSG guide.
    - * Relevant files are: {@code PJ_sterea.c}, {@code pj_gauss.c}, - * {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}
    • - *
+ *

This projection involves two steps: first a conversion of geodetic coordinates to conformal + * coordinates (i.e. latitudes and longitudes on a conformal sphere), then a spherical stereographic projection. + * For this reason this projection method is sometime known as "Double Stereographic".

+ * + *
Note: + * there is another method known as "Oblique Stereographic Alternative" or sometime just + * "Stereographic". That alternative uses a simplified conversion computing the conformal latitude + * of each point on the ellipsoid. Both methods are considered valid but produce slightly different results. + * For this reason EPSG considers them as different projection methods.
* * @author Rémi Maréchal (Geomatys) * @author Martin Desruisseaux (Geomatys) - * @since 0.6 + * @since 0.7 * @version 0.7 * @module */ @@ -59,68 +66,32 @@ public class ObliqueStereographic extend private static final long serialVersionUID = -1454098847621943639L; /** - * Conformal latitude of origin only use - * into {@link #inverseTransform(double[], int, double[], int) }. + * Conformal latitude of origin (χ₀), together with its sine and cosine. + * In the spherical case, χ₀ = φ₀ (the geodetic latitude of origin). */ - private final double χ0; + final double χ0, sinχ0, cosχ0; /** - * Value of sin(χ0) only use - * into {@link #transform(double[], int, double[], int, boolean) }. - * - * @see #χ0 + * Parameters used in the conformal sphere definition. Those parameters are used in both the + * {@linkplain #transform(double[], int, double[], int, boolean) forward} and + * {@linkplain #inverseTransform(double[], int, double[], int) inverse} projection. + * If the user-supplied ellipsoid is already a sphere, then those parameters are equal to 1. */ - private final double sinχ0; + private final double c, n; /** - * Value of cos(χ0) only use - * into {@link #transform(double[], int, double[], int, boolean) }. - * - * @see #χ0 + * Parameters used in the {@linkplain #inverseTransform(double[], int, double[], int) inverse} projection. + * More precisely g and h are used to compute intermediate parameters i + * and j, which are themselves used to compute conformal latitude and longitude. */ - private final double cosχ0; + private final double g, h; /** - * c, internaly parameter used to define conformal sphere, used - * into {@link #transform(double[], int, double[], int, boolean) } - * and {@link #inverseTransform(double[], int, double[], int) }. - */ - private final double c; - - /** - * n, internaly parameter used to define conformal sphere, used - * into {@link #transform(double[], int, double[], int, boolean) } - * and {@link #inverseTransform(double[], int, double[], int) }. - */ - private final double n; - - /** - * g, internaly parameter used to define conformal sphere coordinate conversion, - * during {@link #inverseTransform(double[], int, double[], int) }. - * More precisely g is used to compute i and j parameters and i and j, - * are used to compute only conformal longitude. - */ - private final double g; - - /** - * h, internaly parameter used to define conformal sphere coordinate conversion, - * during {@link #inverseTransform(double[], int, double[], int) }. - * More precisely h is used to compute i and j parameters and i and j, - * are used to compute only conformal longitude. - */ - private final double h; - - /** - * A convenient computing for 1 - {@link #excentricitySquared}. - */ - private final double eS1; - - /** - * Creates a Oblique Stereographic projection from the given parameters. + * Creates an Oblique Stereographic projection from the given parameters. * The {@code method} argument can be the description of one of the following: * *
    - *
  • "Oblique Stereographic".
    • + *
  • "Oblique Stereographic", also known as "Roussilhe".
    • *
* * @param method Description of the projection parameters. @@ -150,47 +121,52 @@ public class ObliqueStereographic extend */ private ObliqueStereographic(final Initializer initializer) { super(initializer); - - eS1 = 1 - excentricitySquared; - - final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN)); - + final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN)); + final double sinφ0 = sin(φ0); final double cosφ0 = cos(φ0); final double cos4_φ0 = pow(cosφ0, 4); - n = sqrt((1 + ((excentricitySquared * cos4_φ0) / eS1))); - - final double sinφ0 = sin(φ0); - final double esinφ0 = excentricity * sinφ0; - - final double s1 = (1 + sinφ0) / (1 - sinφ0); - final double s2 = (1 - esinφ0) / (1 + esinφ0); - final double w1 = pow(s1 * pow(s2, excentricity), n); - + final double ℯsinφ0 = excentricity * sinφ0; + n = sqrt(1 + ((excentricitySquared * cos4_φ0) / (1 - excentricitySquared))); /* - * Original formula : sinχ0 = ... - * To avoid confusion with χ0 conformal latitude of origin, - * renamed sinχ0 into sinχc. + * Following variables use upper-case because they are written that way in the EPSG guide. + */ + final double S1 = (1 + sinφ0) / (1 - sinφ0); + final double S2 = (1 - ℯsinφ0) / (1 + ℯsinφ0); + final double w1 = pow(S1 * pow(S2, excentricity), n); + /* + * The χ₁ variable below was named χ₀ in the EPSG guide. We use the χ₁ name in order to avoid confusion with + * the conformal latitude of origin, which is also named χ₀ in the EPSG guide. Mathematically, χ₀ and χ₁ are + * computed in the same way except that χ₁ is computed with w₁ and χ₀ is computed with w₀. + */ + final double sinχ1 = (w1 - 1) / (w1 + 1); + c = ((n + sinφ0) * (1 - sinχ1)) / + ((n - sinφ0) * (1 + sinχ1)); + /* + * Convert the geodetic latitude of origin φ₀ to the conformal latitude of origin χ₀. */ - final double sinχc = (w1 - 1) / (w1 + 1); - c = (n + sinφ0) * (1 - sinχc) / ((n - sinφ0) * (1 + sinχc)); - - //-- for invert formula final double w2 = c * w1; - χ0 = asin((w2 - 1) / (w2 + 1)); - - sinχ0 = sin(χ0); + sinχ0 = (w2 - 1) / (w2 + 1); + χ0 = asin(sinχ0); cosχ0 = cos(χ0); - - final double R = initializer.radiusOfConformalSphere(sinφ0); - - g = tan(PI / 4 - χ0 / 2); - h = 2 * tan(χ0) + g; - + /* + * Following variables are used only by the inverse projection. + */ + g = tan(PI/4 - χ0/2); + h = 2*tan(χ0) + g; + /* + * One of the first steps performed by the stereographic projection is to multiply the longitude by n. + * Since this is a linear operation, we can combine it with other linear operations performed by the + * normalization matrix. + */ final MatrixSIS normalize = context.getMatrix(true); final MatrixSIS denormalize = context.getMatrix(false); normalize.convertAfter(0, n, null); - - final double R2 = 2 * R; + /* + * One of the last steps performed by the stereographic projection is to multiply the easting and northing + * by 2 times the radius of the conformal sphere. Since this is a linear operation, we combine it with other + * linear operations performed by the denormalization matrix. + */ + final double R2 = 2 * initializer.radiusOfConformalSphere(sinφ0); denormalize.convertBefore(0, R2, null); denormalize.convertBefore(1, R2, null); } @@ -207,7 +183,28 @@ public class ObliqueStereographic extend n = other.n; g = other.g; h = other.h; - eS1 = other.eS1; + } + + /** + * Returns the names of additional internal parameters which need to be taken in account when + * comparing two {@code ObliqueStereographic} projections or formatting them in debug mode. + * + *

We could report any of the internal parameters. But since they are all derived from φ₀ and + * the {@linkplain #excentricity excentricity} and since the excentricity is already reported by + * the super-class, we report only χ₀ is a representative of the internal parameters.

+ */ + @Override + final String[] getInternalParameterNames() { + return new String[] {"χ₀"}; + } + + /** + * Returns the values of additional internal parameters which need to be taken in account when + * comparing two {@code ObliqueStereographic} projections or formatting them in debug mode. + */ + @Override + final double[] getInternalParameterValues() { + return new double[] {χ0}; } /** @@ -226,7 +223,7 @@ public class ObliqueStereographic extend public MathTransform createMapProjection(final MathTransformFactory factory) throws FactoryException { ObliqueStereographic kernel = this; if (excentricity == 0) { -// kernel = new Spherical(this); // not implemented yet + kernel = new Spherical(this); } return context.completeTransform(factory, kernel); } @@ -240,135 +237,196 @@ public class ObliqueStereographic extend * @throws ProjectionException if the coordinate can not be converted. */ @Override - public Matrix transform(double[] srcPts, int srcOff, double[] dstPts, int dstOff, boolean derivate) throws ProjectionException { - final double φ = srcPts[srcOff + 1]; - final double λ = srcPts[srcOff]; - - final double sinφ = sin(φ); - final double esinφ = excentricity * sinφ; - final double v1_sinφ = 1 - sinφ; - final double v1esinφ = 1 + esinφ; - final double Sa = (1 + sinφ) / v1_sinφ; - final double Sb = (1 - esinφ) / v1esinφ; - final double sbpowex = pow(Sb, excentricity); - final double sasbpowex = Sa * sbpowex; - final double w = c * pow(sasbpowex, n); - final double w1 = w + 1; - final double w1_w1 = (w - 1) / w1; + public Matrix transform(final double[] srcPts, final int srcOff, + final double[] dstPts, final int dstOff, + final boolean derivate) throws ProjectionException + { + final double Λ = srcPts[srcOff ]; // Λ = λ⋅n (see below), ignoring longitude of origin. + final double φ = srcPts[srcOff+1]; + final double sinφ = sin(φ); + final double ℯsinφ = excentricity * sinφ; + final double Sa = (1 + sinφ) / (1 - sinφ); + final double Sb = (1 - ℯsinφ) / (1 + ℯsinφ); + final double w = c * pow(Sa * pow(Sb, excentricity), n); /* - * Sometimes to compute projection coordinates values, computing pass by a - * "conformal sphere" to approximate as better, destination projection coordinates. + * Convert the geodetic coordinates (φ,λ) to conformal coordinates (χ,Λ) before to apply the + * actual stereographic projection. The geodetic and conformal coordinates will be the same + * if the ellipsoid is already a sphere. */ - //-- latitude coordinate into conformal sphere space. - final double χ = asin(w1_w1); + final double χ = asin((w - 1) / (w + 1)); final double cosχ = cos(χ); final double sinχ = sin(χ); /* - * Longitude coordinate into conformal sphere space is Λ = n(λ–ΛO)+ ΛO. - * But in our case, all of this linears computing are delegate into - * normalize matrix. See contructor for more precisions. - * We work directly with λ. + * The conformal longitude is Λ = n⋅(λ - λ₀) + Λ₀ where λ is the geodetic longitude. + * But in Apache SIS implementation, the multiplication by n has been merged in the + * constructor with other linear operations performed by the "normalization" matrix. + * Consequently the value obtained at srcPts[srcOff] is already Λ - Λ₀, not λ - λ₀. */ - final double cosλ = cos(λ); - final double sinλ = sin(λ); + final double cosΛ = cos(Λ); + final double sinΛ = sin(Λ); /* - * Now transform conformal sphere coordinates - * into projection destination space + * Now apply the stereographic projection on the conformal sphere using the (χ,Λ) coordinates. + * Note that the formulas below are the same than the formulas in the Spherical inner class. + * The only significant difference is that the spherical case does not contain all the above + * code which converted (φ,λ) into (χ,Λ). */ final double sinχsinχ0 = sinχ * sinχ0; final double cosχcosχ0 = cosχ * cosχ0; - final double cosχsinλ = cosχ * sinλ; - - final double B = 1 + sinχsinχ0 + cosχcosχ0 * cosλ; - - final double y = (sinχ * cosχ0 - cosχ * sinχ0 * cosλ) / B; - final double x = cosχsinλ / B; - + final double B = 1 + sinχsinχ0 + cosχcosχ0 * cosΛ; if (dstPts != null) { - dstPts[dstOff ] = x; - dstPts[dstOff+1] = y; + dstPts[dstOff ] = cosχ * sinΛ / B; // Easting (x) + dstPts[dstOff+1] = (sinχ * cosχ0 - cosχ * sinχ0 * cosΛ) / B; // Northing (y) } - if (!derivate) { return null; } + /* + * Now compute the derivative, if the user asked for it. + * Notes: + * + * ∂Sa/∂λ = 0 + * ∂Sb/∂λ = 0 + * ∂w/∂λ = 0 + * ∂χ/∂λ = 0 + * ∂Sa/∂φ = 2⋅cosφ / (1 - sinφ)² + * ∂Sb/∂φ = -2⋅ℯ⋅cosφ / (1 - ℯ⋅sinφ)² + * ∂w/∂φ = 2⋅n⋅w⋅(1/cosφ - ℯ²⋅cosφ/(1 - ℯ²⋅sin²φ)); + */ + final double cosφ = cos(φ); + final double dχ_dφ = (1/cosφ - cosφ*excentricitySquared/(1 - ℯsinφ*ℯsinφ)) * 2*n*sqrt(w) / (w + 1); - final double cosφ = cos(φ); - final double B2 = B * B; - - //-- derivative code - //-- dSa_dλ = 0; - final double dSa_dφ = 2 * cosφ / (v1_sinφ * v1_sinφ); - - //-- dSb_dλ = 0; - final double dSb_dφ = - 2 * excentricity * cosφ / (v1esinφ * v1esinφ); - - //-- dsasbpowex_dλ = 0; - final double dsasbpowex_dφ = dSa_dφ * sbpowex + Sa * excentricity * dSb_dφ * pow(Sb, excentricity - 1); - - //-- dw_dλ = 0; - final double dw_dφ = c * n * dsasbpowex_dφ * pow(sasbpowex, n - 1); - - //-- dχ_dλ = 0; - final double dχ_dφ = dw_dφ / (w1 * sqrt(w)); - - final double addsinχsinχ0 = sinχ + sinχ0; - - //-- Jacobian coefficients - final double dx_dλ = cosχ * (cosλ * (1 + sinχsinχ0) + cosχcosχ0) / B2; - - final double dx_dφ = - dχ_dφ * sinλ * addsinχsinχ0 / B2; - - final double dy_dλ = cosχsinλ * addsinχsinχ0 / B2; - - final double dy_dφ = dχ_dφ * (cosχcosχ0 + cosλ * (sinχsinχ0 + 1)) / B2; - - return new Matrix2(dx_dλ, dx_dφ, - dy_dλ, dy_dφ); + final double B2 = B * B; + final double d = (cosχcosχ0 + cosΛ * (sinχsinχ0 + 1)) / B2; // Matrix diagonal + final double t = sinΛ * (sinχ + sinχ0) / B2; // Matrix anti-diagonal + /* ┌ ┐ + * │ ∂x/∂λ, ∂x/∂φ │ + * Jacobian = │ │ + * (Proj(λ,φ)) │ ∂y/∂λ, ∂y/∂φ │ + * └ ┘ + */ + return new Matrix2(d*cosχ, -t*dχ_dφ, + t*cosχ, d*dχ_dφ); } /** - * Transforms the specified (x, y) coordinates and stores the result in {@code dstPts} (angles in radians). + * Converts the specified (x,y) coordinates and stores the result in {@code dstPts} + * (angles in radians). * * @throws ProjectionException if the point can not be converted. */ @Override protected void inverseTransform(double[] srcPts, int srcOff, double[] dstPts, int dstOff) throws ProjectionException { - final double x = srcPts[srcOff]; - final double y = srcPts[srcOff + 1]; - + final double x = srcPts[srcOff ]; + final double y = srcPts[srcOff+1]; final double i = atan(x / (h + y)); final double j = atan(x / (g - y)) - i; /* - * Longitude coordinate into conformal sphere space is Λ = j + 2 * i - * Where λ = Λ + Λ0, but Λ0 is added into normalize matrix which regroup all linears operations. - * Also in our particularity case Geodetic longitude λ is the same. + * The conformal longitude is Λ = j + 2i + Λ₀. In the particular case of stereographic projection, + * the geodetic longitude λ is equals to Λ. Furthermore in Apache SIS implementation, Λ₀ is added by + * the denormalization matrix and shall not be handled here. The only remaining part is λ = j + 2i. */ - final double λ = j + 2*i; - - //-- latitude coordinate into conformal sphere space. - final double χ = χ0 + 2*atan((y - x*tan(j/2))); - final double sinχ = sin(χ); - - final double ψ = log((1 + sinχ) / (c * (1 - sinχ))) / (2 * n); - - double φi_1 = 2*atan(exp(ψ)) - PI/2; - + final double λ = j + 2*i; + /* + * Calculation of geodetic latitude φ involves first the calculation of conformal latitude χ, + * then calculation of isometric latitude ψ, and finally calculation of φ by an iterative method. + */ + final double sinχ = sin(χ0 + 2*atan(y - x*tan(j/2))); + final double ψ = log((1 + sinχ) / ((1 - sinχ)*c)) / (2*n); + double φ = 2*atan(exp(ψ)) - PI/2; // First approximation for (int it = 0; it < MAXIMUM_ITERATIONS; it++) { - final double sinφi_1 = sin(φi_1); - final double esinφi_1 = excentricity*sinφi_1; + final double ℯsinφ = excentricity * sin(φ); + final double ψi = log(tan(φ/2 + PI/4) * pow((1 - ℯsinφ) / (1 + ℯsinφ), excentricity/2)); + final double Δφ = (ψ - ψi) * cos(φ) * (1 - ℯsinφ*ℯsinφ) / (1 - excentricitySquared); + φ += Δφ; + if (abs(Δφ) <= ITERATION_TOLERANCE) { + dstPts[dstOff ] = λ; + dstPts[dstOff+1] = φ; + return; + } + } + throw new ProjectionException(Errors.Keys.NoConvergence); + } + + /** + * Provides the transform equations for the spherical case of the Oblique Stereographic projection. + * This implementation can be used when {@link #excentricity} = 0. + * + * @author Rémi Maréchal (Geomatys) + * @author Martin Desruisseaux (Geomatys) + * @since 0.7 + * @version 0.7 + * @module + */ + static final class Spherical extends ObliqueStereographic { + /** + * For cross-version compatibility. + */ + private static final long serialVersionUID = -1454098847621943639L; - double ψi_1 = log(tan(φi_1/2 + PI/4) * pow((1 - esinφi_1) / (1 + esinφi_1), excentricity / 2)); + /** + * Constructs a new map projection from the supplied parameters. + * + * @param parameters The parameters of the projection to be created. + */ + protected Spherical(ObliqueStereographic other) { + super(other); + } - final double φi = φi_1 - (ψi_1 - ψ) * cos(φi_1) * (1 - esinφi_1 * esinφi_1) / eS1; + /** + * {@inheritDoc} + */ + @Override + public Matrix transform(double[] srcPts, int srcOff, double[] dstPts, int dstOff, boolean derivate) throws ProjectionException { + final double λ = srcPts[srcOff ]; + final double φ = srcPts[srcOff+1]; + /* + * In the spherical case, φ = χ and λ = Λ. + */ + final double sinφ = sin(φ); + final double cosφ = cos(φ); + final double sinλ = sin(λ); + final double cosλ = cos(λ); + final double sinφsinφ0 = sinφ * sinχ0; + final double cosφcosφ0 = cosφ * cosχ0; + final double cosφsinλ = cosφ * sinλ; + final double B = 1 + sinφsinφ0 + cosφcosφ0 * cosλ; + if (dstPts != null) { + dstPts[dstOff ] = cosφsinλ / B; + dstPts[dstOff+1] = (sinφ*cosχ0 - cosφ*sinχ0 *cosλ) / B; + } + if (!derivate) { + return null; + } + final double B2 = B * B; + final double d = (cosφcosφ0 + cosλ * (sinφsinφ0 + 1)) / B2; + final double t = sinλ * (sinφ + sinχ0) / B2; + return new Matrix2(d*cosφ, -t, + t*cosφ, d); + } - if (abs(φi - φi_1) <= ITERATION_TOLERANCE) { - dstPts[dstOff] = λ; - dstPts[dstOff + 1] = φi; - return; + /** + * {@inheritDoc} + */ + @Override + protected void inverseTransform(double[] srcPts, int srcOff, double[] dstPts, int dstOff) throws ProjectionException { + final double x = srcPts[srcOff ]; + final double y = srcPts[srcOff+1]; + final double ρ = hypot(x, y); + final double λ, φ; + if (abs(ρ) < ANGULAR_TOLERANCE) { + φ = χ0; + λ = 0.0; + } else { + final double c = 2*atan(ρ); + final double cosc = cos(c); + final double sinc = sin(c); + final double ct = ρ * cosχ0*cosc - y*sinχ0*sinc; + final double t = x * sinc; + φ = asin(cosc*sinχ0 + y*sinc*cosχ0 / ρ); + λ = atan2(t, ct); } - φi_1 = φi; + dstPts[dstOff] = λ; + dstPts[dstOff+1] = φ; } - throw new ProjectionException(Errors.Keys.NoConvergence); } } Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -41,6 +41,13 @@ import static org.apache.sis.internal.ut /** * Polar Stereographic projection (EPSG codes 9810, 9829, 9830). + * This is a special case of {@link ObliqueStereographic} when the projection origin is at a pole. + * + *

EPSG defines three variants for this projection, A, B and C, + * which differ by the way the parameters are specified. The "Polar Stereographic (variant B)" + * projection includes a "Latitude of standard parallel" parameter where is effective the scale factor + * (normally 1). The "Polar Stereographic (variant A)" forces its "Latitude of natural origin" + * parameter to ±90°, depending on the hemisphere.

* * @author Gerald Evenden (USGS) * @author André Gosselin (MPO) @@ -51,7 +58,6 @@ import static org.apache.sis.internal.ut * @version 0.6 * @module * - * @see EquatorialStereographic * @see ObliqueStereographic */ public class PolarStereographic extends ConformalProjection { @@ -316,7 +322,7 @@ public class PolarStereographic extends } /** - * Transforms the specified (x,y) coordinates and stores the result in {@code dstPts} (angles in radians). + * Converts the specified (x,y) coordinates and stores the result in {@code dstPts} (angles in radians). * * @throws ProjectionException if the point can not be converted. */ Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -34,6 +34,7 @@ import org.apache.sis.util.Workaround; import static java.lang.Math.*; import static org.apache.sis.math.MathFunctions.asinh; import static org.apache.sis.math.MathFunctions.atanh; +import static org.apache.sis.internal.referencing.provider.TransverseMercator.*; /** @@ -115,10 +116,10 @@ public class TransverseMercator extends xOffset = ParameterRole.FALSE_WESTING; yOffset = ParameterRole.FALSE_SOUTHING; } - roles.put(ParameterRole.CENTRAL_MERIDIAN, org.apache.sis.internal.referencing.provider.TransverseMercator.LONGITUDE_OF_ORIGIN); - roles.put(ParameterRole.SCALE_FACTOR, org.apache.sis.internal.referencing.provider.TransverseMercator.SCALE_FACTOR); - roles.put(xOffset, org.apache.sis.internal.referencing.provider.TransverseMercator.FALSE_EASTING); - roles.put(yOffset, org.apache.sis.internal.referencing.provider.TransverseMercator.FALSE_NORTHING); + roles.put(ParameterRole.CENTRAL_MERIDIAN, LONGITUDE_OF_ORIGIN); + roles.put(ParameterRole.SCALE_FACTOR, SCALE_FACTOR); + roles.put(xOffset, FALSE_EASTING); + roles.put(yOffset, FALSE_NORTHING); return new Initializer(method, parameters, roles, isSouth ? (byte) 1 : (byte) 0); } @@ -129,8 +130,7 @@ public class TransverseMercator extends @Workaround(library="JDK", version="1.7") private TransverseMercator(final Initializer initializer) { super(initializer); - final double φ0 = toRadians(initializer.getAndStore( - org.apache.sis.internal.referencing.provider.TransverseMercator.LATITUDE_OF_ORIGIN)); + final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN)); /* * Opportunistically use double-double arithmetic for computation of B since we will store * it in the denormalization matrix, and there is no sine/cosine functions involved here. @@ -383,33 +383,30 @@ public class TransverseMercator extends sinh_8η0 = sinh_4η0 * cosh_4η0; assert identityEquals(sinh_8η0, sinh(8*η0) / 8) : η0; } /* - * Assuming that (λ, φ) ↦ Proj((λ, φ)) - * where Proj is defined by: Proj((λ, φ)) : (η(λ, φ), ξ(λ, φ)). - * - * => (λ, φ) ↦ (η(λ, φ), ξ(λ, φ)). + * The projection of (λ,φ) is given by (η⋅B, ξ⋅B+M₀) — ignoring scale factors and false easting/northing. + * But the B and M₀ parameters have been merged by the constructor with other linear operations in the + * "denormalization" matrix. Consequently we only need to compute (η,ξ) below. */ - //-- ξ(λ, φ) - final double ξ = cf8 * sin_8ξ0 * cosh_8η0 - + cf6 * sin_6ξ0 * cosh_6η0 - + cf4 * sin_4ξ0 * cosh_4η0 - + cf2 * sin_2ξ0 * cosh_2η0 - + ξ0; - - //-- η(λ, φ) - final double η = cf8 * cos_8ξ0 * sinh_8η0 - + cf6 * cos_6ξ0 * sinh_6η0 - + cf4 * cos_4ξ0 * sinh_4η0 - + cf2 * cos_2ξ0 * sinh_2η0 - + η0; - if (dstPts != null) { - dstPts[dstOff ] = η; - dstPts[dstOff+1] = ξ; + // η(λ,φ) + dstPts[dstOff ] = cf8 * cos_8ξ0 * sinh_8η0 + + cf6 * cos_6ξ0 * sinh_6η0 + + cf4 * cos_4ξ0 * sinh_4η0 + + cf2 * cos_2ξ0 * sinh_2η0 + + η0; + // ξ(λ,φ) + dstPts[dstOff+1] = cf8 * sin_8ξ0 * cosh_8η0 + + cf6 * sin_6ξ0 * cosh_6η0 + + cf4 * sin_4ξ0 * cosh_4η0 + + cf2 * sin_2ξ0 * cosh_2η0 + + ξ0; } if (!derivate) { return null; } - + /* + * Now compute the derivative, if the user asked for it. + */ final double cosλ = cos(λ); //-- λ final double cosφ = cos(φ); //-- φ final double cosh2Q = coshQ * coshQ; //-- Q @@ -428,14 +425,14 @@ public class TransverseMercator extends final double dξ0_dλ = sinhQ * sinhη0 * cosλ / (cosh2Q_sin2λ * sqrt1_thQchη0); final double dξ0_dφ = (dQ_dφ * coshη0 / cosh2Q + dη0_dφ * sinhη0 * tanhQ) / sqrt1_thQchη0; /* - * Assuming that Jac(Proj((λ, φ))) is the Jacobian matrix of Proj((λ, φ)) function. + * Jac(Proj(λ,φ)) is the Jacobian matrix of Proj(λ,φ) function. + * So the derivative of Proj(λ,φ) is defined by: * - * So the derivative of Proj((λ, φ)) is defined by: - * ┌ ┐ - * │ dη(λ, φ) / dλ, dη(λ, φ) / dφ │ - * Jac = │ │ - * (Proj(λ, φ)) │ dξ(λ, φ) / dλ, dξ(λ, φ) / dφ │ - * └ ┘ + * ┌ ┐ + * │ ∂η(λ,φ)/∂λ, ∂η(λ,φ)/∂φ │ + * Jac = │ │ + * (Proj(λ,φ)) │ ∂ξ(λ,φ)/∂λ, ∂ξ(λ,φ)/∂φ │ + * └ ┘ */ //-- dξ(λ, φ) / dλ final double dξ_dλ = dξ0_dλ Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -17,19 +17,17 @@ package org.apache.sis.referencing.operation.projection; import java.util.EnumMap; - import org.opengis.parameter.ParameterDescriptor; import org.opengis.parameter.ParameterValueGroup; import org.opengis.referencing.operation.OperationMethod; - import org.apache.sis.internal.referencing.Formulas; import org.apache.sis.internal.referencing.provider.ObliqueStereographic; import org.apache.sis.parameter.Parameters; +import org.opengis.test.TestCase; import org.junit.Test; -import static org.apache.sis.internal.referencing.provider.ObliqueStereographic.*; +import static java.lang.StrictMath.*; import static org.opengis.test.Assert.*; -import org.opengis.test.TestCase; /** @@ -41,40 +39,52 @@ import org.opengis.test.TestCase; * @version 0.7 * @module */ -public strictfp class InitializerTest extends TestCase{ +public final strictfp class InitializerTest extends TestCase{ /** * Tests the {@link Initializer#radiusOfConformalSphere(double)} method. - * This test compute the Radius of conformal Sphere using the values given by the - * EPSG guide for the Stereographic projection. - * - * @see Initializer#radiusOfConformalSphere(double) + * This test computes the Radius of Conformal Sphere using the values given + * by the EPSG guide for + * the Amersfoort / RD New projection (a Stereographic one). */ @Test public void testRadiusOfConformalSphere() { - final EnumMap> roles = - new EnumMap<>(NormalizedProjection.ParameterRole.class); - roles.put(NormalizedProjection.ParameterRole.CENTRAL_MERIDIAN, LONGITUDE_OF_ORIGIN); - roles.put(NormalizedProjection.ParameterRole.SCALE_FACTOR, SCALE_FACTOR); - roles.put(NormalizedProjection.ParameterRole.FALSE_EASTING, FALSE_EASTING); - roles.put(NormalizedProjection.ParameterRole.FALSE_NORTHING, FALSE_NORTHING); - final OperationMethod op = new ObliqueStereographic(); final ParameterValueGroup p = op.getParameters().createValue(); - - //-- Implicit parameters (OGC names). + /* + * Following parameters are not given explicitely by EPSG definitions since they are + * usually inferred from the datum. However in the particular case of this test, we + * need to provide them. The names used below are either OGC names or SIS extensions. + */ p.parameter("semi_major").setValue(6377397.155); p.parameter("inverse_flattening").setValue(299.15281); - - //-- Explicit parameters from EPSG registry + /* + * Following parameters are reproduced verbatim from EPSG registry and EPSG guide. + */ p.parameter("Latitude of natural origin").setValue(52.156160556); p.parameter("Longitude of natural origin").setValue(5.387638889); p.parameter("Scale factor at natural origin").setValue(0.9999079); p.parameter("False easting").setValue(155000.00); p.parameter("False northing").setValue(463000.00); - + /* + * The following lines are a typical way to create an Initializer instance. + * The EnumMap tells to the Initializer constructor which parameters to look for. + * We construct this map here for testing purpose, but users normally do not have + * to do that since this map is provided by the ObliqueStereographic class itself. + */ + final EnumMap> roles = + new EnumMap<>(NormalizedProjection.ParameterRole.class); + roles.put(NormalizedProjection.ParameterRole.CENTRAL_MERIDIAN, ObliqueStereographic.LONGITUDE_OF_ORIGIN); + roles.put(NormalizedProjection.ParameterRole.SCALE_FACTOR, ObliqueStereographic.SCALE_FACTOR); + roles.put(NormalizedProjection.ParameterRole.FALSE_EASTING, ObliqueStereographic.FALSE_EASTING); + roles.put(NormalizedProjection.ParameterRole.FALSE_NORTHING, ObliqueStereographic.FALSE_NORTHING); final Initializer initializer = new Initializer(op, (Parameters) p, roles, (byte) 0); - + /* + * The following lines give an example of how Apache SIS projection constructors + * use the Initializer class. + */ + final double φ0 = toRadians(initializer.getAndStore(ObliqueStereographic.LATITUDE_OF_ORIGIN)); + assertTrue(φ0 > 0); assertEquals("Conformal Sphere Radius", 6382644.571, 6377397.155 * - initializer.radiusOfConformalSphere(Math.sin(Math.toRadians(52.156160556))), Formulas.LINEAR_TOLERANCE); + initializer.radiusOfConformalSphere(sin(φ0)), Formulas.LINEAR_TOLERANCE); } } Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -166,10 +166,10 @@ public final strictfp class LambertConic @Test @DependsOnMethod("testSpecialLatitudes") public void testDerivative() throws TransformException { - if (transform == null) { // May have been initialized by 'testSphericalCase'. - createNormalizedProjection(true, 40); // Elliptical case + if (transform == null) { // May have been initialized by 'testSphericalCase'. + createNormalizedProjection(true, 40); // Elliptical case } - final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. derivativeDeltas = new double[] {delta, delta}; tolerance = 1E-9; verifyDerivative(toRadians( 0), toRadians( 0)); Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -50,6 +50,7 @@ import static org.apache.sis.referencing public final strictfp class MercatorTest extends MapProjectionTestCase { /** * Creates a new instance of {@link Mercator} for a sphere or an ellipsoid. + * The new instance is stored in the inherited {@link #transform} field. * * @param ellipse {@code false} for a sphere, or {@code true} for WGS84 ellipsoid. */ @@ -116,12 +117,12 @@ public final strictfp class MercatorTest @Test @DependsOnMethod("testSpecialLatitudes") public void testDerivative() throws TransformException { - if (transform == null) { // May have been initialized by 'testSphericalCase'. - createNormalizedProjection(true); // Elliptical case + if (transform == null) { // May have been initialized by 'testSphericalCase'. + createNormalizedProjection(true); // Elliptical case } - final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. derivativeDeltas = new double[] {delta, delta}; - tolerance = 1E-9; // More severe than Formulas.LINEAR_TOLERANCE. + tolerance = 1E-9; // More severe than Formulas.LINEAR_TOLERANCE. verifyDerivative(toRadians(15), toRadians( 30)); verifyDerivative(toRadians(10), toRadians(-60)); } Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NormalizedProjectionTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NormalizedProjectionTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NormalizedProjectionTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NormalizedProjectionTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -39,7 +39,8 @@ import static org.apache.sis.internal.jd @DependsOn({ // Following dependency is where the basic parameters (e.g. SEMI_MAJOR) are tested. // Those parameters are needed by NoOp pseudo-projection, which is used in this package. - org.apache.sis.internal.referencing.provider.MapProjectionTest.class + org.apache.sis.internal.referencing.provider.MapProjectionTest.class, + InitializerTest.class }) public final strictfp class NormalizedProjectionTest extends TransformTestCase { /** Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -16,24 +16,26 @@ */ package org.apache.sis.referencing.operation.projection; +import javax.measure.unit.SI; import org.opengis.parameter.ParameterValueGroup; -import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.Matrix; import org.opengis.referencing.operation.OperationMethod; import org.opengis.referencing.operation.TransformException; import org.opengis.util.FactoryException; - -import org.apache.sis.internal.referencing.Formulas; import org.apache.sis.parameter.Parameters; import org.apache.sis.referencing.operation.transform.ContextualParameters; - -import org.junit.Assert; +import org.apache.sis.referencing.operation.matrix.Matrix2; +import org.apache.sis.internal.referencing.Formulas; +import org.apache.sis.test.DependsOnMethod; +import org.apache.sis.test.DependsOn; import org.junit.Test; -import static java.lang.Math.sqrt; +import static java.lang.StrictMath.*; +import static org.junit.Assert.*; /** - * Tests {@link ObliqueStereographic} projection. + * Tests the {@link ObliqueStereographic} class. * * @author Rémi Marechal (Geomatys) * @author Martin Desruisseaux (Geomatys) @@ -41,123 +43,169 @@ import static java.lang.Math.sqrt; * @version 0.7 * @module */ +@DependsOn({ + InitializerTest.class, + NormalizedProjectionTest.class +}) public final strictfp class ObliqueStereographicTest extends MapProjectionTestCase { /** - * Parameter values provided by the EPSG guide - * for testing {@link ObliqueStereographic} transform conformity. + * Parameter values provided by the EPSG guide + * for testing {@link ObliqueStereographic} transform conformity. The test uses the parameters for + * the Amersfoort / RD New projection: + * + *
    + *
  • Semi-major axis length: a = 6377397.155 metres
    • + *
  • Inverse flattening: 1/f = 299.15281
    • + *
  • Latitude of natural origin: φ₀ = 52°09'22.178"N
    • + *
  • Longitude of natural origin: λ₀ = 5°23'15.500"E
    • + *
  • Scale factor at natural origin: k₀ = 0.9999079
    • + *
  • False easting: FE = 155000.00 metres
    • + *
  • False northing: FN = 463000.00 metres
    • + *
+ * + * Other parameters (n, R, g, h) are computed from the above. + * Those parameters fall in three groups: + * + *
    + *
  • Parameters used in linear operations (additions or multiplications) performed before the + * non-linear part (the "kernel") of the map projection. Those parameters are λ₀ and n + * and their values are stored in the normalization matrix given by + * {@linkplain ContextualParameters#getMatrix(boolean) ContextualParameters.getMatrix}(true).
    • + * + *
  • Parameters used in linear operations (additions or multiplications) performed after the + * non-linear part (the "kernel") of the map projection. Those parameters are k₀, R, + * FE and FN and their values are stored in the denormalization matrix given by + * {@linkplain ContextualParameters#getMatrix(boolean) ContextualParameters.getMatrix}(false).
    • * - * @see ContextualParameters#getMatrix(boolean) where boolean is true for value n. - * @see ContextualParameters#getMatrix(boolean) where boolean is false for k0, a, FE, FN and R. - * @see Initializer#radiusOfConformalSphere(double) for value R + *
  • Other parameters are either used in the non-linear "kernel" of the map projection or used for computing the + * above-cited parameters.
    • + *
* + *

Note 1: value of R is computed by {@link Initializer#radiusOfConformalSphere(double)}.

+ * + *

Note 2: we do not follow the Java naming convention here (constants in upper cases) in order to use + * as much as possible the exact same symbols as in the EPSG guide.

*/ - private final static double eSQUARED = 0.08169683 * 0.08169683, // Excentricity squared. - φ0 = 0.910296727, // Latitude of natural origin (rad) - - //-- Some attributs are considered as linear and put into normalize matrix and apply before transform - n = sqrt(1 + (eSQUARED * Math.pow(Math.cos(φ0), 4)) / (1 - eSQUARED)), - - //-- Some attributs are considered as linear and put into denormalize matrix and apply just after - k0 = 0.9999079, - a = 6377397.155, - FE = 155000.00, - FN = 463000.00, - R = 6382644.571 / a; + private static final double φ0 = 0.910296727, // Latitude of natural origin (rad) + /* Before kernel */ λ0 = 0.094032038, // Longitude of natural origin (rad) + /* After kernel */ R = 6382644.571, // Radius of conformal sphere (m) + a = 6377397.155, // Semi-major axis length (m) + ivf = 299.15281, // Inverse flattening factor + e = 0.08169683, // Excentricity + /* Before kernel */ n = 1.000475857, // Coefficient computed from excentricity and φ₀. + /* After kernel */ k0 = 0.9999079, // Scale factor + /* After kernel */ FE = 155000.00, // False Easting (m) + /* After kernel */ FN = 463000.00; // False Northing (m) /** - * Tested {@link MathTransform} projection. + * Compares the n value given in the EPSG guide with the value computed from the formula. */ - private final NormalizedProjection ObliqueStereographic; + @Test + public void testN() { + assertEquals(n, sqrt(1 + (e*e * pow(cos(φ0), 4)) / (1 - e*e)), 0.5E-9); + } /** - * Buid tested {@link ObliqueStereographic} {@link MathTransform}. + * Creates a new instance of {@link ObliqueStereographic} for a sphere or an ellipsoid. + * The new instance is stored in the inherited {@link #transform} field. + * + * @param ellipse {@code false} for the spherical case, or {@code true} for the ellipsoidal case. */ - public ObliqueStereographicTest() { + private void createNormalizedProjection(final boolean ellipse) { final OperationMethod op = new org.apache.sis.internal.referencing.provider.ObliqueStereographic(); - final ParameterValueGroup p = op.getParameters().createValue(); + /* + * Following parameters are not given explicitely by EPSG definitions since they are + * usually inferred from the datum. However in the particular case of this test, we + * need to provide them. The names used below are either OGC names or SIS extensions. + */ + if (!ellipse) { + p.parameter("semi_major").setValue(R); + p.parameter("semi_minor").setValue(R); + } else { + p.parameter("semi_major").setValue(a); + p.parameter("inverse_flattening").setValue(ivf); + } + /* + * Following parameters are reproduced verbatim from EPSG registry and EPSG guide. + */ + p.parameter("Latitude of natural origin") .setValue(φ0, SI.RADIAN); + p.parameter("Longitude of natural origin") .setValue(λ0, SI.RADIAN); + p.parameter("Scale factor at natural origin").setValue(k0); + p.parameter("False easting") .setValue(FE, SI.METRE); + p.parameter("False northing") .setValue(FN, SI.METRE); - //-- implicit names from OGC. - p.parameter("semi_major").setValue(6377397.155); - p.parameter("inverse_flattening").setValue(299.15281); - - //-- Name parameters from Epsg registry - p.parameter("Latitude of natural origin").setValue(52.156160556); - p.parameter("Longitude of natural origin").setValue(5.387638889); - p.parameter("Scale factor at natural origin").setValue(0.9999079); - p.parameter("False easting").setValue(155000.00); - p.parameter("False northing").setValue(463000.00); - - ObliqueStereographic = new ObliqueStereographic(op, (Parameters) p); + transform = new ObliqueStereographic(op, (Parameters) p); } /** - * {@link MathTransform#transform(org.opengis.geometry.DirectPosition, org.opengis.geometry.DirectPosition) } - * test with expected values from - * EPSG guide + * The point given in the EPSG guide for testing the map projection. + * (φt, λt) is the source geographic coordinate in degrees and + * (xt, yt) is the target projected coordinate in metres. + */ + private static final double φt = 53, // Latitude in degrees + λt = 6, // Longitude in degrees + Et = 196105.283, // Easting in metres + Nt = 557057.739; // Northing in metres + + /** + * Tests {@link ObliqueStereographic#transform(double[], int, double[], int, boolean)} + * with the values given by the EPSG guide. * - * @throws FactoryException if an error occurred while creating the map projection. - * @throws TransformException if an error occurred while projecting a coordinate. + * @throws TransformException if an error occurred while projecting the coordinate. */ @Test - public void testEPSGTransform() throws FactoryException, TransformException { - - final double[] srcPts = new double[]{6, 53}; //-- deg - srcPts[0] = Math.toRadians(srcPts[0] - 5.387638889) ; - srcPts[1] = Math.toRadians(srcPts[1]); - + public void testTransform() throws TransformException { + final double[] srcPts = new double[] {λt, φt}; // in degrees final double[] dstPts = new double[2]; - srcPts[0] = srcPts[0] * n; + // Linear operations (normalization) applied before NormalizedTransform. + srcPts[0] = toRadians(srcPts[0]) - λ0; + srcPts[1] = toRadians(srcPts[1]); + srcPts[0] *= n; + + // The non-linear part of map projection (the "kernel"). + createNormalizedProjection(true); + transform.transform(srcPts, 0, dstPts, 0, 1); + + // Linear operations (denormalization) applied after NormalizedTransform. + dstPts[0] *= (k0 * 2*R); + dstPts[1] *= (k0 * 2*R); + dstPts[0] += FE; + dstPts[1] += FN; - ObliqueStereographic.transform(srcPts, 0, dstPts, 0, 1); - - final double destE = dstPts[0] * k0 * a * 2 * R + FE; - final double destN = dstPts[1] * k0 * a * 2 * R + FN; - - Assert.assertEquals("destination East coordinate", 196105.283, destE, Formulas.LINEAR_TOLERANCE); - Assert.assertEquals("destination North coordinate", 557057.739, destN, Formulas.LINEAR_TOLERANCE); + assertEquals("Easting", Et, dstPts[0], Formulas.LINEAR_TOLERANCE); + assertEquals("Northing", Nt, dstPts[1], Formulas.LINEAR_TOLERANCE); } - - /** - * Test method {@link ObliqueStereographic#inverseTransform(double[], int, double[], int)} - * test with expected values from - * EPSG guide + /** + * Tests {@link ObliqueStereographic#inverseTransform(double[], int, double[], int)} + * with the values given by the EPSG guide. * - * @throws org.apache.sis.referencing.operation.projection.ProjectionException + * @throws TransformException if an error occurred while projecting the coordinate. */ @Test - public void testEPSGinvertTransform() throws ProjectionException { - - double srcEast = 196105.28; - double srcNorth = 557057.74; - - srcEast -= FE; - srcNorth -= FN; - srcEast /= k0; - srcNorth /= k0; - srcEast /= a; - srcNorth /= a; - srcEast /= (2 * R); - srcNorth /= (2 * R); - - //-- tcheck transform - final double[] srcPts = new double[]{srcEast, srcNorth}; //-- meter - + public void testInverseTransform() throws TransformException { + final double[] srcPts = new double[] {Et, Nt}; // in metres final double[] dstPts = new double[2]; - ObliqueStereographic.inverseTransform(srcPts, 0, dstPts, 0); - final double λO = 0.094032038; + // Linear operations (normalization) applied before NormalizedTransform. + srcPts[0] -= FE; + srcPts[1] -= FN; + srcPts[0] /= (k0 * 2*R); + srcPts[1] /= (k0 * 2*R); + + // The non-linear part of map projection (the "kernel"). + createNormalizedProjection(true); + ((NormalizedProjection) transform).inverseTransform(srcPts, 0, dstPts, 0); + + // Linear operations (denormalization) applied after NormalizedTransform. + dstPts[0] /= n; + dstPts[0] = toDegrees(dstPts[0] + λ0); + dstPts[1] = toDegrees(dstPts[1]); - double destλ = dstPts[0] / n + λO; - double destφ = dstPts[1]; - - destλ = Math.toDegrees(destλ); - destφ = Math.toDegrees(destφ); - - Assert.assertEquals("destination East coordinate", 6, destλ, Formulas.ANGULAR_TOLERANCE); - Assert.assertEquals("destination North coordinate", 53, destφ, Formulas.ANGULAR_TOLERANCE); + assertEquals("Longitude", λt, dstPts[0], Formulas.ANGULAR_TOLERANCE); + assertEquals("Latitude", φt, dstPts[1], Formulas.ANGULAR_TOLERANCE); } /** @@ -170,7 +218,142 @@ public final strictfp class ObliqueStere * @see org.opengis.test.referencing.ParameterizedTransformTest#testObliqueStereographic() */ @Test - public void testGeoapi() throws FactoryException, TransformException { + @DependsOnMethod({"testTransform", "testInverseTransform"}) + public void testObliqueStereographic() throws FactoryException, TransformException { createGeoApiTest(new org.apache.sis.internal.referencing.provider.ObliqueStereographic()).testObliqueStereographic(); } + + /** + * Verifies the consistency of spherical formulas with the elliptical formulas. + * This test transforms the point given in the EPSG guide and takes the result + * of the elliptical implementation as a reference. + * + * @throws TransformException if an error occurred while projecting the coordinate. + */ + @Test + @DependsOnMethod("testTransform") + public void testSphericalTransform() throws TransformException { + final double[] srcPts = new double[] {λt, φt}; // in degrees + final double[] dstPts = new double[2]; + final double[] refPts = new double[2]; + + // Linear operations (normalization) applied before NormalizedTransform. + srcPts[0] = toRadians(srcPts[0]) - λ0; + srcPts[1] = toRadians(srcPts[1]); + srcPts[0] *= n; + + // The non-linear part of map projection (the "kernel"). + createNormalizedProjection(false); + transform.transform(srcPts, 0, refPts, 0, 1); + + // Linear operations (denormalization) applied after NormalizedTransform. + refPts[0] *= (k0 * 2*R); + refPts[1] *= (k0 * 2*R); + refPts[0] += FE; + refPts[1] += FN; + + // Transform the same point, now using the spherical implementation. + ObliqueStereographic spherical = (ObliqueStereographic) transform; + spherical = new ObliqueStereographic.Spherical(spherical); + spherical.transform(srcPts, 0, dstPts, 0, 1); + + // Linear operations (denormalization) applied after NormalizedTransform. + dstPts[0] *= (k0 * 2*R); + dstPts[1] *= (k0 * 2*R); + dstPts[0] += FE; + dstPts[1] += FN; + + // Use a smaller tolerance because spherical and elliptical formulas should be equivalent in this case. + assertArrayEquals("Spherical projection", refPts, dstPts, Formulas.ANGULAR_TOLERANCE / 1E6); + } + + /** + * Verifies the consistency of spherical formulas with the elliptical formulas. + * This test computes the inverse transform of the point given in the EPSG guide + * and takes the result of the elliptical implementation as a reference. + * + * @throws TransformException if an error occurred while projecting the coordinate. + */ + @Test + @DependsOnMethod("testInverseTransform") + public void testSphericalInverseTransform() throws TransformException { + final double[] srcPts = new double[] {Et, Nt}; // in metres + final double[] dstPts = new double[2]; + final double[] refPts = new double[2]; + + // Linear operations (normalization) applied before NormalizedTransform. + srcPts[0] -= FE; + srcPts[1] -= FN; + srcPts[0] /= (k0 * 2*R); + srcPts[1] /= (k0 * 2*R); + + // The non-linear part of map projection (the "kernel"). + createNormalizedProjection(false); + ((NormalizedProjection) transform).inverseTransform(srcPts, 0, refPts, 0); + + // Linear operations (denormalization) applied after NormalizedTransform. + refPts[0] /= n; + refPts[0] = toDegrees(refPts[0] + λ0); + refPts[1] = toDegrees(refPts[1]); + + // Transform the same point, now using the spherical implementation. + ObliqueStereographic spherical = (ObliqueStereographic) transform; + spherical = new ObliqueStereographic.Spherical(spherical); + spherical.inverseTransform(srcPts, 0, dstPts, 0); + + // Linear operations (denormalization) applied after NormalizedTransform. + dstPts[0] /= n; + dstPts[0] = toDegrees(dstPts[0] + λ0); + dstPts[1] = toDegrees(dstPts[1]); + + // Use a smaller tolerance because spherical and elliptical formulas should be equivalent in this case. + assertArrayEquals("Spherical inverse projection", refPts, dstPts, Formulas.ANGULAR_TOLERANCE / 1E6); + } + + /** + * Verifies the consistency of spherical formulas with the elliptical formulas. + * This test computes the derivative at a point and takes the result of the elliptical + * implementation as a reference. + * + * @throws TransformException if an error occurred while computing the derivative. + */ + @Test + @DependsOnMethod("testDerivative") + public void testSphericalDerivative() throws TransformException { + final double[] srcPts = new double[] {λt, φt}; // in degrees + srcPts[0] = toRadians(srcPts[0]) - λ0; + srcPts[1] = toRadians(srcPts[1]); + srcPts[0] *= n; + + // Using elliptical implementation. + createNormalizedProjection(false); + final Matrix reference = ((NormalizedProjection) transform).transform(srcPts, 0, null, 0, true); + + // Using spherical implementation. + ObliqueStereographic spherical = (ObliqueStereographic) transform; + spherical = new ObliqueStereographic.Spherical(spherical); + final Matrix derivative = spherical.transform(srcPts, 0, null, 0, true); + + tolerance = 1E-12; + assertMatrixEquals("Spherical derivative", reference, derivative, + new Matrix2(tolerance, tolerance, + tolerance, tolerance)); + } + + /** + * Creates a projection and derivates a few points. + * + * @throws TransformException Should never happen. + */ + @Test + public void testDerivative() throws TransformException { + createNormalizedProjection(true); + tolerance = 1E-9; + + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + derivativeDeltas = new double[] {delta, delta}; + verifyDerivative(toRadians( 0), toRadians( 0)); + verifyDerivative(toRadians(-3), toRadians(30)); + verifyDerivative(toRadians(+6), toRadians(60)); + } } Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -70,7 +70,7 @@ public final strictfp class PolarStereog @Test public void testSphericalCaseSouth() throws FactoryException, TransformException { createNormalizedProjection(new PolarStereographicSouth()); - final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. derivativeDeltas = new double[] {delta, delta}; verifyInDomain(CoordinateDomain.GEOGRAPHIC_RADIANS_SOUTH, 56763886); } Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java?rev=1709462&r1=1709461&r2=1709462&view=diff ============================================================================== --- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java [UTF-8] (original) +++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java [UTF-8] Mon Oct 19 18:17:37 2015 @@ -117,7 +117,7 @@ public final strictfp class TransverseMe createNormalizedProjection(false, 0); tolerance = 1E-9; - final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. derivativeDeltas = new double[] {delta, delta}; verifyDerivative(toRadians( 0), toRadians( 0)); verifyDerivative(toRadians(-3), toRadians(30)); @@ -134,7 +134,7 @@ public final strictfp class TransverseMe createNormalizedProjection(true, 0); tolerance = 1E-9; - final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. derivativeDeltas = new double[] {delta, delta}; verifyDerivative(toRadians( 0), toRadians( 0)); verifyDerivative(toRadians(-3), toRadians(30));