ashatora/openlayers
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Merge pull request #7105 from tschaub/sphere

Browse Source 
Functions for spherical calculations
This commit is contained in:
Tim Schaub
2017-08-10 20:47:38 -06:00
committed by GitHub
parent 35f79f401b 5f794ab562
commit f2c51cb39e
7 changed files with 339 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
examples/measure.html
View File

@@ -1,20 +1,24 @@
---
layout: example.html
title: Measure
shortdesc: Example of using the ol.interaction.Draw interaction to create a simple measuring application.
shortdesc: Using a draw interaction to measure lengths and areas.
docs: >
<p><i>NOTE: By default, length and area are calculated using the projected coordinates. This is not accurate for projections like Mercator where the projected meters do not correspond to meters on the ground. To get a standarized measurement across all projections, use the geodesic measures.</i></p>
<p>The <code>ol.Sphere.getLength()</code> and <code>ol.Sphere.getArea()</code>
functions calculate spherical lengths and areas for geometries. Lengths are
calculated by assuming great circle segments between geometry coordinates.
Areas are calculated as if edges of polygons were great circle segments.</p>
<p>Note that the <code>geometry.getLength()</code> and <code>geometry.getArea()</code>
methods return measures of projected (planar) geometries. These can be very
different than on-the-ground measures in certain situations — in northern
and southern latitudes using Web Mercator for example. For better results,
use the functions on <code>ol.Sphere</code>.</p>
tags: "draw, edit, measure, vector"
---
<div id="map" class="map"></div>
<form class="form-inline">
<label>Measurement type &nbsp;</label>
<select id="type">
<option value="length">Length (LineString)</option>
<option value="area">Area (Polygon)</option>
</select>
<label class="checkbox">
<input type="checkbox" id="geodesic">
use geodesic measures
</label>
<select id="type">
<option value="length">Length (LineString)</option>
<option value="area">Area (Polygon)</option>
</select>
</form>

29
examples/measure.js
View File

