diff --git a/src/ol/geom/flat/geodesicflatgeom.js b/src/ol/geom/flat/geodesicflatgeom.js new file mode 100644 index 0000000000..53d3aa8925 --- /dev/null +++ b/src/ol/geom/flat/geodesicflatgeom.js @@ -0,0 +1,180 @@ +goog.provide('ol.geom.flat.geodesic'); + +goog.require('goog.asserts'); +goog.require('goog.math'); +goog.require('goog.object'); +goog.require('ol.TransformFunction'); +goog.require('ol.math'); +goog.require('ol.proj'); + + +/** + * @private + * @param {function(number): ol.Coordinate} interpolate Interpolate function. + * @param {ol.TransformFunction} transform Transform from longitude/latitude to + * projected coordinates. + * @param {number} squaredTolerance Squared tolerance. + * @return {Array.} Flat coordinates. + */ +ol.geom.flat.geodesic.line_ = + function(interpolate, transform, squaredTolerance) { + // FIXME reduce garbage generation + // FIXME optimize stack operations + + /** @type {Array.} */ + var flatCoordinates = []; + + var geoA = interpolate(0); + var geoB = interpolate(1); + + var a = transform(geoA); + var b = transform(geoB); + + /** @type {Array.} */ + var geoStack = [geoB, geoA]; + /** @type {Array.} */ + var stack = [b, a]; + /** @type {Array.} */ + var fractionStack = [1, 0]; + + /** @type {Object.} */ + var fractions = {}; + + var maxIterations = 1e5; + var geoM, m, fracA, fracB, fracM, key; + + while (--maxIterations > 0 && fractionStack.length > 0) { + // Pop the a coordinate off the stack + fracA = fractionStack.pop(); + geoA = geoStack.pop(); + a = stack.pop(); + // Add the a coordinate if it has not been added yet + key = fracA.toString(); + if (!goog.object.containsKey(fractions, key)) { + flatCoordinates.push(a[0], a[1]); + goog.object.set(fractions, key, true); + } + // Pop the b coordinate off the stack + fracB = fractionStack.pop(); + geoB = geoStack.pop(); + b = stack.pop(); + // Find the m point between the a and b coordinates + fracM = (fracA + fracB) / 2; + geoM = interpolate(fracM); + m = transform(geoM); + if (ol.math.squaredSegmentDistance(m[0], m[1], a[0], a[1], + b[0], b[1]) < squaredTolerance) { + // If the m point is sufficiently close to the straight line, then we + // discard it. Just use the b coordinate and move on to the next line + // segment. + flatCoordinates.push(b[0], b[1]); + key = fracB.toString(); + goog.asserts.assert(!goog.object.containsKey(fractions, key)); + goog.object.set(fractions, key, true); + } else { + // Otherwise, we need to subdivide the current line segment. Split it + // into two and push the two line segments onto the stack. + fractionStack.push(fracB, fracM, fracM, fracA); + stack.push(b, m, m, a); + geoStack.push(geoB, geoM, geoM, geoA); + } + } + goog.asserts.assert(maxIterations > 0); + + return flatCoordinates; +}; + + +/** +* Generate a great-circle arcs between two lat/lon points. +* @param {number} lon1 Longitude 1 in degrees. +* @param {number} lat1 Latitude 1 in degrees. +* @param {number} lon2 Longitude 2 in degrees. +* @param {number} lat2 Latitude 2 in degrees. + * @param {ol.proj.Projection} projection Projection. +* @param {number} squaredTolerance Squared tolerance. +* @return {Array.} Flat coordinates. +*/ +ol.geom.flat.geodesic.greatCircleArc = function( + lon1, lat1, lon2, lat2, projection, squaredTolerance) { + + var geoProjection = ol.proj.get('EPSG:4326'); + + var cosLat1 = Math.cos(goog.math.toRadians(lat1)); + var sinLat1 = Math.sin(goog.math.toRadians(lat1)); + var cosLat2 = Math.cos(goog.math.toRadians(lat2)); + var sinLat2 = Math.sin(goog.math.toRadians(lat2)); + var cosDeltaLon = Math.cos(goog.math.toRadians(lon2 - lon1)); + var sinDeltaLon = Math.sin(goog.math.toRadians(lon2 - lon1)); + var d = sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosDeltaLon; + + return ol.geom.flat.geodesic.line_( + /** + * @param {number} frac Fraction. + * @return {ol.Coordinate} Coordinate. + */ + function(frac) { + if (1 <= d) { + return [lon2, lat2]; + } + var D = frac * Math.acos(d); + var cosD = Math.cos(D); + var sinD = Math.sin(D); + var y = sinDeltaLon * cosLat2; + var x = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosDeltaLon; + var theta = Math.atan2(y, x); + var lat = Math.asin(sinLat1 * cosD + cosLat1 * sinD * Math.cos(theta)); + var lon = goog.math.toRadians(lon1) + + Math.atan2(Math.sin(theta) * sinD * cosLat1, + cosD - sinLat1 * Math.sin(lat)); + return [goog.math.toDegrees(lon), goog.math.toDegrees(lat)]; + }, ol.proj.getTransform(geoProjection, projection), squaredTolerance); +}; + + +/** + * Generate a meridian (line at constant longitude). + * @param {number} lon Longitude. + * @param {number} lat1 Latitude 1. + * @param {number} lat2 Latitude 2. + * @param {ol.proj.Projection} projection Projection. + * @param {number} squaredTolerance Squared tolerance. + * @return {Array.} Flat coordinates. + */ +ol.geom.flat.geodesic.meridian = + function(lon, lat1, lat2, projection, squaredTolerance) { + var epsg4326Projection = ol.proj.get('EPSG:4326'); + return ol.geom.flat.geodesic.line_( + /** + * @param {number} frac Fraction. + * @return {ol.Coordinate} Coordinate. + */ + function(frac) { + return [lon, lat1 + ((lat2 - lat1) * frac)]; + }, + ol.proj.getTransform(epsg4326Projection, projection), squaredTolerance); +}; + + +/** + * Generate a parallel (line at constant latitude). + * @param {number} lat Latitude. + * @param {number} lon1 Longitude 1. + * @param {number} lon2 Longitude 2. + * @param {ol.proj.Projection} projection Projection. + * @param {number} squaredTolerance Squared tolerance. + * @return {Array.} Flat coordinates. + */ +ol.geom.flat.geodesic.parallel = + function(lat, lon1, lon2, projection, squaredTolerance) { + var epsg4326Projection = ol.proj.get('EPSG:4326'); + return ol.geom.flat.geodesic.line_( + /** + * @param {number} frac Fraction. + * @return {ol.Coordinate} Coordinate. + */ + function(frac) { + return [lon1 + ((lon2 - lon1) * frac), lat]; + }, + ol.proj.getTransform(epsg4326Projection, projection), squaredTolerance); +};