Changeset 9951 in josm


Ignore:
Timestamp:
2016-03-07T23:29:23+01:00 (3 years ago)
Author:
simon04
Message:

see #11516 - Compute closed way area using Sinusoidal projection

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/org/openstreetmap/josm/tools/Geometry.java

    r9231 r9951  
    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/**
     
    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)
     
    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
     
    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;
     
    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);
  • trunk/test/unit/org/openstreetmap/josm/tools/GeometryTest.java

    r9878 r9951  
    7171            DataSet ds = OsmReader.parseDataSet(in, null);
    7272            Way closedWay = (Way) Utils.filter(ds.allPrimitives(), SearchCompiler.compile("landuse=forest")).iterator().next();
    73             Assert.assertEquals(5721923.660644531, Geometry.closedWayArea(closedWay), 1e-3);
     73            Assert.assertEquals(5760015.7353515625, Geometry.closedWayArea(closedWay), 1e-3);
    7474        }
    7575    }
Note: See TracChangeset for help on using the changeset viewer.