You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2015/07/13 18:22:26 UTC
svn commit: r1690747 - in /sis/branches/JDK8/core/sis-referencing/src:
main/java/org/apache/sis/referencing/operation/projection/
test/java/org/apache/sis/referencing/operation/projection/
Author: desruisseaux
Date: Mon Jul 13 16:22:25 2015
New Revision: 1690747
URL: http://svn.apache.org/r1690747
Log:
Review Lambert Conic Conformal formulas in light of §1.3.1.1 in
IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015.
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/GeneralLambert.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConformal.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/GeneralLambertTest.java
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConformalTest.java
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/GeneralLambert.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/GeneralLambert.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/GeneralLambert.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/GeneralLambert.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -167,7 +167,7 @@ abstract class GeneralLambert extends No
+ c4χ * sin(4*φ)
+ c2χ * sin(2*φ);
/*
- * Note: a previous version checked if the value of the smallest term c8χ·sin(8φ) was smaller than
+ * Note: a previous version checked if the value of the smallest term c8χ⋅sin(8φ) was smaller than
* the iteration tolerance. But this was not reliable enough. We use now a hard coded threshold
* determined empirically by MercatorMethodComparison.
*/
Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConformal.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConformal.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConformal.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConformal.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -43,6 +43,7 @@ import org.apache.sis.util.Debug;
import static java.lang.Math.*;
import static java.lang.Double.*;
+import static org.apache.sis.math.MathFunctions.isPositive;
/**
@@ -88,6 +89,8 @@ public class LambertConformal extends Ge
/**
* Internal coefficients for computation, depending only on values of standards parallels.
+ * This is defined as {@literal n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂)} in §1.3.1.1 of
+ * IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015.
*/
final double n;
@@ -177,6 +180,33 @@ public class LambertConformal extends Ge
throw new IllegalArgumentException(Errors.format(Errors.Keys.LatitudesAreOpposite_2,
new Latitude(φ1), new Latitude(φ2)));
}
+ /*
+ * Whether to favorize precision in the North hemisphere (true) or in the South hemisphere (false).
+ * The IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015 uses the
+ * following formula:
+ *
+ * t = tan(π/4 – φ/2) / [(1 – ℯ⋅sinφ)/(1 + ℯ⋅sinφ)] ^ (ℯ/2)
+ *
+ * while our 'expOfNorthing' function is defined like above, but with tan(π/4 + φ/2) instead of
+ * tan(π/4 - φ/2). Those two expressions are the reciprocal of each other if we reverse the sign of
+ * φ (see 'expOfNorthing' for trigonometric identities), but their accuracies are not equivalent:
+ * the hemisphere having values closer to zero is favorized. The EPSG formulas favorize the North
+ * hemisphere.
+ *
+ * Since Apache SIS's formula uses the + operator instead of -, we need to reverse the sign of φ
+ * values in order to match the EPSG formulas, but we will do that only if the map projection is
+ * for the North hemisphere.
+ *
+ * TEST: whether 'isNorth' is true of false does not change the formulas "correctness": it is only
+ * a small accuracy improvement. One can safely force this boolean value to 'true' or 'false' for
+ * testing purpose.
+ */
+ final boolean isNorth = isPositive(φ0);
+ if (isNorth) {
+ φ0 = -φ0;
+ φ1 = -φ1;
+ φ2 = -φ2;
+ }
φ0 = toRadians(φ0);
φ1 = toRadians(φ1);
φ2 = toRadians(φ2);
@@ -189,6 +219,11 @@ public class LambertConformal extends Ge
final double sinφ1 = sin(φ1);
final double m1 = cos(φ1) / rν(sinφ1);
final double t1 = expOfNorthing(φ1, excentricity*sinφ1);
+ /*
+ * Computes n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as ln(m₁/m₂) / ln(t₁/t₂)
+ * since division is less at risk of precision lost than subtraction. Note that this equation
+ * tends toward 0/0 if φ₁ ≈ φ₂, which force us to do a special check for the SP1 case.
+ */
if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) { // Should be 'true' for 2SP case.
final double sinφ2 = sin(φ2);
final double m2 = cos(φ2) / rν(sinφ2);
@@ -198,18 +233,26 @@ public class LambertConformal extends Ge
n = -sinφ1;
}
/*
+ * Computes F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
* Following constants will be stored in the denormalization matrix, to be applied after
* the non-linear formulas implemented by this LambertConformal class. Opportunistically
* use double-double arithmetic since we the matrix coefficients will be stored in this
* format anyway. This makes a change in the 2 or 3 last digits.
*/
- final DoubleDouble F = new DoubleDouble(-pow(t1, -n), 0);
- F.multiply(m1, 0);
- F.divide(n, 0);
- final DoubleDouble ρ0 = new DoubleDouble(); // Initialized to zero.
+ final DoubleDouble F = new DoubleDouble(n, 0);
+ F.multiply(pow(t1, n), 0);
+ F.inverseDivide(m1, 0);
+ if (!isNorth) {
+ F.negate();
+ }
+ /*
+ * Compute r = a⋅F⋅tⁿ from EPSG notes where (in our case) a=1 and t is our 'expOfNorthing' function.
+ * Note that Snyder calls this term "ρ0".
+ */
+ final DoubleDouble r0 = new DoubleDouble(); // Initialized to zero.
if (φ0 != copySign(PI/2, -n)) { // For avoiding the rounding error documented in expOfNorthing(+π/2).
- ρ0.value = pow(expOfNorthing(φ0, excentricity*sin(φ0)), n);
- ρ0.multiply(F);
+ r0.value = pow(expOfNorthing(φ0, excentricity*sin(φ0)), n);
+ r0.multiply(F);
}
/*
* At this point, all parameters have been processed. Now store
@@ -228,13 +271,16 @@ public class LambertConformal extends Ge
* - Multiply by the scale factor (done by the super-class constructor).
* - Add false easting and false northing (done by the super-class constructor).
*/
- context.getMatrix(true).convertAfter(0, new DoubleDouble(-n, 0), // Multiplication factor for longitudes.
+ final MatrixSIS normalize = context.getMatrix(true);
+ normalize.convertAfter(0, new DoubleDouble(isNorth ? n : -n, 0), // Multiplication factor for longitudes.
(type == BELGIUM) ? new DoubleDouble(-BELGE_A, 0) : null); // Longitude translation for Belgium.
-
+ if (isNorth) {
+ normalize.convertAfter(1, new DoubleDouble(-1, 0), null);
+ }
final MatrixSIS denormalize = context.getMatrix(false);
denormalize.convertBefore(0, F, null);
F.negate();
- denormalize.convertBefore(1, F, ρ0);
+ denormalize.convertBefore(1, F, r0);
}
/**
@@ -301,7 +347,7 @@ public class LambertConformal extends Ge
}
/**
- * Converts the specified (λ,φ) coordinate (units in radians) and stores the result in {@code dstPts}
+ * Converts the specified (θ,φ) coordinate (units in radians) and stores the result in {@code dstPts}
* (linear distance on a unit sphere). In addition, opportunistically computes the projection derivative
* if {@code derivate} is {@code true}.
*
@@ -319,24 +365,24 @@ public class LambertConformal extends Ge
* 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];
- final double φ = srcPts[srcOff + 1];
+ final double θ = srcPts[srcOff]; // θ = λ⋅n
+ final double φ = srcPts[srcOff + 1]; // Sign may be reversed
final double absφ = abs(φ);
- final double sinλ = sin(λ);
- final double cosλ = cos(λ);
+ final double sinθ = sin(θ);
+ final double cosθ = cos(θ);
final double sinφ;
- final double ρ; // Snyder p. 108
+ final double r; // From EPSG guide. Note that Snyder p. 108 calls this term "ρ".
if (absφ < PI/2) {
sinφ = sin(φ);
- ρ = pow(expOfNorthing(φ, excentricity*sinφ), n);
+ r = pow(expOfNorthing(φ, excentricity*sinφ), n);
} else if (absφ < PI/2 + ANGULAR_TOLERANCE) {
sinφ = 1;
- ρ = (φ*n >= 0) ? POSITIVE_INFINITY : 0;
+ r = (φ*n >= 0) ? POSITIVE_INFINITY : 0;
} else {
- ρ = sinφ = NaN;
+ r = sinφ = NaN;
}
- final double x = ρ * sinλ;
- final double y = ρ * cosλ;
+ final double x = r * sinθ;
+ final double y = r * cosθ;
if (dstPts != null) {
dstPts[dstOff ] = x;
dstPts[dstOff + 1] = y;
@@ -349,12 +395,12 @@ public class LambertConformal extends Ge
//
final double dρ;
if (sinφ != 1) {
- dρ = n * dy_dφ(sinφ, cos(φ)) * ρ;
+ dρ = n * dy_dφ(sinφ, cos(φ)) * r;
} else {
- dρ = ρ;
+ dρ = r;
}
- return new Matrix2(y, dρ*sinλ, // ∂x/∂λ , ∂x/∂φ
- -x, dρ*cosλ); // ∂y/∂λ , ∂y/∂φ
+ return new Matrix2(y, dρ*sinθ, // ∂x/∂λ , ∂x/∂φ
+ -x, dρ*cosθ); // ∂y/∂λ , ∂y/∂φ
}
/**
@@ -431,31 +477,31 @@ public class LambertConformal extends Ge
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- final double λ = srcPts[srcOff];
- final double φ = srcPts[srcOff + 1];
+ final double θ = srcPts[srcOff]; // θ = λ⋅n
+ final double φ = srcPts[srcOff + 1]; // Sign may be reversed
final double absφ = abs(φ);
- final double sinλ = sin(λ);
- final double cosλ = cos(λ);
- final double ρ;
+ final double sinθ = sin(θ);
+ final double cosθ = cos(θ);
+ final double r; // Snyder p. 108 calls this term "ρ", but we use "r" for consistency with EPSG guide.
if (absφ < PI/2) {
- ρ = pow(tan(PI/4 + 0.5*φ), n);
+ r = pow(tan(PI/4 + 0.5*φ), n);
} else if (absφ < PI/2 + ANGULAR_TOLERANCE) {
- ρ = (φ*n >= 0) ? POSITIVE_INFINITY : 0;
+ r = (φ*n >= 0) ? POSITIVE_INFINITY : 0;
} else {
- ρ = NaN;
+ r = NaN;
}
- final double x = ρ * sinλ;
- final double y = ρ * cosλ;
+ final double x = r * sinθ;
+ final double y = r * cosθ;
Matrix derivative = null;
if (derivate) {
final double dρ;
if (absφ < PI/2) {
- dρ = n*ρ / cos(φ);
+ dρ = n*r / cos(φ);
} else {
dρ = NaN;
}
- derivative = new Matrix2(y, dρ*sinλ, // ∂x/∂λ , ∂x/∂φ
- -x, dρ*cosλ); // ∂y/∂λ , ∂y/∂φ
+ derivative = new Matrix2(y, dρ*sinθ, // ∂x/∂λ , ∂x/∂φ
+ -x, dρ*cosθ); // ∂y/∂λ , ∂y/∂φ
}
// Following part is common to all spherical projections: verify, store and return.
assert Assertions.checkDerivative(derivative, super.transform(srcPts, srcOff, dstPts, dstOff, derivate))
@@ -479,7 +525,7 @@ public class LambertConformal extends Ge
double y = srcPts[srcOff+1];
final double ρ = hypot(x, y);
x = atan2(x, y); // Really (x,y), not (y,x)
- y = 2 * atan(pow(1/ρ, -1/n)) - PI/2;
+ y = 2*atan(pow(1/ρ, -1/n)) - PI/2;
assert checkInverseTransform(srcPts, srcOff, dstPts, dstOff, x, y);
dstPts[dstOff ] = x;
dstPts[dstOff + 1] = y;
@@ -491,11 +537,11 @@ public class LambertConformal extends Ge
*/
private boolean checkInverseTransform(final double[] srcPts, final int srcOff,
final double[] dstPts, final int dstOff,
- final double λ, final double φ)
+ final double θ, final double φ)
throws ProjectionException
{
super.inverseTransform(srcPts, srcOff, dstPts, dstOff);
- return Assertions.checkInverseTransform(dstPts, dstOff, λ, φ);
+ return Assertions.checkInverseTransform(dstPts, dstOff, θ, φ);
}
}
Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -153,6 +153,9 @@ public abstract class NormalizedProjecti
/**
* Maximum number of iterations for iterative computations.
+ * The iterative methods used in subclasses should converge quickly (in 3 or 4 iterations)
+ * when used for a planet with an excentricity similar to Earth. But we allow a high limit
+ * in case someone uses SIS for some planet with higher excentricity.
*/
static final int MAXIMUM_ITERATIONS = 15;
@@ -865,6 +868,9 @@ public abstract class NormalizedProjecti
* of the value returned by this method. This function is <em>almost</em> the converse of
* {@link GeneralLambert#φ(double)}.
*
+ * <p>In IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015,
+ * a function closely related to this one has the letter <var>t</var>.</p>
+ *
*
* <div class="section">Properties</div>
* This function is used with φ values in the [-π/2 … π/2] range and has a periodicity of 2π.
@@ -914,6 +920,10 @@ public abstract class NormalizedProjecti
final double expOfNorthing(final double φ, final double ℯsinφ) {
/*
* Note: tan(π/4 - φ/2) = 1 / tan(π/4 + φ/2)
+ *
+ * A + sign in the equation favorises slightly the accuracy in South hemisphere, while a - sign
+ * favorises slightly the North hemisphere (but the differences are very small). In Apache SIS,
+ * we handle that by changing the sign of some terms in the (de)normalisation matrices.
*/
return tan(PI/4 + 0.5*φ) * pow((1 - ℯsinφ) / (1 + ℯsinφ), 0.5*excentricity);
}
Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/GeneralLambertTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/GeneralLambertTest.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/GeneralLambertTest.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/GeneralLambertTest.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -130,7 +130,7 @@ public final strictfp class GeneralLambe
final double byIterativeMethod = comparator.byIterativeMethod(t);
final double bySeriesExpansion = comparator.bySeriesExpansion(t);
assertEquals(bySeriesExpansion, byIterativeMethod, 1E-11);
- assertEquals(bySeriesExpansion, projection.φ(t), 1E-15);
+ assertEquals(bySeriesExpansion, projection.φ(t), 1E-15); // Tolerance threshold close to 1 ULP of 2π.
}
}
}
Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConformalTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConformalTest.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConformalTest.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConformalTest.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -93,25 +93,25 @@ public final strictfp class LambertConfo
@Test
public void testSpecialLatitudes() throws ProjectionException {
if (transform == null) { // May have been initialized by 'testSphericalCase'.
- initialize(true, 40); // Elliptical case
+ initialize(true, 40); // Elliptical case
}
final double INF = POSITIVE_INFINITY;
assertEquals ("Not a number", NaN, transform(NaN), NORMALIZED_TOLERANCE);
assertEquals ("Out of range", NaN, transform(+2), NORMALIZED_TOLERANCE);
assertEquals ("Out of range", NaN, transform(-2), NORMALIZED_TOLERANCE);
assertEquals ("Forward 0°N", 1, transform(0), NORMALIZED_TOLERANCE);
- assertEquals ("Forward 90°N", 0, transform(+PI/2), NORMALIZED_TOLERANCE);
- assertEquals ("Forward 90°S", INF, transform(-PI/2), NORMALIZED_TOLERANCE);
- assertEquals ("Forward (90+ε)°N", 0, transform(+nextUp(+PI/2)), NORMALIZED_TOLERANCE);
- assertEquals ("Forward (90+ε)°S", INF, transform(-nextUp( PI/2)), NORMALIZED_TOLERANCE);
- assertEquals ("Forward (90-ε)°N", 0, transform(-nextUp(-PI/2)), 1E-10);
+ assertEquals ("Forward 90°S", 0, transform(-PI/2), NORMALIZED_TOLERANCE);
+ assertEquals ("Forward 90°N", INF, transform(+PI/2), NORMALIZED_TOLERANCE);
+ assertEquals ("Forward (90+ε)°S", 0, transform(-nextUp( PI/2)), NORMALIZED_TOLERANCE);
+ assertEquals ("Forward (90+ε)°N", INF, transform(+nextUp(+PI/2)), NORMALIZED_TOLERANCE);
+ assertEquals ("Forward (90-ε)°S", 0, transform(+nextUp(-PI/2)), 1E-10);
assertEquals ("Not a number", NaN, inverseTransform(NaN), NORMALIZED_TOLERANCE);
- assertEquals ("Inverse 0", +PI/2, inverseTransform( 0), NORMALIZED_TOLERANCE);
+ assertEquals ("Inverse 0", -PI/2, inverseTransform( 0), NORMALIZED_TOLERANCE);
assertEquals ("Inverse +1", 0, inverseTransform(+1), NORMALIZED_TOLERANCE);
assertEquals ("Inverse -1", 0, inverseTransform(-1), NORMALIZED_TOLERANCE);
- assertEquals ("Inverse +∞", -PI/2, inverseTransform(INF), NORMALIZED_TOLERANCE);
- assertEquals ("Inverse -∞", -PI/2, inverseTransform(-INF), NORMALIZED_TOLERANCE);
+ assertEquals ("Inverse +∞", +PI/2, inverseTransform(INF), NORMALIZED_TOLERANCE);
+ assertEquals ("Inverse -∞", +PI/2, inverseTransform(-INF), NORMALIZED_TOLERANCE);
// Like the north case, but with sign inversed.
initialize(((LambertConformal) transform).excentricity != 0, -40);
Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java?rev=1690747&r1=1690746&r2=1690747&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java [UTF-8] Mon Jul 13 16:22:25 2015
@@ -26,7 +26,7 @@ import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.internal.metadata.ReferencingServices;
-import static java.lang.Math.*;
+import static java.lang.Math.*; // Not StrictMath in this particular case.
/**
@@ -215,7 +215,7 @@ public final class MercatorMethodCompari
}
out.println("Comparison of different ways to compute φ for excentricity " + excentricity + '.');
out.println("Values are in units of " + NormalizedProjection.ITERATION_TOLERANCE + " radians (about "
- + Math.round(toDegrees(NormalizedProjection.ITERATION_TOLERANCE) * 60 * ReferencingServices.NAUTICAL_MILE * 1000)
+ + round(toDegrees(NormalizedProjection.ITERATION_TOLERANCE) * 60 * ReferencingServices.NAUTICAL_MILE * 1000)
+ " mm on Earth).");
final StatisticsFormat format = StatisticsFormat.getInstance();
format.setBorderWidth(1);