Ticket #11516: _11516___Sinusoidal_projection-v2.patch

data_nodist/projection/josm-epsg

@@ -417 +417 @@
 <6811> +proj=omerc +lat_0=45.91666666666666 +lonc=-123 +alpha=295 +k=1 +x_0=7000000.00000248 +y_0=-2999999.999988 +no_uoff +gamma=295 +ellps=GRS80 +units=ft +no_defs +bounds=-125,45,-121,47 <>
 # NAD83(2011) / Oregon Columbia River West zone (m)
 <6810> +proj=omerc +lat_0=45.91666666666666 +lonc=-123 +alpha=295 +k=1 +x_0=7000000 +y_0=-3000000 +no_uoff +gamma=295 +ellps=GRS80 +units=m +no_defs +bounds=-125,45,-121,47 <>
+##
+## Following entries use Sinusoidal projection
+##
+# ESRI:54008 / World Sinusoidal
+<54008> +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs <>

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

@@ -40 +40 @@
 import org.openstreetmap.josm.data.projection.proj.PolarStereographic;
 import org.openstreetmap.josm.data.projection.proj.Proj;
 import org.openstreetmap.josm.data.projection.proj.ProjFactory;
+import org.openstreetmap.josm.data.projection.proj.Sinusoidal;
 import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator;
 import org.openstreetmap.josm.data.projection.proj.TransverseMercator;
 import org.openstreetmap.josm.gui.preferences.projection.ProjectionChoice;
@@ -94 +95 @@
         registerBaseProjection("merc", Mercator.class, "core");
         registerBaseProjection("omerc", ObliqueMercator.class, "core");
         registerBaseProjection("somerc", SwissObliqueMercator.class, "core");
+        registerBaseProjection("sinu", Sinusoidal.class, "core");
         registerBaseProjection("stere", PolarStereographic.class, "core");
         registerBaseProjection("sterea", DoubleStereographic.class, "core");
         registerBaseProjection("tmerc", TransverseMercator.class, "core");

new file src/org/openstreetmap/josm/data/projection/proj/Sinusoidal.java

@@ +1 @@
+// License: GPL. For details, see LICENSE file.
+package org.openstreetmap.josm.data.projection.proj;
+
+import static java.lang.Math.abs;
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import static java.lang.Math.sqrt;
+import static org.openstreetmap.josm.tools.I18n.tr;
+
+import org.openstreetmap.josm.data.Bounds;
+
+/**
+ * Sinusoidal projection (aka. Sanson–Flamsteed, Mercator equal-area projection)
+ * <p/>
+ * This class has been derived from the implementation of the <a href="https://github.com/geotools/geotools">Geotools</a> project;
+ * git 577dd2d, org.geotools.referencing.operation.projection.Sinusoidal at the time of migration.
+ */
+public class Sinusoidal extends AbstractProj {
+
+    @Override
+    public String getName() {
+        return tr("Sinusoidal");
+    }
+
+    @Override
+    public String getProj4Id() {
+        return "sinu";
+    }
+
+    @Override
+    public double[] project(final double phi, final double lambda) {
+        if (spherical) {
+            return new double[]{lambda * cos(phi), phi};
+        } else {
+            final double s = sin(phi);
+            return new double[]{lambda * cos(phi) / sqrt(1. - e2 * s * s), mlfn(phi, s, cos(phi))};
+        }
+    }
+
+    @Override
+    public double[] invproject(final double east, final double north) {
+        if (spherical) {
+            return new double[]{east / cos(north), north};
+        } else {
+            final double phi = inv_mlfn(north);
+            double s = abs(phi);
+            final double lambda;
+            if (abs(s - Math.PI / 2) < 1e-10) {
+                lambda = 0.;
+            } else if (s < Math.PI / 2) {
+                s = sin(phi);
+                lambda = (east * sqrt(1. - e2 * s * s) / cos(phi)) % Math.PI;
+            } else {
+                return new double[]{0., 0.}; // this is an error and should be handled somehow
+            }
+            return new double[]{phi, lambda};
+        }
+    }
+
+    @Override
+    public Bounds getAlgorithmBounds() {
+        return new Bounds(-90, -180, 90, 180, false);
+    }
+}

