Ticket #11516: _11516___Sinusoidal_projection-v2.patch

File _11516___Sinusoidal_projection-v2.patch, 9.0 KB (added by simon04, 3 years ago)
  • data_nodist/projection/josm-epsg

    commit c2cd233d2f709294012268f6d7298d7bb98142e7
    Author: Simon Legner <Simon.Legner@gmail.com>
    Date:   Sun Mar 6 21:54:02 2016 +0100
    
        11516
    
    diff --git a/data_nodist/projection/josm-epsg b/data_nodist/projection/josm-epsg
    index 4ac4385..486a52d 100644
    a b  
    417417<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 <>
    418418# NAD83(2011) / Oregon Columbia River West zone (m)
    419419<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 <>
     420##
     421## Following entries use Sinusoidal projection
     422##
     423# ESRI:54008 / World Sinusoidal
     424<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

    diff --git a/src/org/openstreetmap/josm/data/projection/Projections.java b/src/org/openstreetmap/josm/data/projection/Projections.java
    index a39c9a1..931533a 100644
    a b  
    4040import org.openstreetmap.josm.data.projection.proj.PolarStereographic;
    4141import org.openstreetmap.josm.data.projection.proj.Proj;
    4242import org.openstreetmap.josm.data.projection.proj.ProjFactory;
     43import org.openstreetmap.josm.data.projection.proj.Sinusoidal;
    4344import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator;
    4445import org.openstreetmap.josm.data.projection.proj.TransverseMercator;
    4546import org.openstreetmap.josm.gui.preferences.projection.ProjectionChoice;
    public ProjectionDefinition(String code, String name, String definition) { 
    9495        registerBaseProjection("merc", Mercator.class, "core");
    9596        registerBaseProjection("omerc", ObliqueMercator.class, "core");
    9697        registerBaseProjection("somerc", SwissObliqueMercator.class, "core");
     98        registerBaseProjection("sinu", Sinusoidal.class, "core");
    9799        registerBaseProjection("stere", PolarStereographic.class, "core");
    98100        registerBaseProjection("sterea", DoubleStereographic.class, "core");
    99101        registerBaseProjection("tmerc", TransverseMercator.class, "core");
  • new file src/org/openstreetmap/josm/data/projection/proj/Sinusoidal.java

    diff --git a/src/org/openstreetmap/josm/data/projection/proj/Sinusoidal.java b/src/org/openstreetmap/josm/data/projection/proj/Sinusoidal.java
    new file mode 100644
    index 0000000..c7b8570
    - +  
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection.proj;
     3
     4import static java.lang.Math.abs;
     5import static java.lang.Math.cos;
     6import static java.lang.Math.sin;
     7import static java.lang.Math.sqrt;
     8import static org.openstreetmap.josm.tools.I18n.tr;
     9
     10import org.openstreetmap.josm.data.Bounds;
     11
     12/**
     13 * Sinusoidal projection (aka. Sanson–Flamsteed, Mercator equal-area projection)
     14 * <p/>
     15 * This class has been derived from the implementation of the <a href="https://github.com/geotools/geotools">Geotools</a> project;
     16 * git 577dd2d, org.geotools.referencing.operation.projection.Sinusoidal at the time of migration.
     17 */
     18public class Sinusoidal extends AbstractProj {
     19
     20    @Override
     21    public String getName() {
     22        return tr("Sinusoidal");
     23    }
     24
     25    @Override
     26    public String getProj4Id() {
     27        return "sinu";
     28    }
     29
     30    @Override
     31    public double[] project(final double phi, final double lambda) {
     32        if (spherical) {
     33            return new double[]{lambda * cos(phi), phi};
     34        } else {
     35            final double s = sin(phi);
     36            return new double[]{lambda * cos(phi) / sqrt(1. - e2 * s * s), mlfn(phi, s, cos(phi))};
     37        }
     38    }
     39
     40    @Override
     41    public double[] invproject(final double east, final double north) {
     42        if (spherical) {
     43            return new double[]{east / cos(north), north};
     44        } else {
     45            final double phi = inv_mlfn(north);
     46            double s = abs(phi);
     47            final double lambda;
     48            if (abs(s - Math.PI / 2) < 1e-10) {
     49                lambda = 0.;
     50            } else if (s < Math.PI / 2) {
     51                s = sin(phi);
     52                lambda = (east * sqrt(1. - e2 * s * s) / cos(phi)) % Math.PI;
     53            } else {
     54                return new double[]{0., 0.}; // this is an error and should be handled somehow
     55            }
     56            return new double[]{phi, lambda};
     57        }
     58    }
     59
     60    @Override
     61    public Bounds getAlgorithmBounds() {
     62        return new Bounds(-90, -180, 90, 180, false);
     63    }
     64}
  • src/org/openstreetmap/josm/tools/Geometry.java

    diff --git a/src/org/openstreetmap/josm/tools/Geometry.java b/src/org/openstreetmap/josm/tools/Geometry.java
    index 0f62e1b..0a8a4ac 100644
    a b  
    3030import org.openstreetmap.josm.data.osm.Relation;
    3131import org.openstreetmap.josm.data.osm.RelationMember;
    3232import org.openstreetmap.josm.data.osm.Way;
     33import org.openstreetmap.josm.data.projection.Projection;
     34import org.openstreetmap.josm.data.projection.Projections;
    3335
    3436/**
    3537 * Some tools for geometry related tasks.
    public static boolean nodeInsidePolygon(Node point, List<Node> polygonNodes) { 
    629631
    630632    /**
    631633     * Returns area of a closed way in square meters.
    632      * (approximate(?), but should be OK for small areas)
    633      *
    634      * Relies on the current projection: Works correctly, when
    635      * one unit in projected coordinates corresponds to one meter.
    636      * This is true for most projections, but not for WGS84 and
    637      * Mercator (EPSG:3857).
    638634     *
    639635     * @param way Way to measure, should be closed (first node is the same as last node)
    640636     * @return area of the closed way.
    641637     */
    642638    public static double closedWayArea(Way way) {
    643 
    644         //http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
    645         double area = 0;
    646         Node lastN = null;
    647         for (Node n : way.getNodes()) {
    648             if (lastN != null) {
    649                 area += (calcX(n) * calcY(lastN)) - (calcY(n) * calcX(lastN));
    650             }
    651             lastN = n;
    652         }
    653         return Math.abs(area/2);
    654     }
    655 
    656     protected static double calcX(Node p1) {
    657         double lat1, lon1, lat2, lon2;
    658         double dlon, dlat;
    659 
    660         lat1 = p1.getCoor().lat() * Math.PI / 180.0;
    661         lon1 = p1.getCoor().lon() * Math.PI / 180.0;
    662         lat2 = lat1;
    663         lon2 = 0;
    664 
    665         dlon = lon2 - lon1;
    666         dlat = lat2 - lat1;
    667 
    668         double a = Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2);
    669         double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    670         return 6367000 * c;
    671     }
    672 
    673     protected static double calcY(Node p1) {
    674         double lat1, lon1, lat2, lon2;
    675         double dlon, dlat;
    676 
    677         lat1 = p1.getCoor().lat() * Math.PI / 180.0;
    678         lon1 = p1.getCoor().lon() * Math.PI / 180.0;
    679         lat2 = 0;
    680         lon2 = lon1;
    681 
    682         dlon = lon2 - lon1;
    683         dlat = lat2 - lat1;
    684 
    685         double a = Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2);
    686         double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    687         return 6367000 * c;
     639        return getAreaAndPerimeter(way.getNodes(), Projections.getProjectionByCode("EPSG:54008")).getArea();
    688640    }
    689641
    690642    /**
    public double getPerimeter() { 
    985937     * @return area and perimeter
    986938     */
    987939    public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes) {
     940        return getAreaAndPerimeter(nodes, null);
     941    }
     942
     943    /**
     944     * Calculate area and perimeter length of a polygon in the given projection.
     945     *
     946     * @param nodes the list of nodes representing the polygon
     947     * @param projection the projection to use for the calculation, {@code null} defaults to {@link Main#getProjection()}
     948     * @return area and perimeter
     949     */
     950    public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes, Projection projection) {
     951        CheckParameterUtil.ensureParameterNotNull(nodes, "nodes");
    988952        double area = 0;
    989953        double perimeter = 0;
    990954        if (!nodes.isEmpty()) {
    991955            boolean closed = nodes.get(0) == nodes.get(nodes.size() - 1);
    992956            int numSegments = closed ? nodes.size() - 1 : nodes.size();
    993             EastNorth p1 = nodes.get(0).getEastNorth();
     957            EastNorth p1 = projection == null ? nodes.get(0).getEastNorth() : projection.latlon2eastNorth(nodes.get(0).getCoor());
    994958            for (int i = 1; i <= numSegments; i++) {
    995                 EastNorth p2 = nodes.get(i == numSegments ? 0 : i).getEastNorth();
     959                final Node node = nodes.get(i == numSegments ? 0 : i);
     960                final EastNorth p2 = projection == null ? node.getEastNorth() : projection.latlon2eastNorth(node.getCoor());
    996961                area += p1.east() * p2.north() - p2.east() * p1.north();
    997962                perimeter += p1.distance(p2);
    998963                p1 = p2;