Merge pull request #3222 from bartvde/geodesic-measure

Add geodesic option for measure

examples/measure.html

@@ -75,10 +75,11 @@
                 <option value="length">Length</option>
                 <option value="area">Area</option>
               </select>
+              <label class="checkbox"><input type="checkbox" id="geodesic"/>use geodesic measures</label>
           </form>
 
           <div id="docs">
-              <p><i>NOTE: Measure is done in simple way on projected plane. Earth
+              <p><i>NOTE: If use geodesic measures is not checked, measure is done in simple way on projected plane. Earth
               curvature is not taken into account</i></p>
             <p>See the <a href="measure.js" target="_blank">measure.js source</a> to see how this is done.</p>
           </div>

examples/measure.js

@@ -1,5 +1,6 @@
 goog.require('ol.Map');
 goog.require('ol.Overlay');
+goog.require('ol.Sphere');
 goog.require('ol.View');
 goog.require('ol.geom.LineString');
 goog.require('ol.geom.Polygon');
@@ -7,6 +8,7 @@ goog.require('ol.interaction');
 goog.require('ol.interaction.Draw');
 goog.require('ol.layer.Tile');
 goog.require('ol.layer.Vector');
+goog.require('ol.proj');
 goog.require('ol.source.MapQuest');
 goog.require('ol.source.Vector');
 goog.require('ol.style.Circle');
@@ -14,6 +16,9 @@ goog.require('ol.style.Fill');
 goog.require('ol.style.Stroke');
 goog.require('ol.style.Style');
 
+
+var wgs84Sphere = new ol.Sphere(6378137);
+
 var raster = new ol.layer.Tile({
   source: new ol.source.MapQuest({layer: 'sat'})
 });
@@ -135,6 +140,7 @@ var map = new ol.Map({
 map.on('pointermove', pointerMoveHandler);
 
 var typeSelect = document.getElementById('type');
+var geodesicCheckbox = document.getElementById('geodesic');
 
 var draw; // global so we can remove it later
 function addInteraction() {
@@ -238,7 +244,19 @@ typeSelect.onchange = function(e) {
  * @return {string}
  */
 var formatLength = function(line) {
-  var length = Math.round(line.getLength() * 100) / 100;
+  var length;
+  if (geodesicCheckbox.checked) {
+    var coordinates = line.getCoordinates();
+    length = 0;
+    var sourceProj = map.getView().getProjection();
+    for (var i = 0, ii = coordinates.length - 1; i < ii; ++i) {
+      var c1 = ol.proj.transform(coordinates[i], sourceProj, 'EPSG:4326');
+      var c2 = ol.proj.transform(coordinates[i + 1], sourceProj, 'EPSG:4326');
+      length += wgs84Sphere.haversineDistance(c1, c2);
+    }
+  } else {
+    length = Math.round(line.getLength() * 100) / 100;
+  }
   var output;
   if (length > 100) {
     output = (Math.round(length / 1000 * 100) / 100) +
@@ -257,7 +275,16 @@ var formatLength = function(line) {
  * @return {string}
  */
 var formatArea = function(polygon) {
-  var area = polygon.getArea();
+  var area;
+  if (geodesicCheckbox.checked) {
+    var sourceProj = map.getView().getProjection();
+    var geom = /** @type {ol.geom.Polygon} */(polygon.clone().transform(
+        sourceProj, 'EPSG:4326'));
+    var coordinates = geom.getLinearRing(0).getCoordinates();
+    area = Math.abs(wgs84Sphere.geodesicArea(coordinates));
+  } else {
+    area = polygon.getArea();
+  }
   var output;
   if (area > 10000) {
     output = (Math.round(area / 1000000 * 100) / 100) +

src/ol/sphere/sphere.js

@@ -57,6 +57,36 @@ ol.Sphere.prototype.cosineDistance = function(c1, c2) {
 };
 
 
+/**
+ * Returns the geodesic area for a list of coordinates.
+ *
+ * [Reference](http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409)
+ * Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
+ * Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
+ * Laboratory, Pasadena, CA, June 2007
+ *
+ * @param {Array.<ol.Coordinate>} coordinates List of coordinates of a linear
+ * ring. If the ring is oriented clockwise, the area will be positive,
+ * otherwise it will be negative.
+ * @return {number} Area.
+ * @api
+ */
+ol.Sphere.prototype.geodesicArea = function(coordinates) {
+  var area = 0, len = coordinates.length;
+  var x1 = coordinates[len - 1][0];
+  var y1 = coordinates[len - 1][1];
+  for (var i = 0; i < len; i++) {
+    var x2 = coordinates[i][0], y2 = coordinates[i][1];
+    area += goog.math.toRadians(x2 - x1) *
+        (2 + Math.sin(goog.math.toRadians(y1)) +
+        Math.sin(goog.math.toRadians(y2)));
+    x1 = x2;
+    y1 = y2;
+  }
+  return area * this.radius * this.radius / 2.0;
+};
+
+
 /**
  * Returns the distance of c3 from the great circle path defined by c1 and c2.
  *
@@ -110,6 +140,7 @@ ol.Sphere.prototype.finalBearing = function(c1, c2) {
  * @param {ol.Coordinate} c1 Coordinate 1.
  * @param {ol.Coordinate} c2 Coordinate 2.
  * @return {number} Haversine distance.
+ * @api
  */
 ol.Sphere.prototype.haversineDistance = function(c1, c2) {
   var lat1 = goog.math.toRadians(c1[1]);

test/spec/ol/sphere/sphere.test.js

@@ -264,8 +264,30 @@ describe('ol.Sphere', function() {
 
   });
 
+  describe('Vincenty area', function() {
+    var geometry;
+    before(function(done) {
+      afterLoadText('spec/ol/format/wkt/illinois.wkt', function(wkt) {
+        try {
+          var format = new ol.format.WKT();
+          geometry = format.readGeometry(wkt);
+        } catch (e) {
+          done(e);
+        }
+        done();
+      });
+    });
+
+    it('results match the expected area of Ilinois', function() {
+      var coords = geometry.getPolygon(0).getLinearRing(0).getCoordinates();
+      expect(ol.sphere.WGS84.geodesicArea(coords)).to.equal(145978332359.37125);
+    });
+  });
+
 });
 
 
 goog.require('goog.math');
 goog.require('ol.Sphere');
+goog.require('ol.sphere.WGS84');
+goog.require('ol.format.WKT');