diff --git a/lib/OpenLayers/Geometry/LineString.js b/lib/OpenLayers/Geometry/LineString.js index 72ecbb2549..28fc9ab2ed 100644 --- a/lib/OpenLayers/Geometry/LineString.js +++ b/lib/OpenLayers/Geometry/LineString.js @@ -640,3 +640,139 @@ OpenLayers.Geometry.LineString = OpenLayers.Class(OpenLayers.Geometry.Curve, { CLASS_NAME: "OpenLayers.Geometry.LineString" }); + + +/** + * Function: OpenLayers.Geometry.LineString.geodesic + * + * Parameters: + * interpolate - {function(number): OpenLayers.Geometry.Point} Interpolate + * function. + * transform - {function(OpenLayers.Geometry.Point): OpenLayers.Geometry.Point} + * Transform from longitude/latitude to projected coordinates. + * squaredTolerance - {number} Squared tolerance. + * + * Returns: + * {OpenLayers.Geometry.LineString} + */ +OpenLayers.Geometry.LineString.geodesic = + function(interpolate, transform, squaredTolerance) { + // FIXME reduce garbage generation + // FIXME optimize stack operations + + var components = []; + + var geoA = interpolate(0); + var geoB = interpolate(1); + + var a = transform(geoA); + var b = transform(geoB); + + var geoStack = [geoB, geoA]; + var stack = [b, a]; + var fractionStack = [1, 0]; + + 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 (!(key in fractions)) { + components.push(a); + 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 (OpenLayers.Geometry.distanceSquaredToSegment(m, {x1: a.x, y1: a.y, + x2: b.x, y2: b.y}).distance < 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. + components.push(b); + key = fracB.toString(); + 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); + } + } + + return new OpenLayers.Geometry.LineString(components); +}; + + +/** + * Function: OpenLayers.Geometry.LineString.geodesicMeridian + * Generate a meridian (line at constant longitude). + * + * Parameters: + * lon - {number} Longitude. + * lat1 - {number} Latitude 1. + * lat2 - {number} Latitude 2. + * projection - {OpenLayers.Projection} Projection. + * squaredTolerance - {number} Squared tolerance. + * + * Returns: + * {OpenLayers.Geometry.LineString} Line geometry for the meridian at . + */ +OpenLayers.Geometry.LineString.geodesicMeridian = + function(lon, lat1, lat2, projection, squaredTolerance) { + var epsg4326Projection = new OpenLayers.Projection('EPSG:4326'); + return OpenLayers.Geometry.LineString.geodesic( + function(frac) { + return new OpenLayers.Geometry.Point( + lon, lat1 + ((lat2 - lat1) * frac)); + }, + function(point) { + return point.transform(epsg4326Projection, projection); + }, + squaredTolerance + ); +}; + + +/** + * Function: OpenLayers.Geometry.LineString.geodesicParallel + * Generate a parallel (line at constant latitude). + * + * Parameters: + * lat - {number} Latitude. + * lon1 - {number} Longitude 1. + * lon2 - {number} Longitude 2. + * projection {OpenLayers.Projection} Projection. + * squaredTolerance - {number} Squared tolerance. + * + * Returns: + * {OpenLayers.Geometry.LineString} Line geometry for the parallel at . + */ +OpenLayers.Geometry.LineString.geodesicParallel = + function(lat, lon1, lon2, projection, squaredTolerance) { + var epsg4326Projection = new OpenLayers.Projection('EPSG:4326'); + return OpenLayers.Geometry.LineString.geodesic( + function(frac) { + return new OpenLayers.Geometry.Point( + lon1 + ((lon2 - lon1) * frac), lat); + }, + function(point) { + return point.transform(epsg4326Projection, projection); + }, + squaredTolerance + ); +}; +