Function for getting great circle lengths

This commit is contained in:
Tim Schaub
2017-08-10 17:45:22 -06:00
parent 445c157ee3
commit 92c62e5432
5 changed files with 212 additions and 18 deletions

View File

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

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 * @private
* @type {ol.Sphere} * @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) { if (ol.ENABLE_PROJ4JS) {
@@ -90,9 +90,9 @@ ol.proj.getPointResolution = function(projection, resolution, point, opt_units)
point[0], point[1] + resolution / 2 point[0], point[1] + resolution / 2
]; ];
vertices = toEPSG4326(vertices, vertices, 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)); 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)); vertices.slice(4, 6), vertices.slice(6, 8));
pointResolution = (width + height) / 2; pointResolution = (width + height) / 2;
var metersPerUnit = opt_units ? var metersPerUnit = opt_units ?

View File

@@ -8,6 +8,7 @@
goog.provide('ol.Sphere'); goog.provide('ol.Sphere');
goog.require('ol.math'); goog.require('ol.math');
goog.require('ol.geom.GeometryType');
/** /**
@@ -75,14 +76,7 @@ ol.Sphere.prototype.geodesicArea = function(coordinates) {
* @api * @api
*/ */
ol.Sphere.prototype.haversineDistance = function(c1, c2) { ol.Sphere.prototype.haversineDistance = function(c1, c2) {
var lat1 = ol.math.toRadians(c1[1]); return ol.Sphere.getDistance_(c1, c2, this.radius);
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));
}; };
@@ -107,3 +101,108 @@ ol.Sphere.prototype.offset = function(c1, distance, bearing) {
Math.cos(dByR) - Math.sin(lat1) * Math.sin(lat)); Math.cos(dByR) - Math.sin(lat1) * Math.sin(lat));
return [ol.math.toDegrees(lon), ol.math.toDegrees(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.SphereLengthOptions} 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));
};

View File

@@ -221,21 +221,21 @@ describe('ol.proj', function() {
}); });
it('returns the correct point resolution for EPSG:4326 with custom units', 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'); 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'); 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() { it('returns the correct point resolution for EPSG:3857', function() {
var pointResolution = ol.proj.getPointResolution('EPSG:3857', 1, [0, 0]); 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])); 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() { it('returns the correct point resolution for EPSG:3857 with custom units', function() {
var pointResolution = ol.proj.getPointResolution('EPSG:3857', 1, [0, 0], 'degrees'); 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'); pointResolution = ol.proj.getPointResolution('EPSG:4326', 1, ol.proj.fromLonLat([0, 52]), 'degrees');
expect (pointResolution).to.be(1); expect(pointResolution).to.be(1);
}); });
}); });

View File

@@ -108,3 +108,72 @@ 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);
});
});
});