Changeset 9565 in josm for trunk/src

Timestamp:

2016-01-22T16:50:40+01:00 (9 years ago)

Author:

bastiK

Message:

add 2 standard parallel & non-spherical variants for Mercator projection (see #12186)
(imports pieces of code from the Geotools project)

Location:

trunk/src/org/openstreetmap/josm/data/projection

Files:

• CustomProjection.java (modified) (2 diffs)
• Ellipsoid.java (modified) (7 diffs)
• Projections.java (modified) (1 diff)
• proj/AbstractProj.java (modified) (2 diffs)
• proj/IScaleFactorProvider.java (added)
• proj/Mercator.java (modified) (4 diffs)
• proj/ProjParameters.java (modified) (1 diff)

trunk/src/org/openstreetmap/josm/data/projection/CustomProjection.java

@@ -27 +27 @@
 import org.openstreetmap.josm.data.projection.datum.WGS84Datum;
 import org.openstreetmap.josm.data.projection.proj.ICentralMeridianProvider;
+import org.openstreetmap.josm.data.projection.proj.IScaleFactorProvider;
 import org.openstreetmap.josm.data.projection.proj.Mercator;
 import org.openstreetmap.josm.data.projection.proj.Proj;
@@ -274 +275 @@
                 this.k0 = parseDouble(s, Param.k_0.key);
             }
+            if (proj instanceof IScaleFactorProvider) {
+                this.k0 *= ((IScaleFactorProvider) proj).getScaleFactor();
+            }
             s = parameters.get(Param.bounds.key);
             if (s != null) {

trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java

@@ -132 +132 @@
 
     /**
+     * if ellipsoid is spherical, i.e. the major and minor semiaxis are
+     * the same
+     */
+    public final boolean spherical;
+
+    /**
      * private constructur - use one of the create_* methods
      *
@@ -139 +145 @@
      * @param e2 first eccentricity squared
      * @param eb2 square of the second eccentricity
-     */
-    private Ellipsoid(double a, double b, double e, double e2, double eb2) {
+     * @param sperical if the ellipsoid is sphere
+     */
+    private Ellipsoid(double a, double b, double e, double e2, double eb2, boolean sperical) {
         this.a = a;
         this.b = b;
@@ -146 +153 @@
         this.e2 = e2;
         this.eb2 = eb2;
+        this.spherical = sperical;
     }
 
@@ -159 +167 @@
         double e = Math.sqrt(e2);
         double eb2 = e2 / (1.0 - e2);
-        return new Ellipsoid(a, b, e, e2, eb2);
+        return new Ellipsoid(a, b, e, e2, eb2, a == b);
     }
 
@@ -173 +181 @@
         double e = Math.sqrt(es);
         double eb2 = es / (1.0 - es);
-        return new Ellipsoid(a, b, e, es, eb2);
+        return new Ellipsoid(a, b, e, es, eb2, es == 0);
     }
 
@@ -188 +196 @@
         double e = Math.sqrt(e2);
         double eb2 = e2 / (1.0 - e2);
-        return new Ellipsoid(a, b, e, e2, eb2);
+        return new Ellipsoid(a, b, e, e2, eb2, f == 0);
     }
 
@@ -334 +342 @@
         return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg));
     }
-    
+
     /**
      * convert ellipsoidal coordinates to cartesian coordinates

trunk/src/org/openstreetmap/josm/data/projection/Projections.java

@@ -89 +89 @@
         registerBaseProjection("aea", AlbersEqualArea.class, "core");
         registerBaseProjection("cass", CassiniSoldner.class, "core");
-        registerBaseProjection("josm:smerc", Mercator.class, "core");
         registerBaseProjection("laea", LambertAzimuthalEqualArea.class, "core");
         registerBaseProjection("lcc", LambertConformalConic.class, "core");
         registerBaseProjection("lonlat", LonLat.class, "core");
+        registerBaseProjection("merc", Mercator.class, "core");
         registerBaseProjection("omerc", ObliqueMercator.class, "core");
         registerBaseProjection("somerc", SwissObliqueMercator.class, "core");

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

@@ -74 +74 @@
      */
     protected double e2;
+    
+    /**
+     * is ellipsoid spherical?
+     * @see Ellisoid.spherical
+     */
+    protected boolean spherical;
 
     @Override
@@ -79 +85 @@
         e2 = params.ellps.e2;
         e = params.ellps.e;
+        spherical = params.ellps.spherical;
         //  Compute constants for the mlfn
         double t;

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

@@ -2 +2 @@
 package org.openstreetmap.josm.data.projection.proj;
 
-import static java.lang.Math.PI;
-import static java.lang.Math.atan;
-import static java.lang.Math.log;
-import static java.lang.Math.sinh;
-import static java.lang.Math.tan;
 import static org.openstreetmap.josm.tools.I18n.tr;
 
