Ticket #11516: _11516___Sinusoidal_projectionv2.patch
File _11516___Sinusoidal_projectionv2.patch, 9.0 KB (added by , 3 years ago) 


data_nodist/projection/josmepsg
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/josmepsg b/data_nodist/projection/josmepsg index 4ac4385..486a52d 100644
a b 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 <> 418 418 # NAD83(2011) / Oregon Columbia River West zone (m) 419 419 <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 40 40 import org.openstreetmap.josm.data.projection.proj.PolarStereographic; 41 41 import org.openstreetmap.josm.data.projection.proj.Proj; 42 42 import org.openstreetmap.josm.data.projection.proj.ProjFactory; 43 import org.openstreetmap.josm.data.projection.proj.Sinusoidal; 43 44 import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator; 44 45 import org.openstreetmap.josm.data.projection.proj.TransverseMercator; 45 46 import org.openstreetmap.josm.gui.preferences.projection.ProjectionChoice; … … public ProjectionDefinition(String code, String name, String definition) { 94 95 registerBaseProjection("merc", Mercator.class, "core"); 95 96 registerBaseProjection("omerc", ObliqueMercator.class, "core"); 96 97 registerBaseProjection("somerc", SwissObliqueMercator.class, "core"); 98 registerBaseProjection("sinu", Sinusoidal.class, "core"); 97 99 registerBaseProjection("stere", PolarStereographic.class, "core"); 98 100 registerBaseProjection("sterea", DoubleStereographic.class, "core"); 99 101 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. 2 package org.openstreetmap.josm.data.projection.proj; 3 4 import static java.lang.Math.abs; 5 import static java.lang.Math.cos; 6 import static java.lang.Math.sin; 7 import static java.lang.Math.sqrt; 8 import static org.openstreetmap.josm.tools.I18n.tr; 9 10 import org.openstreetmap.josm.data.Bounds; 11 12 /** 13 * Sinusoidal projection (aka. Sanson–Flamsteed, Mercator equalarea 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 */ 18 public 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) < 1e10) { 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 30 30 import org.openstreetmap.josm.data.osm.Relation; 31 31 import org.openstreetmap.josm.data.osm.RelationMember; 32 32 import org.openstreetmap.josm.data.osm.Way; 33 import org.openstreetmap.josm.data.projection.Projection; 34 import org.openstreetmap.josm.data.projection.Projections; 33 35 34 36 /** 35 37 * Some tools for geometry related tasks. … … public static boolean nodeInsidePolygon(Node point, List<Node> polygonNodes) { 629 631 630 632 /** 631 633 * 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, when635 * one unit in projected coordinates corresponds to one meter.636 * This is true for most projections, but not for WGS84 and637 * Mercator (EPSG:3857).638 634 * 639 635 * @param way Way to measure, should be closed (first node is the same as last node) 640 636 * @return area of the closed way. 641 637 */ 642 638 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(1a)); 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(1a)); 687 return 6367000 * c; 639 return getAreaAndPerimeter(way.getNodes(), Projections.getProjectionByCode("EPSG:54008")).getArea(); 688 640 } 689 641 690 642 /** … … public double getPerimeter() { 985 937 * @return area and perimeter 986 938 */ 987 939 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"); 988 952 double area = 0; 989 953 double perimeter = 0; 990 954 if (!nodes.isEmpty()) { 991 955 boolean closed = nodes.get(0) == nodes.get(nodes.size()  1); 992 956 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()); 994 958 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()); 996 961 area += p1.east() * p2.north()  p2.east() * p1.north(); 997 962 perimeter += p1.distance(p2); 998 963 p1 = p2;