@@ -8,7 +8,6 @@ goog.require('ol.geom.Polygon');
goog.require('ol.interaction.Draw');
goog.require('ol.layer.Tile');
goog.require('ol.layer.Vector');
goog.require('ol.proj');
goog.require('ol.source.OSM');
goog.require('ol.source.Vector');
goog.require('ol.style.Circle');
@@ -17,8 +16,6 @@ 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.OSM()
});
@@ -137,7 +134,6 @@ map.getViewport().addEventListener('mouseout', function() {
});
var typeSelect = document.getElementById('type');
var geodesicCheckbox = document.getElementById('geodesic');
var draw; // global so we can remove it later
@@ -148,19 +144,7 @@ var draw; // global so we can remove it later
* @return {string} The formatted length.
*/
var formatLength = function(line) {
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 length = ol.Sphere.getLength(line);
var output;
if (length > 100) {
output = (Math.round(length / 1000 * 100) / 100) +
@@ -179,16 +163,7 @@ var formatLength = function(line) {
* @return {string} Formatted area.
*/
var formatArea = function(polygon) {
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 area = ol.Sphere.getArea(polygon);
var output;
if (area > 10000) {
output = (Math.round(area / 1000000 * 100) / 100) +

27
externs/olx.js
View File

@@ -436,6 +436,33 @@ olx.MapOptions.prototype.target;
olx.MapOptions.prototype.view;
/**
* Object literal with options for the {@link ol.Sphere.getLength} or
* {@link ol.Sphere.getArea} functions.
* @typedef {{projection: (ol.ProjectionLike|undefined),
* radius: (number|undefined)}}
*/
olx.SphereMetricOptions;
/**
* Projection of the geometry. By default, the geometry is assumed to be in
* EPSG:3857 (Web Mercator).
* @type {(ol.ProjectionLike|undefined)}
* @api
*/
olx.SphereMetricOptions.prototype.projection;
/**
* Sphere radius. By default, the radius of the earth is used (Clarke 1866
* Authalic Sphere).
* @type {(number|undefined)}
* @api
*/
olx.SphereMetricOptions.prototype.radius;
/**
* Object literal with options for the {@link ol.Map#forEachFeatureAtPixel} and
* {@link ol.Map#hasFeatureAtPixel} methods.

8
src/ol/proj.js
View File

@@ -22,11 +22,11 @@ ol.proj.METERS_PER_UNIT = ol.proj.Units.METERS_PER_UNIT;
/**
* A place to store the radius of the Clarke 1866 Authalic Sphere.
* A place to store the mean radius of the Earth.
* @private
* @type {ol.Sphere}
*/
ol.proj.AUTHALIC_SPHERE_ = new ol.Sphere(6370997);
ol.proj.SPHERE_ = new ol.Sphere(ol.Sphere.DEFAULT_RADIUS);
if (ol.ENABLE_PROJ4JS) {
@@ -90,9 +90,9 @@ ol.proj.getPointResolution = function(projection, resolution, point, opt_units)
point[0], point[1] + resolution / 2
];
vertices = toEPSG4326(vertices, vertices, 2);
var width = ol.proj.AUTHALIC_SPHERE_.haversineDistance(
var width = ol.proj.SPHERE_.haversineDistance(
vertices.slice(0, 2), vertices.slice(2, 4));
var height = ol.proj.AUTHALIC_SPHERE_.haversineDistance(
var height = ol.proj.SPHERE_.haversineDistance(
vertices.slice(4, 6), vertices.slice(6, 8));
pointResolution = (width + height) / 2;
var metersPerUnit = opt_units ?

217
src/ol/sphere.js
View File

@@ -8,6 +8,7 @@
goog.provide('ol.Sphere');
goog.require('ol.math');
goog.require('ol.geom.GeometryType');
/**
@@ -51,18 +52,7 @@ ol.Sphere = function(radius) {
* @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 += ol.math.toRadians(x2 - x1) *
(2 + Math.sin(ol.math.toRadians(y1)) +
Math.sin(ol.math.toRadians(y2)));
x1 = x2;
y1 = y2;
}
return area * this.radius * this.radius / 2.0;
return ol.Sphere.getArea_(coordinates, this.radius);
};
@@ -75,14 +65,7 @@ ol.Sphere.prototype.geodesicArea = function(coordinates) {
* @api
*/
ol.Sphere.prototype.haversineDistance = function(c1, c2) {
var lat1 = ol.math.toRadians(c1[1]);
var lat2 = ol.math.toRadians(c2[1]);
var deltaLatBy2 = (lat2 - lat1) / 2;
var deltaLonBy2 = ol.math.toRadians(c2[0] - c1[0]) / 2;
var a = Math.sin(deltaLatBy2) * Math.sin(deltaLatBy2) +
Math.sin(deltaLonBy2) * Math.sin(deltaLonBy2) *
Math.cos(lat1) * Math.cos(lat2);
return 2 * this.radius * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return ol.Sphere.getDistance_(c1, c2, this.radius);
};
@@ -107,3 +90,197 @@ ol.Sphere.prototype.offset = function(c1, distance, bearing) {
Math.cos(dByR) - Math.sin(lat1) * Math.sin(lat));
return [ol.math.toDegrees(lon), ol.math.toDegrees(lat)];
};
/**
* The mean Earth radius (1/3 * (2a + b)) for the WGS84 ellipsoid.
* https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
* @type {number}
*/
ol.Sphere.DEFAULT_RADIUS = 6371008.8;
/**
* Get the spherical length of a geometry. This length is the sum of the
* great circle distances between coordinates. For polygons, the length is
* the sum of all rings. For points, the length is zero. For multi-part
* geometries, the length is the sum of the length of each part.
* @param {ol.geom.Geometry} geometry A geometry.
* @param {olx.SphereMetricOptions=} opt_options Options for the length
* calculation. By default, geometries are assumed to be in 'EPSG:3857'.
* You can change this by providing a `projection` option.
* @return {number} The spherical length (in meters).
*/
ol.Sphere.getLength = function(geometry, opt_options) {
var options = opt_options || {};
var radius = options.radius || ol.Sphere.DEFAULT_RADIUS;
var projection = options.projection || 'EPSG:3857';
geometry = geometry.clone().transform(projection, 'EPSG:4326');
var type = geometry.getType();
var length = 0;
var coordinates, coords, i, ii, j, jj;
switch (type) {
case ol.geom.GeometryType.POINT:
case ol.geom.GeometryType.MULTI_POINT: {
break;
}
case ol.geom.GeometryType.LINE_STRING:
case ol.geom.GeometryType.LINEAR_RING: {
coordinates = /** @type {ol.geom.SimpleGeometry} */ (geometry).getCoordinates();
length = ol.Sphere.getLength_(coordinates, radius);
break;
}
case ol.geom.GeometryType.MULTI_LINE_STRING:
case ol.geom.GeometryType.POLYGON: {
coordinates = /** @type {ol.geom.SimpleGeometry} */ (geometry).getCoordinates();
for (i = 0, ii = coordinates.length; i < ii; ++i) {
length += ol.Sphere.getLength_(coordinates[i], radius);
}
break;
}
case ol.geom.GeometryType.MULTI_POLYGON: {
coordinates = /** @type {ol.geom.SimpleGeometry} */ (geometry).getCoordinates();
for (i = 0, ii = coordinates.length; i < ii; ++i) {
coords = coordinates[i];
for (j = 0, jj = coords.length; j < jj; ++j) {
length += ol.Sphere.getLength_(coords[j], radius);
}
}
break;
}
case ol.geom.GeometryType.GEOMETRY_COLLECTION: {
var geometries = /** @type {ol.geom.GeometryCollection} */ (geometry).getGeometries();
for (i = 0, ii = geometries.length; i < ii; ++i) {
length += ol.Sphere.getLength(geometries[i], opt_options);
}
break;
}
default: {
throw new Error('Unsupported geometry type: ' + type);
}
}
return length;
};
/**
* Get the cumulative great circle length of linestring coordinates (geographic).
* @param {Array} coordinates Linestring coordinates.
* @param {number} radius The sphere radius to use.
* @return {number} The length (in meters).
*/
ol.Sphere.getLength_ = function(coordinates, radius) {
var length = 0;
for (var i = 0, ii = coordinates.length; i < ii - 1; ++i) {
length += ol.Sphere.getDistance_(coordinates[i], coordinates[i + 1], radius);
}
return length;
};
/**
* Get the great circle distance between two geographic coordinates.
* @param {Array} c1 Starting coordinate.
* @param {Array} c2 Ending coordinate.
* @param {number} radius The sphere radius to use.
* @return {number} The great circle distance between the points (in meters).
*/
ol.Sphere.getDistance_ = function(c1, c2, radius) {
var lat1 = ol.math.toRadians(c1[1]);
var lat2 = ol.math.toRadians(c2[1]);
var deltaLatBy2 = (lat2 - lat1) / 2;
var deltaLonBy2 = ol.math.toRadians(c2[0] - c1[0]) / 2;
var a = Math.sin(deltaLatBy2) * Math.sin(deltaLatBy2) +
Math.sin(deltaLonBy2) * Math.sin(deltaLonBy2) *
Math.cos(lat1) * Math.cos(lat2);
return 2 * radius * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
};
/**
* Get the spherical area of a geometry. This is the area (in meters) assuming
* that polygon edges are segments of great circles on a sphere.
* @param {ol.geom.Geometry} geometry A geometry.
* @param {olx.SphereMetricOptions=} opt_options Options for the area
* calculation. By default, geometries are assumed to be in 'EPSG:3857'.
* You can change this by providing a `projection` option.
* @return {number} The spherical area (in square meters).
*/
ol.Sphere.getArea = function(geometry, opt_options) {
var options = opt_options || {};
var radius = options.radius || ol.Sphere.DEFAULT_RADIUS;
var projection = options.projection || 'EPSG:3857';
geometry = geometry.clone().transform(projection, 'EPSG:4326');
var type = geometry.getType();
var area = 0;
var coordinates, coords, i, ii, j, jj;
switch (type) {
case ol.geom.GeometryType.POINT:
case ol.geom.GeometryType.MULTI_POINT:
case ol.geom.GeometryType.LINE_STRING:
case ol.geom.GeometryType.MULTI_LINE_STRING:
case ol.geom.GeometryType.LINEAR_RING: {
break;
}
case ol.geom.GeometryType.POLYGON: {
coordinates = /** @type {ol.geom.Polygon} */ (geometry).getCoordinates();
area = Math.abs(ol.Sphere.getArea_(coordinates[0], radius));
for (i = 1, ii = coordinates.length; i < ii; ++i) {
area -= Math.abs(ol.Sphere.getArea_(coordinates[i], radius));
}
break;
}
case ol.geom.GeometryType.MULTI_POLYGON: {
coordinates = /** @type {ol.geom.SimpleGeometry} */ (geometry).getCoordinates();
for (i = 0, ii = coordinates.length; i < ii; ++i) {
coords = coordinates[i];
area += Math.abs(ol.Sphere.getArea_(coords[0], radius));
for (j = 1, jj = coords.length; j < jj; ++j) {
area -= Math.abs(ol.Sphere.getArea_(coords[j], radius));
}
}
break;
}
case ol.geom.GeometryType.GEOMETRY_COLLECTION: {
var geometries = /** @type {ol.geom.GeometryCollection} */ (geometry).getGeometries();
for (i = 0, ii = geometries.length; i < ii; ++i) {
area += ol.Sphere.getArea(geometries[i], opt_options);
}
break;
}
default: {
throw new Error('Unsupported geometry type: ' + type);
}
}
return area;
};
/**
* Returns the spherical area for a list of coordinates.
*
* [Reference](https://trs-new.jpl.nasa.gov/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.
* @param {number} radius The sphere radius.
* @return {number} Area (in square meters).
*/
ol.Sphere.getArea_ = function(coordinates, radius) {
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 += ol.math.toRadians(x2 - x1) *
(2 + Math.sin(ol.math.toRadians(y1)) +
Math.sin(ol.math.toRadians(y2)));
x1 = x2;
y1 = y2;
}
return area * radius * radius / 2.0;
};

12
test/spec/ol/proj/index.test.js
View File

@@ -221,21 +221,21 @@ describe('ol.proj', function() {
});
it('returns the correct point resolution for EPSG:4326 with custom units', function() {
var pointResolution = ol.proj.getPointResolution('EPSG:4326', 1, [0, 0], 'm');
expect (pointResolution).to.roughlyEqual(111194.874284, 1e-5);
expect(pointResolution).to.roughlyEqual(111195.0802335329, 1e-5);
pointResolution = ol.proj.getPointResolution('EPSG:4326', 1, [0, 52], 'm');
expect (pointResolution).to.roughlyEqual(89826.367538, 1e-5);
expect(pointResolution).to.roughlyEqual(89826.53390979706, 1e-5);
});
it('returns the correct point resolution for EPSG:3857', function() {
var pointResolution = ol.proj.getPointResolution('EPSG:3857', 1, [0, 0]);
expect (pointResolution).to.be(1);
expect(pointResolution).to.be(1);
pointResolution = ol.proj.getPointResolution('EPSG:3857', 1, ol.proj.fromLonLat([0, 52]));
expect (pointResolution).to.roughlyEqual(0.615661, 1e-5);
expect(pointResolution).to.roughlyEqual(0.615661, 1e-5);
});
it('returns the correct point resolution for EPSG:3857 with custom units', function() {
var pointResolution = ol.proj.getPointResolution('EPSG:3857', 1, [0, 0], 'degrees');
expect (pointResolution).to.be(1);
expect(pointResolution).to.be(1);
pointResolution = ol.proj.getPointResolution('EPSG:4326', 1, ol.proj.fromLonLat([0, 52]), 'degrees');
expect (pointResolution).to.be(1);
expect(pointResolution).to.be(1);
});
});

89
test/spec/ol/sphere/index.test.js → test/spec/ol/sphere.test.js
View File

@@ -108,3 +108,92 @@ describe('ol.Sphere', function() {
});
});
describe('ol.Sphere.getLength()', function() {
var cases = [{
geometry: new ol.geom.Point([0, 0]),
length: 0
}, {
geometry: new ol.geom.MultiPoint([[0, 0], [1, 1]]),
length: 0
}, {
geometry: new ol.geom.LineString([
[12801741.441226462, -3763310.627144653],
[14582853.293918837, -2511525.2348457114],
[15918687.18343812, -2875744.624352243],
[16697923.618991036, -4028802.0261344076]
]),
length: 4407939.124914191
}, {
geometry: new ol.geom.LineString([
[115, -32],
[131, -22],
[143, -25],
[150, -34]
]),
options: {projection: 'EPSG:4326'},
length: 4407939.124914191
}, {
geometry: new ol.geom.MultiLineString([
[
[115, -32],
[131, -22],
[143, -25],
[150, -34]
], [
[115, -32],
[131, -22],
[143, -25],
[150, -34]
]
]),
options: {projection: 'EPSG:4326'},
length: 2 * 4407939.124914191
}, {
geometry: new ol.geom.GeometryCollection([
new ol.geom.LineString([
[115, -32],
[131, -22],
[143, -25],
[150, -34]
]),
new ol.geom.LineString([
[115, -32],
[131, -22],
[143, -25],
[150, -34]
])
]),
options: {projection: 'EPSG:4326'},
length: 2 * 4407939.124914191
}];
cases.forEach(function(c, i) {
it('works for case ' + i, function() {
var c = cases[i];
var length = ol.Sphere.getLength(c.geometry, c.options);
expect(length).to.equal(c.length);
});
});
});
describe('ol.Sphere.getArea()', 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('calculates the area of Ilinois', function() {
var area = ol.Sphere.getArea(geometry, {projection: 'EPSG:4326'});
expect(area).to.equal(145652224192.4434);
});
});