This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
commit 45813913866a31d4a14d812b5a4bc3723f9ca49f
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Tue May 3 14:30:51 2022 +0200
Slight modification of `AuthalicConversion` parameters for making clearer
that
this is a conversion between geodetic latitude φ and authalic latitude β.
---
.../operation/projection/AlbersEqualArea.java | 12 +++----
.../operation/projection/AuthalicConversion.java | 37 ++++++++++++++++------
.../operation/projection/CylindricalEqualArea.java | 2 +-
.../projection/LambertAzimuthalEqualArea.java | 22 +++++--------
.../projection/AuthalicConversionTest.java | 6 ++--
5 files changed, 46 insertions(+), 33 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index d32c47c518..4fbe7b7155 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -50,7 +50,7 @@ import static
org.apache.sis.internal.referencing.provider.AlbersEqualArea.*;
*
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.1
+ * @version 1.2
* @since 0.8
* @module
*/
@@ -274,12 +274,12 @@ public class AlbersEqualArea extends AuthalicConversion {
final double x = srcPts[srcOff ];
final double y = srcPts[srcOff+1];
/*
- * Note: Snyder suggests to reverse the sign of x, y and ρ₀ if n is
negative. It should not done in Apache SIS
- * implementation because (x,y) are premultiplied by n (by the
normalization affine transform) before to enter
- * in this method, so if n was negative those values have already
their sign reverted.
+ * Note: Snyder suggests to reverse the sign of x, y and ρ₀ if n is
negative. It should not be done
+ * in SIS implementation because (x,y) are premultiplied by n (by the
normalization affine transform)
+ * before to enter in this method, so if n was negative those values
have already their sign reverted.
*/
dstPts[dstOff ] = atan2(x, y);
- dstPts[dstOff+1] = φ((C - (x*x + y*y)) / nm);
+ dstPts[dstOff+1] = φ((C - (x*x + y*y)) / (nm*qmPolar));
/*
* Note: Snyder 14-19 gives q = (C - ρ²n²/a²)/n where ρ = √(x² +
(ρ₀ - y)²).
* But in Apache SIS implementation, ρ₀ has already been subtracted by
the matrix before we reach this point.
@@ -290,7 +290,7 @@ public class AlbersEqualArea extends AuthalicConversion {
* q = (C - (x² + y²)) / n
*
* We divide by nm instead of n, so a (1-ℯ²) term is missing. But that
missing term will be cancelled with
- * the missing (1-ℯ²) term in qmPolar (the divisor applied by the
φ(double) method that we invoke).
+ * the missing (1-ℯ²) term in `qmPolar`. The division by `qmPolar` is
for converting y to sin(β).
*/
}
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
index 38a373aace..e96a9d3c4a 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
@@ -184,7 +184,7 @@ abstract class AuthalicConversion extends
NormalizedProjection {
* We pay the cost of checking for the spherical case in each method
invocation because otherwise,
* users creating their own map projection subclasses could get a
non-working implementation.
*
- * @param sinφ the sine of the latitude <var>q</var> is calculated for.
+ * @param sinφ the sine of the geodetic latitude for which <var>q</var>
is calculated.
* @return <var>q</var> from Snyder equation (3-12).
*/
final double qm(final double sinφ) {
@@ -194,7 +194,7 @@ abstract class AuthalicConversion extends
NormalizedProjection {
*/
if (isSpherical) return 2*sinφ;
final double ℯsinφ = eccentricity * sinφ;
- return sinφ / (1 - ℯsinφ*ℯsinφ) + atanh(ℯsinφ) / eccentricity;
+ return sinφ/(1 - ℯsinφ*ℯsinφ) + atanh(ℯsinφ)/eccentricity;
}
/**
@@ -211,6 +211,20 @@ abstract class AuthalicConversion extends
NormalizedProjection {
return 2*cosφ / (t*t);
}
+ /**
+ * Converts the sine of geodetic latitude to the sin of authalic latitude.
+ * This is defined by {@code qm(sinφ) / qmPolar}.
+ *
+ * @param sinφ the sine of the geodetic latitude.
+ * @return the sine of the authalic latitude.
+ */
+ final double sinβ(double sinφ) {
+ // Edited copy of `qm(double)` method.
+ if (isSpherical) return sinφ;
+ sinφ *= eccentricity; // Become `ℯsinφ` from here.
+ return (sinφ/(1 - sinφ*sinφ) + atanh(sinφ)) / (eccentricity * qmPolar);
+ }
+
/**
* Computes the latitude using equation 3-18 from Snyder, followed by
iterative resolution of Snyder 3-16.
* In theory, the series expansion given by equation 3-18 (φ ≈ c₂⋅sin(2β)
+ c₄⋅sin(4β) + c₈⋅sin(8β)) should
@@ -218,15 +232,18 @@ abstract class AuthalicConversion extends
NormalizedProjection {
* have a sufficient amount of terms for achieving the centimetric
precision, so we "finish" it by the
* iterative method. The series expansion is nevertheless useful for
reducing the number of iterations.
*
+ * <h4>Relationship with northing</h4>
+ * The simplest projection using this formula is the {@link
CylindricalEqualArea} projection.
+ * In that case, sin(β) = <var>y</var> / {@link #qmPolar}.
+ *
* <h4>Spherical case</h4>
- * In the spherical case, this method returns {@code asin(y/2)}. This
method does not use that
- * simplification for the spherical case. This optimization is left to the
caller if desired.
+ * In the spherical case, this method returns {@code β = asin(sinβ)}. This
method does not check for
+ * that simplification for the spherical case. This optimization is left
to the caller if desired.
*
- * @param y in the cylindrical case, this is northing on the normalized
ellipsoid.
- * @return the latitude in radians.
+ * @param sinβ sine of the authalic latitude.
+ * @return the geodetic latitude in radians.
*/
- final double φ(final double y) throws ProjectionException {
- final double sinβ = y / qmPolar;
+ final double φ(final double sinβ) throws ProjectionException {
final double sinβ2 = sinβ * sinβ;
final double β = asin(sinβ);
/*
@@ -244,6 +261,7 @@ abstract class AuthalicConversion extends
NormalizedProjection {
* as q = (C - ρ²)/n (see comment in
AlbersEqualArea.inverseTransform(…) for the mathematic).
* The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) =
q/(1-ℯ²), the desired value.
*/
+ final double y = qmPolar * sinβ;
for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
final double sinφ = sin(φ);
final double cosφ = cos(φ);
@@ -277,9 +295,10 @@ abstract class AuthalicConversion extends
NormalizedProjection {
*/
final double as = abs(sinβ);
if (abs(as - 1) < ANGULAR_TOLERANCE) {
+ final double y = qmPolar * sinβ; // Do not use β because it
may be NaN.
return copySign(PI/2, y); // Value is at a pole.
}
- if (as >= 1 || Double.isNaN(y)) {
+ if (!(as < 1)) { // Use `!` for catching
NaN.
return Double.NaN; // Value "after" the pole.
}
// Value should have converged but did not.
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
index ed85a6e3d7..24c2877a3c 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
@@ -294,7 +294,7 @@ public class CylindricalEqualArea extends
AuthalicConversion {
{
final double y = srcPts[srcOff+1]; // Must be before
writing x.
dstPts[dstOff ] = srcPts[srcOff ]; // Must be before
writing y.
- dstPts[dstOff+1] = φ(y);
+ dstPts[dstOff+1] = φ(y / qmPolar);
/*
* Equation 10-26 of Snyder gives β = asin(2y⋅k₀/(a⋅qPolar)).
* In our case it simplifies to sinβ = (y/qmPolar) because:
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
index c56446e768..4901bfeeb8 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
@@ -108,7 +108,7 @@ public class LambertAzimuthalEqualArea extends
AuthalicConversion {
final double sinφ0 = sin(φ0);
final double cosφ0 = cos(φ0);
polar = abs(sinφ0) == 1; // sin(φ) == 1 implies a tolerance
∆φ≈1E-6°.
- sinß0 = qm(sinφ0) / qmPolar;
+ sinß0 = sinβ(sinφ0);
cosß0 = sqrt(1 - sinß0*sinß0);
/*
* In the polar case we have cos(φ₀) ≈ 0 and cos(ß₀) ≈ 0, which cause
D = 0/0.
@@ -180,16 +180,11 @@ public class LambertAzimuthalEqualArea extends
AuthalicConversion {
final double sinλ = sin(λ);
final double cosλ = cos(λ);
final double sinφ = sin(φ);
- final double qm = qm(sinφ);
if (!polar) {
/*
* Note: in the spherical case, ß = φ (ß is the authalic radius).
- * We could do an optimization, but it would be a cost for the
ellipsoidal
- * case without saving a lot for the spherical case. The general
path is:
- *
- * sinß = qm(sinφ) / qmPolar = 2*sinφ / 2 = sinφ
*/
- final double sinß = qm / qmPolar;
+ final double sinß = sinβ(sinφ);
final double cosß = sqrt(1 - sinß*sinß);
final double c = sinß0*sinß + cosß0*cosß*cosλ + 1;
/*
@@ -232,7 +227,7 @@ public class LambertAzimuthalEqualArea extends
AuthalicConversion {
* affine transform (which cause the sign of q to be also reversed).
After the transform, sign of
* y will be reversed by the denormalize affine transform, so it
should not be reversed here.
*/
- double ρ = sqrt(qmPolar + qm);
+ double ρ = sqrt(qmPolar + qm(sinφ));
if (sinφ == 1) {
ρ = Double.NaN; // Antipodal point.
}
@@ -262,25 +257,24 @@ public class LambertAzimuthalEqualArea extends
AuthalicConversion {
{
double x = srcPts[srcOff ];
double y = srcPts[srcOff+1];
- double sinß_qp; // sin(ß) × qmPolar
+ double sinß;
if (polar) {
- sinß_qp = (x*x + y*y) - qmPolar;
+ sinß = (x*x + y*y)/qmPolar - 1;
} else {
final double ρ = fastHypot(x, y);
final double C = 2*asin(ρ / 2);
final double cosC = cos(C);
final double sinC = sin(C);
final double ysinC = y*sinC;
- sinß_qp = cosC*sinß0;
+ sinß = cosC*sinß0;
if (cosC != 1) { // cos(C) == 1 implies y/ρ =
0/0.
- sinß_qp += ysinC*cosß0/ρ;
+ sinß += ysinC*cosß0/ρ;
}
- sinß_qp *= qmPolar;
y = ρ*cosC*cosß0 - ysinC*sinß0;
x *= sinC;
}
dstPts[dstOff ] = atan2(x, y);
- dstPts[dstOff+1] = isSpherical ? asin(sinß_qp/2) : φ(sinß_qp);
+ dstPts[dstOff+1] = isSpherical ? asin(sinß) : φ(sinß);
}
/*
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AuthalicConversionTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AuthalicConversionTest.java
index 213b7decd8..54a8e7e363 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AuthalicConversionTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AuthalicConversionTest.java
@@ -34,7 +34,7 @@ import static org.junit.Assert.assertEquals;
* Tests {@link AuthalicConversion}.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.2
* @since 1.0
* @module
*/
@@ -89,7 +89,7 @@ public final strictfp class AuthalicConversionTest extends
MapProjectionTestCase
for (int i=0; i<100; i++) {
final double y = random.nextDouble() * 3 - 1.5;
final double reference = reference(projection, y);
- final double actual = projection.φ(y);
+ final double actual = projection.φ(y / projection.qmPolar);
assertEquals(reference, actual,
NormalizedProjection.ITERATION_TOLERANCE);
}
}
@@ -116,7 +116,7 @@ public final strictfp class AuthalicConversionTest extends
MapProjectionTestCase
final CylindricalEqualArea projection = new
CylindricalEqualArea(provider, parameters);
for (double y = -1.25; y <= 1.25; y += 0.01) {
final double reference = reference(projection, y);
- final double actual = projection.φ(y);
+ final double actual = projection.φ(y / projection.qmPolar);
if (abs(actual - reference) >
NormalizedProjection.ANGULAR_TOLERANCE) {
System.out.println("Error exceeds tolerance threshold at
eccentricity " + e);
return;
