Permalink
Browse files

adding douglas puecker simplification algorithm

  • Loading branch information...
1 parent 80b17ef commit 5bb90869fa23c67f5ae41c3d165241c75e3fcb5e @maxogden maxogden committed Apr 19, 2011
Showing with 164 additions and 55 deletions.
  1. +164 −55 geojson-utils.js
View
@@ -8,98 +8,207 @@
// adapted from http://www.kevlindev.com/gui/math/intersection/Intersection.js
gju.lineStringsIntersect = function(l1, l2) {
- var intersects = [];
- for (var i = 0; i <= l1.coordinates.length - 2; ++i) {
- for (var j = 0; j <= l2.coordinates.length - 2; ++j) {
- var a1 = {x: l1.coordinates[i][1], y: l1.coordinates[i][0]},
- a2 = {x: l1.coordinates[i+1][1], y: l1.coordinates[i+1][0]},
- b1 = {x: l2.coordinates[j][1], y: l2.coordinates[j][0]},
- b2 = {x: l2.coordinates[j+1][1], y: l2.coordinates[j+1][0]},
- ua_t = (b2.x - b1.x) * (a1.y - b1.y) - (b2.y - b1.y) * (a1.x - b1.x),
- ub_t = (a2.x - a1.x) * (a1.y - b1.y) - (a2.y - a1.y) * (a1.x - b1.x),
- u_b = (b2.y - b1.y) * (a2.x - a1.x) - (b2.x - b1.x) * (a2.y - a1.y);
- if ( u_b != 0 ) {
- var ua = ua_t / u_b,
- ub = ub_t / u_b;
- if ( 0 <= ua && ua <= 1 && 0 <= ub && ub <= 1 ) {
- intersects.push({
- 'type': 'Point',
- 'coordinates': [a1.x + ua * (a2.x - a1.x), a1.y + ua * (a2.y - a1.y)]
- });
- }
- }
+ var intersects = [];
+ for (var i = 0; i <= l1.coordinates.length - 2; ++i) {
+ for (var j = 0; j <= l2.coordinates.length - 2; ++j) {
+ var a1 = {x: l1.coordinates[i][1], y: l1.coordinates[i][0]},
+ a2 = {x: l1.coordinates[i+1][1], y: l1.coordinates[i+1][0]},
+ b1 = {x: l2.coordinates[j][1], y: l2.coordinates[j][0]},
+ b2 = {x: l2.coordinates[j+1][1], y: l2.coordinates[j+1][0]},
+ ua_t = (b2.x - b1.x) * (a1.y - b1.y) - (b2.y - b1.y) * (a1.x - b1.x),
+ ub_t = (a2.x - a1.x) * (a1.y - b1.y) - (a2.y - a1.y) * (a1.x - b1.x),
+ u_b = (b2.y - b1.y) * (a2.x - a1.x) - (b2.x - b1.x) * (a2.y - a1.y);
+ if ( u_b != 0 ) {
+ var ua = ua_t / u_b,
+ ub = ub_t / u_b;
+ if ( 0 <= ua && ua <= 1 && 0 <= ub && ub <= 1 ) {
+ intersects.push({
+ 'type': 'Point',
+ 'coordinates': [a1.x + ua * (a2.x - a1.x), a1.y + ua * (a2.y - a1.y)]
+ });
}
}
- if (intersects.length == 0) intersects = false;
+ }
+ }
+ if (intersects.length == 0) intersects = false;
return intersects;
}
// adapted from http://jsfromhell.com/math/is-point-in-poly
gju.pointInPolygon = function(point, polygon) {
- var x = point.coordinates[1],
- y = point.coordinates[0],
- poly = polygon.coordinates[0]; //TODO: support polygons with holes
- for (var c = false, i = -1, l = poly.length, j = l - 1; ++i < l; j = i) {
- var px = poly[i][1], py = poly[i][0],
- jx = poly[j][1], jy = poly[j][0];
- if (((py <= y && y < jy) || (jy <= y && y < py)) && (x < (jx - px) * (y - py) / (jy - py) + px)) {
- c = [point];
- }
+ var x = point.coordinates[1],
+ y = point.coordinates[0],
+ poly = polygon.coordinates[0]; //TODO: support polygons with holes
+ for (var c = false, i = -1, l = poly.length, j = l - 1; ++i < l; j = i) {
+ var px = poly[i][1], py = poly[i][0],
+ jx = poly[j][1], jy = poly[j][0];
+ if (((py <= y && y < jy) || (jy <= y && y < py)) && (x < (jx - px) * (y - py) / (jy - py) + px)) {
+ c = [point];
}
- return c;
+ }
+ return c;
}
gju.numberToRadius = function(number) {
- return number * Math.PI / 180;
+ return number * Math.PI / 180;
}
gju.numberToDegree = function(number) {
- return number * 180 / Math.PI;
+ return number * 180 / Math.PI;
}
// written with help from @tautologe
gju.drawCircle = function(radiusInMeters, centerPoint) {
var center = [centerPoint.coordinates[1], centerPoint.coordinates[0]],
- dist = (radiusInMeters / 1000) / 6371, // convert meters to radiant
- radCenter = [gju.numberToRadius(center[0]), gju.numberToRadius(center[1])],
- steps = 15, // 15 sided circle
- poly = [[center[0], center[1]]];
+ dist = (radiusInMeters / 1000) / 6371, // convert meters to radiant
+ radCenter = [gju.numberToRadius(center[0]), gju.numberToRadius(center[1])],
+ steps = 15, // 15 sided circle
+ poly = [[center[0], center[1]]];
for (var i = 0; i < steps + 1; i++) {
var brng = 2 * Math.PI * i / steps;
var lat = Math.asin(Math.sin(radCenter[0]) * Math.cos(dist) +
- Math.cos(radCenter[0]) * Math.sin(dist) * Math.cos(brng));
+ Math.cos(radCenter[0]) * Math.sin(dist) * Math.cos(brng));
var lng = radCenter[1] + Math.atan2(Math.sin(brng) * Math.sin(dist) *
- Math.cos(radCenter[0]),
- Math.cos(dist) - Math.sin(radCenter[0]) *
- Math.sin(lat));
- poly[i] = [];
- poly[i][1] = gju.numberToDegree(lat);
- poly[i][0] = gju.numberToDegree(lng);
+ Math.cos(radCenter[0]),
+ Math.cos(dist) - Math.sin(radCenter[0]) *
+ Math.sin(lat));
+ poly[i] = [];
+ poly[i][1] = gju.numberToDegree(lat);
+ poly[i][0] = gju.numberToDegree(lng);
}
- return { "type": "Polygon",
- "coordinates": [poly] };
+ return { "type": "Polygon", "coordinates": [poly] };
}
gju.rectangleCentroid = function(rectangle) {
var bbox = rectangle.coordinates[0];
var xmin = bbox[0][0], ymin = bbox[0][1], xmax = bbox[1][0], ymax = bbox[1][1];
var xwidth = xmax - xmin;
var ywidth = ymax - ymin;
- return { 'type': 'Point',
- 'coordinates': [xmin + xwidth/2, ymin + ywidth/2] };
+ return { 'type': 'Point', 'coordinates': [xmin + xwidth/2, ymin + ywidth/2] };
}
// from http://www.movable-type.co.uk/scripts/latlong.html
gju.pointDistance = function(pt1, pt2) {
var lon1 = pt1.coordinates[0], lat1 = pt1.coordinates[1],
- lon2 = pt2.coordinates[0], lat2 = pt2.coordinates[1],
- dLat = gju.numberToRadius(lat2 - lat1),
- dLon = gju.numberToRadius(lon2 - lon1),
- a = Math.sin(dLat/2) * Math.sin(dLat/2) +
- Math.cos(gju.numberToRadius(lat1)) * Math.cos(gju.numberToRadius(lat2)) *
- Math.sin(dLon/2) * Math.sin(dLon/2),
- c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
+ lon2 = pt2.coordinates[0], lat2 = pt2.coordinates[1],
+ dLat = gju.numberToRadius(lat2 - lat1),
+ dLon = gju.numberToRadius(lon2 - lon1),
+ a = Math.sin(dLat/2) * Math.sin(dLat/2) +
+ Math.cos(gju.numberToRadius(lat1)) * Math.cos(gju.numberToRadius(lat2)) *
+ Math.sin(dLon/2) * Math.sin(dLon/2),
+ c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return (6371 * c) * 1000; // returns meters
}
+ gju.simplify = function (source, kink) {
+ /* source[] array of geojson points */
+ /* kink in metres, kinks above this depth kept */
+ /* kink depth is the height of the triangle abc where a-b and b-c are two consecutive line segments */
+ kink = kink || 20;
+ source = source.map(function(o) { return {lng: o.coordinates[0], lat: o.coordinates[1]} });
+
+ var n_source, n_stack, n_dest, start, end, i, sig;
+ var dev_sqr, max_dev_sqr, band_sqr;
+ var x12, y12, d12, x13, y13, d13, x23, y23, d23;
+ var F = ((Math.PI / 180.0) * 0.5 );
+ var index = new Array(); /* aray of indexes of source points to include in the reduced line */
+ var sig_start = new Array(); /* indices of start & end of working section */
+ var sig_end = new Array();
+
+ /* check for simple cases */
+
+ if ( source.length < 3 )return(source); /* one or two points */
+
+ /* more complex case. initialize stack */
+
+ n_source = source.length;
+ band_sqr = kink * 360.0 / (2.0 * Math.PI * 6378137.0); /* Now in degrees */
+ band_sqr *= band_sqr;
+ n_dest = 0;
+ sig_start[0] = 0;
+ sig_end[0] = n_source-1;
+ n_stack = 1;
+
+ /* while the stack is not empty ... */
+ while ( n_stack > 0 ){
+
+ /* ... pop the top-most entries off the stacks */
+
+ start = sig_start[n_stack-1];
+ end = sig_end[n_stack-1];
+ n_stack--;
+
+ if ( (end - start) > 1 ){ /* any intermediate points ? */
+
+ /* ... yes, so find most deviant intermediate point to
+ either side of line joining start & end points */
+
+ x12 = (source[end].lng() - source[start].lng());
+ y12 = (source[end].lat() - source[start].lat());
+ if (Math.abs(x12) > 180.0)
+ x12 = 360.0 - Math.abs(x12);
+ x12 *= Math.cos(F * (source[end].lat() + source[start].lat()));/* use avg lat to reduce lng */
+ d12 = (x12*x12) + (y12*y12);
+
+ for ( i = start + 1, sig = start, max_dev_sqr = -1.0; i < end; i++ ){
+
+ x13 = (source[i].lng() - source[start].lng());
+ y13 = (source[i].lat() - source[start].lat());
+ if (Math.abs(x13) > 180.0)
+ x13 = 360.0 - Math.abs(x13);
+ x13 *= Math.cos (F * (source[i].lat() + source[start].lat()));
+ d13 = (x13*x13) + (y13*y13);
+
+ x23 = (source[i].lng() - source[end].lng());
+ y23 = (source[i].lat() - source[end].lat());
+ if (Math.abs(x23) > 180.0)
+ x23 = 360.0 - Math.abs(x23);
+ x23 *= Math.cos(F * (source[i].lat() + source[end].lat()));
+ d23 = (x23*x23) + (y23*y23);
+
+ if ( d13 >= ( d12 + d23 ) )
+ dev_sqr = d23;
+ else if ( d23 >= ( d12 + d13 ) )
+ dev_sqr = d13;
+ else
+ dev_sqr = (x13 * y12 - y13 * x12) * (x13 * y12 - y13 * x12) / d12;// solve triangle
+
+ if ( dev_sqr > max_dev_sqr ){
+ sig = i;
+ max_dev_sqr = dev_sqr;
+ }
+ }
+
+ if ( max_dev_sqr < band_sqr ){ /* is there a sig. intermediate point ? */
+ /* ... no, so transfer current start point */
+ index[n_dest] = start;
+ n_dest++;
+ }
+ else{
+ /* ... yes, so push two sub-sections on stack for further processing */
+ n_stack++;
+ sig_start[n_stack-1] = sig;
+ sig_end[n_stack-1] = end;
+ n_stack++;
+ sig_start[n_stack-1] = start;
+ sig_end[n_stack-1] = sig;
+ }
+ }
+ else{
+ /* ... no intermediate points, so transfer current start point */
+ index[n_dest] = start;
+ n_dest++;
+ }
+ }
+
+ /* transfer last point */
+ index[n_dest] = n_source-1;
+ n_dest++;
+
+ /* make return array */
+ var r = new Array();
+ for(var i=0; i < n_dest; i++)
+ r.push(source[index[i]]);
+ return r.map(function(o) { return {type: "Point", coordinates: [o.lng, o.lat]} });;
+ }
})();

0 comments on commit 5bb9086

Please sign in to comment.