Changeset 13632 in josm

Timestamp:

2018-04-15T14:46:01+02:00 (7 years ago)

Author:

Don-vip

Message:

see #16129 - add spherical versions of Cassini and Albers projections

Location:

trunk

Files:

• data_nodist/projection/projection-regression-test-data (modified) (1 diff)
• src/org/openstreetmap/josm/data/projection/proj/AlbersEqualArea.java (modified) (5 diffs)
• src/org/openstreetmap/josm/data/projection/proj/CassiniSoldner.java (modified) (2 diffs)

trunk/data_nodist/projection/projection-regression-test-data

@@ -18831 +18831 @@
 EPSG:53028
   ll  -35.72259959292767 0.08011233833985498
-  en  7232.054410449118 -3972172.922872512
-  ll2 -35.72259959292767 0.08011233833986359
+  en  7232.0544104483415 -3972172.9228725126
+  ll2 -35.722603125469135 0.08011233833985497
 EPSG:5303
   ll  59.74559228260094 92.1133124256515

trunk/src/org/openstreetmap/josm/data/projection/proj/AlbersEqualArea.java

@@ -64 +64 @@
      * Constants used by the spherical and elliptical Albers projection.
      */
-    private double n, c, rho0;
+    private double n, n2, c, rho0;
 
     /**
@@ -116 +116 @@
         ec = 1.0 - .5 * (1.0-e2) *
              Math.log((1.0 - e) / (1.0 + e)) / e;
+        n2 = n * n;
         this.n = n;
     }
@@ -121 +122 @@
     @Override
     public double[] project(double y, double x) {
-        x *= n;
-        double rho = c - n * qsfn(Math.sin(y));
+        double theta = n * x;
+        double rho = c - (spherical ? n2 * Math.sin(y) : n * qsfn(Math.sin(y)));
         if (rho < 0.0) {
             if (rho > -EPSILON) {
@@ -132 +133 @@
         rho = Math.sqrt(rho) / n;
         // CHECKSTYLE.OFF: SingleSpaceSeparator
-        y = rho0 - rho * Math.cos(x);
-        x =        rho * Math.sin(x);
+        y = rho0 - rho * Math.cos(theta);
+        x =        rho * Math.sin(theta);
         // CHECKSTYLE.ON: SingleSpaceSeparator
         return new double[] {x, y};
@@ -150 +151 @@
             x = Math.atan2(x, y) / n;
             y = rho * n;
-            y = (c - y*y) / n;
-            if (Math.abs(y) <= ec) {
-                y = phi1(y);
+            if (spherical) {
+                y = aasin((c - y*y) / n2);
             } else {
-                y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0;
+                y = (c - y*y) / n;
+                if (Math.abs(y) <= ec) {
+                    y = phi1(y);
+                } else {
+                    y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0;
+                }
             }
         } else {

trunk/src/org/openstreetmap/josm/data/projection/proj/CassiniSoldner.java

@@ -33 +33 @@
 
     /**
-     * Contants used for the forward and inverse transform for the eliptical
+     * Latitude of origin.
+     */
+    private double phi0;
+
+    /**
+     * Constants used for the forward and inverse transform for the elliptical
      * case of the Cassini-Soldner.
      */
@@ -57 +62 @@
         if (params.lat0 == null)
             throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
-        double latitudeOfOrigin = Utils.toRadians(params.lat0);
-        ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin));
+        phi0 = Utils.toRadians(params.lat0);
+        ml0 = mlfn(phi0, Math.sin(phi0), Math.cos(phi0));
     }
 
     @Override
     public double[] project(double phi, double lam) {
-        double sinphi = Math.sin(phi);
-        double cosphi = Math.cos(phi);
+        if (spherical) {
+            double x = aasin(Math.cos(phi) * Math.sin(lam));
+            double y = Math.atan2(Math.tan(phi), Math.cos(lam));
+            return new double[] {x, y};
+        } else {
+            double sinphi = Math.sin(phi);
+            double cosphi = Math.cos(phi);
 
-        double n = 1.0 / (Math.sqrt(1.0 - e2 * sinphi * sinphi));
-        double tn = Math.tan(phi);
-        double t = tn * tn;
-        double a1 = lam * cosphi;
-        double c = cosphi * cosphi * e2 / (1 - e2);
-        double a2 = a1 * a1;
+            double n = 1.0 / (Math.sqrt(1.0 - e2 * sinphi * sinphi));
+            double tn = Math.tan(phi);
+            double t = tn * tn;
+            double a1 = lam * cosphi;
+            double c = cosphi * cosphi * e2 / (1 - e2);
+            double a2 = a1 * a1;
 
-        double x = n * a1 * (1.0 - a2 * t * (C1 - (8.0 - t + 8.0 * c) * a2 * C2));
-        double y = mlfn(phi, sinphi, cosphi) - ml0 + n * tn * a2 * (0.5 + (5.0 - t + 6.0 * c) * a2 * C3);
-        return new double[] {x, y};
+            double x = n * a1 * (1.0 - a2 * t * (C1 - (8.0 - t + 8.0 * c) * a2 * C2));
+            double y = mlfn(phi, sinphi, cosphi) - ml0 + n * tn * a2 * (0.5 + (5.0 - t + 6.0 * c) * a2 * C3);
+            return new double[] {x, y};
+        }
     }
 
     @Override
     public double[] invproject(double x, double y) {
-        double ph1 = invMlfn(ml0 + y);
-        double tn = Math.tan(ph1);
-        double t = tn * tn;
-        double n = Math.sin(ph1);
-        double r = 1.0 / (1.0 - e2 * n * n);
-        n = Math.sqrt(r);
-        r *= (1.0 - e2) * n;
-        double dd = x / n;
-        double d2 = dd * dd;
-        double phi = ph1 - (n * tn / r) * d2 * (0.5 - (1.0 + 3.0 * t) * d2 * C3);
-        double lam = dd * (1.0 + t * d2 * (-C4 + (1.0 + 3.0 * t) * d2 * C5)) / Math.cos(ph1);
-        return new double[] {phi, lam};
+        if (spherical) {
+            double dd = y + phi0;
+            double phi = aasin(Math.sin(dd * Math.cos(x)));
+            double lam = Math.atan2(Math.tan(x), Math.cos(dd));
+            return new double[] {phi, lam};
+        } else {
+            double ph1 = invMlfn(ml0 + y);
+            double tn = Math.tan(ph1);
+            double t = tn * tn;
+            double n = Math.sin(ph1);
+            double r = 1.0 / (1.0 - e2 * n * n);
+            n = Math.sqrt(r);
+            r *= (1.0 - e2) * n;
+            double dd = x / n;
+            double d2 = dd * dd;
+            double phi = ph1 - (n * tn / r) * d2 * (0.5 - (1.0 + 3.0 * t) * d2 * C3);
+            double lam = dd * (1.0 + t * d2 * (-C4 + (1.0 + 3.0 * t) * d2 * C5)) / Math.cos(ph1);
+            return new double[] {phi, lam};
+        }
     }