@@ -13 +8 @@
 
 /**
- * Mercator Projection.
+ * Mercator Cylindrical Projection. The parallels and the meridians are straight lines and
+ * cross at right angles; this projection thus produces rectangular charts. The scale is true
+ * along the equator (by default) or along two parallels equidistant of the equator (if a scale
+ * factor other than 1 is used). This projection is used to represent areas close to the equator.
+ * It is also often used for maritime navigation because all the straight lines on the chart are
+ * <em>loxodrome</em> lines, i.e. a ship following this line would keep a constant azimuth on its
+ * compass.
+ * <p>
+ * This implementation handles both the 1 and 2 stardard parallel cases.
+ * For 1 SP (EPSG code 9804), the line of contact is the equator.
+ * For 2 SP (EPSG code 9805) lines of contact are symmetrical
+ * about the equator.
+ * <p>
+ * This class has been derived from the implementation of the Geotools project;
+ * git 8cbf52d, org.geotools.referencing.operation.projection.CassiniSoldner
+ * at the time of migration.
+ * <p>
+ * <b>References:</b>
+ * <ul>
+ *   <li>John P. Snyder (Map Projections - A Working Manual,<br>
+ *       U.S. Geological Survey Professional Paper 1395, 1987)</li>
+ *   <li>"Coordinate Conversions and Transformations including Formulas",<br>
+ *       EPSG Guidence Note Number 7, Version 19.</li>
+ * </ul>
+ *
+ * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Mercator projection on MathWorld</A>
+ * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/mercator_1sp.html">"mercator_1sp" on RemoteSensing.org</A>
+ * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/mercator_2sp.html">"mercator_2sp" on RemoteSensing.org</A>
+ *
+ * @author André Gosselin
+ * @author Martin Desruisseaux (PMO, IRD)
+ * @author Rueben Schulz
+ * @author Simone Giannecchini
  */
-public class Mercator implements Proj {
+public class Mercator extends AbstractProj implements IScaleFactorProvider {
+    /**
+     * Maximum difference allowed when comparing real numbers.
+     */
+    private static final double EPSILON = 1E-6;
+
+    protected double scaleFactor;
 
     @Override
@@ -24 +57 @@
     @Override
     public String getProj4Id() {
-        return "josm:smerc"; // "merc" is ellipsoidal Mercator projection in PROJ.4
+        return "merc";
     }
 
     @Override
     public void initialize(ProjParameters params) throws ProjectionConfigurationException {
+        super.initialize(params);
+        scaleFactor = 1;
+        if (params.lat_ts != null) {
+            /*
+             * scaleFactor is not a parameter in the 2 SP case and is computed from
+             * the standard parallel.
+             */
+            double standardParallel = Math.toRadians(params.lat_ts);
+            if (spherical) {
+                scaleFactor *= Math.cos(standardParallel);
+            }  else {
+                scaleFactor *= msfn(Math.sin(standardParallel), Math.cos(standardParallel));
+            }
+        }
+        /*
+         * A correction that allows us to employs a latitude of origin that is not
+         * correspondent to the equator. See Snyder and al. for reference, page 47.
+         */
+        if (params.lat0 != null) {
+            final double lat0 = Math.toRadians(params.lat0);
+            final double sinPhi = Math.sin(lat0);
+            scaleFactor *= (Math.cos(lat0) / (Math.sqrt(1 - e2 * sinPhi * sinPhi)));
+        }
     }
 
     @Override
-    public double[] project(double lat_rad, double lon_rad) {
-        return new double[] {lon_rad, log(tan(PI/4 + lat_rad/2))};
+    public double[] project(double y, double x) {
+        if (Math.abs(y) > (Math.PI/2 - EPSILON)) {
+            return new double[] {0, 0}; // this is an error and should be handled somehow
+        }
+        if (spherical) {
+            y = Math.log(Math.tan(Math.PI/4 + 0.5*y));
+        } else {
+            y = -Math.log(tsfn(y, Math.sin(y)));
+        }
+        return new double[] {x, y};
     }
 
     @Override
-    public double[] invproject(double east, double north) {
-        return new double[] {atan(sinh(north)), east};
+    public double[] invproject(double x, double y) {
+        if (spherical) {
+            y = Math.PI/2 - 2.0*Math.atan(Math.exp(-y));
+        } else {
+            y = Math.exp(-y);
+            y = cphi2(y);
+        }
+        return new double[] {y, x};
     }
 
@@ -45 +115 @@
         return new Bounds(-89, -180, 89, 180, false);
     }
+
+    @Override
+    public double getScaleFactor() {
+        return scaleFactor;
+    }
 }

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

@@ -15 +15 @@
     public Double lat2;
 
-    // Polar Stereographic
+    // Polar Stereographic and Mercator
     public Double lat_ts;