src/org/openstreetmap/josm/tools/Geometry.java

@@ -30 +30 @@
 import org.openstreetmap.josm.data.osm.Relation;
 import org.openstreetmap.josm.data.osm.RelationMember;
 import org.openstreetmap.josm.data.osm.Way;
+import org.openstreetmap.josm.data.projection.Projection;
+import org.openstreetmap.josm.data.projection.Projections;
 
 /**
  * Some tools for geometry related tasks.
@@ -629 +631 @@
 
     /**
      * Returns area of a closed way in square meters.
-     * (approximate(?), but should be OK for small areas)
-     *
-     * Relies on the current projection: Works correctly, when
-     * one unit in projected coordinates corresponds to one meter.
-     * This is true for most projections, but not for WGS84 and
-     * Mercator (EPSG:3857).
      *
      * @param way Way to measure, should be closed (first node is the same as last node)
      * @return area of the closed way.
      */
     public static double closedWayArea(Way way) {
-
-        //http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
-        double area = 0;
-        Node lastN = null;
-        for (Node n : way.getNodes()) {
-            if (lastN != null) {
-                area += (calcX(n) * calcY(lastN)) - (calcY(n) * calcX(lastN));
-            }
-            lastN = n;
-        }
-        return Math.abs(area/2);
-    }
-
-    protected static double calcX(Node p1) {
-        double lat1, lon1, lat2, lon2;
-        double dlon, dlat;
-
-        lat1 = p1.getCoor().lat() * Math.PI / 180.0;
-        lon1 = p1.getCoor().lon() * Math.PI / 180.0;
-        lat2 = lat1;
-        lon2 = 0;
-
-        dlon = lon2 - lon1;
-        dlat = lat2 - lat1;
-
-        double a = Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2);
-        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
-        return 6367000 * c;
-    }
-
-    protected static double calcY(Node p1) {
-        double lat1, lon1, lat2, lon2;
-        double dlon, dlat;
-
-        lat1 = p1.getCoor().lat() * Math.PI / 180.0;
-        lon1 = p1.getCoor().lon() * Math.PI / 180.0;
-        lat2 = 0;
-        lon2 = lon1;
-
-        dlon = lon2 - lon1;
-        dlat = lat2 - lat1;
-
-        double a = Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2);
-        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
-        return 6367000 * c;
+        return getAreaAndPerimeter(way.getNodes(), Projections.getProjectionByCode("EPSG:54008")).getArea();
     }
 
     /**
@@ -985 +937 @@
      * @return area and perimeter
      */
     public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes) {
+        return getAreaAndPerimeter(nodes, null);
+    }
+
+    /**
+     * Calculate area and perimeter length of a polygon in the given projection.
+     *
+     * @param nodes the list of nodes representing the polygon
+     * @param projection the projection to use for the calculation, {@code null} defaults to {@link Main#getProjection()}
+     * @return area and perimeter
+     */
+    public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes, Projection projection) {
+        CheckParameterUtil.ensureParameterNotNull(nodes, "nodes");
         double area = 0;
         double perimeter = 0;
         if (!nodes.isEmpty()) {
             boolean closed = nodes.get(0) == nodes.get(nodes.size() - 1);
             int numSegments = closed ? nodes.size() - 1 : nodes.size();
-            EastNorth p1 = nodes.get(0).getEastNorth();
+            EastNorth p1 = projection == null ? nodes.get(0).getEastNorth() : projection.latlon2eastNorth(nodes.get(0).getCoor());
             for (int i = 1; i <= numSegments; i++) {
-                EastNorth p2 = nodes.get(i == numSegments ? 0 : i).getEastNorth();
+                final Node node = nodes.get(i == numSegments ? 0 : i);
+                final EastNorth p2 = projection == null ? node.getEastNorth() : projection.latlon2eastNorth(node.getCoor());
                 area += p1.east() * p2.north() - p2.east() * p1.north();
                 perimeter += p1.distance(p2);
                 p1 = p2;