Sort latitude and longitude coordinates into clockwise ordered quadrilateral

Question

Problem

Users can provide up to four latitude and longitude coordinates, in any order. They do so with Google Maps. Using Google's Polygon API (v3), the coordinates they select should highlight the selected area between the four coordinates.

Question

How do you sort an array of latitude and longitude coordinates in (counter-)clockwise order?

Solutions and Searches

StackOverflow Questions

• Drawing resizable (not intersecting) polygons
• How to sort points in a Google maps polygon so that lines do not cross?
• Sort Four Points in Clockwise Order

Related Sites

• http://www.daftlogic.com/projects-google-maps-area-calculator-tool.htm
• http://en.literateprograms.org/Quickhull_%28Javascript%29
• http://www.geocodezip.com/map-markers_ConvexHull_Polygon.asp
• http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm

Known Algorithms

• Graham's scan (too complicated)
• Jarvis March algorithm (handles N points)
• Recursive Convex Hull (removes a point)

Code

Here is what I have so far:

// Ensures the markers are sorted: NW, NE, SE, SW
function sortMarkers() {
  var ns = markers.slice( 0 );
  var ew = markers.slice( 0 );

  ew.sort( function( a, b ) {
    if( a.position.lat() < b.position.lat() ) {
      return -1;
    }
    else if( a.position.lat() > b.position.lat() ) {
      return 1;
    }

    return 0;
  });

  ns.sort( function( a, b ) {
    if( a.position.lng() < b.position.lng() ) {
      return -1;
    }
    else if( a.position.lng() > b.position.lng() ) {
      return 1;
    }

    return 0;
  });

  var nw;
  var ne;
  var se;
  var sw;

  if( ew.indexOf( ns[0] ) > 1 ) {
    nw = ns[0];
  }
  else {
    ne = ns[0];
  }

  if( ew.indexOf( ns[1] ) > 1 ) {
    nw = ns[1];
  }
  else {
    ne = ns[1];
  }

  if( ew.indexOf( ns[2] ) > 1 ) {
    sw = ns[2];
  }
  else {
    se = ns[2];
  }

  if( ew.indexOf( ns[3] ) > 1 ) {
    sw = ns[3];
  }
  else {
    se = ns[3];
  }

  markers[0] = nw;
  markers[1] = ne;
  markers[2] = se;
  markers[3] = sw;
}

Thank you.

Solution 1

Given the points:

   4  +        [d]            [g]                 
      |                             
   3 [a]            [e]                 
      |                             
   2  +                  [f]       [h]    
      |                             
   1  +   [b]                             
      |                             
   0  +----+---[c]---+----+----+----+
      0    1    2    3    4    5    6

you want to find the following bound walk:

   4  +     ___[d]------------[g]                 
      |  __/                     \    
   3 [a]/           [e]__         \       
      | \             \_ ```---    \  
   2  +  \              `[f]   \___[h]    
      |   \           __/            
   1  +   [b]      __/                   
      |      \    /                
   0  +----+--`[c]---+----+----+----+
      0    1    2    3    4    5    6

?

If this is correct, here's a way:

• find the upper most point, P_top, in the set of points. In case of a tie, pick the point with the smallest x coordinate
• sort all points by comparing the slopes m_i and m_j of the lines each pair of points (excluding P_top!) P_i and P_j make when passing through P_top
  • if m_i and m_j are equal, let the point P_i or P_j closest to P_top come first
  • if m_i is positive and m_j is negative (or zero), P_j comes first
  • if both m_i and m_j are either positive or negative, let the point belonging to the line with the largest slope come first

Here's a quick demo for the map:

(I know little JavaScript, so I might, or probably have, violated some JavaScript code conventions...):

var points = [
    new Point("Stuttgard", 48.7771056, 9.1807688),
    new Point("Rotterdam", 51.9226899, 4.4707867),
    new Point("Paris", 48.8566667, 2.3509871),
    new Point("Hamburg", 53.5538148, 9.9915752),
    new Point("Praha", 50.0878114, 14.4204598),
    new Point("Amsterdam", 52.3738007, 4.8909347),
    new Point("Bremen", 53.074981, 8.807081),
    new Point("Calais", 50.9580293, 1.8524129),
];
var upper = upperLeft(points);

print("points :: " + points);
print("upper  :: " + upper);
points.sort(pointSort);
print("sorted :: " + points);

// A representation of a 2D Point.
function Point(label, lat, lon) {

    this.label = label;
    this.x = (lon + 180) * 360;
    this.y = (lat + 90) * 180;

    this.distance=function(that) {
        var dX = that.x - this.x;
        var dY = that.y - this.y;
        return Math.sqrt((dX*dX) + (dY*dY));
    }

    this.slope=function(that) {
        var dX = that.x - this.x;
        var dY = that.y - this.y;
        return dY / dX;
    }

    this.toString=function() {
        return this.label;
    }
}

// A custom sort function that sorts p1 and p2 based on their slope
// that is formed from the upper most point from the array of points.
function pointSort(p1, p2) {
    // Exclude the 'upper' point from the sort (which should come first).
    if(p1 == upper) return -1;
    if(p2 == upper) return 1;

    // Find the slopes of 'p1' and 'p2' when a line is 
    // drawn from those points through the 'upper' point.
    var m1 = upper.slope(p1);
    var m2 = upper.slope(p2);

    // 'p1' and 'p2' are on the same line towards 'upper'.
    if(m1 == m2) {
        // The point closest to 'upper' will come first.
        return p1.distance(upper) < p2.distance(upper) ? -1 : 1;
    }

    // If 'p1' is to the right of 'upper' and 'p2' is the the left.
    if(m1 <= 0 && m2 > 0) return -1;

    // If 'p1' is to the left of 'upper' and 'p2' is the the right.
    if(m1 > 0 && m2 <= 0) return 1;

    // It seems that both slopes are either positive, or negative.
    return m1 > m2 ? -1 : 1;
}

// Find the upper most point. In case of a tie, get the left most point.
function upperLeft(points) {
    var top = points[0];
    for(var i = 1; i < points.length; i++) {
        var temp = points[i];
        if(temp.y > top.y || (temp.y == top.y && temp.x < top.x)) {
            top = temp;
        }
    }
    return top;
}

Note: your should double, or triple check the conversions from lat,lon to x,y as I am a novice if it comes to GIS!!! But perhaps you don't even need to convert anything. If you don't, the upperLeft function might just return the lowest point instead of the highest, depending on the locations of the points in question. Again: triple check these assumptions!

When executing the snippet above, the following gets printed:

points :: Stuttgard,Rotterdam,Paris,Hamburg,Praha,Amsterdam,Bremen,Calais
upper  :: Hamburg
sorted :: Hamburg,Praha,Stuttgard,Paris,Bremen,Calais,Rotterdam,Amsterdam

Alternate Distance Function

function distance(lat1, lng1, lat2, lng2) {
  var R = 6371; // km
  var dLat = (lat2-lat1).toRad();
  var dLon = (lng2-lng1).toRad();
  var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
          Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) *
          Math.sin(dLon/2) * Math.sin(dLon/2);
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
  return R * c;
}

Solution 2

Algorithm idea: average the four points to get a point inside the polygon. Then calculate the angle of the ray between that center point and each point, using inverse trigonometric functions, like explained here. Then sort by the angles. That should give you a (counter-)clockwise ordering, depending on the sort order and what you consider "zero degrees".

UPDATE: here's some code. Mostly untested, but it's the idea.

function sorted_points(points) {
    points = points.slice(0); // copy the array, since sort() modifies it
    var stringify_point = function(p) { return p.x + ',' + p.y; };

    // finds a point in the interior of `pts`
    var avg_points = function(pts) {
        var x = 0;
        y = 0;
        for(i = 0; i < pts.length; i++) {
            x += pts[i].x;
            y += pts[i].y;
        }
        return {x: x/pts.length, y:y/pts.length};
    }
    var center = avg_points(points);

    // calculate the angle between each point and the centerpoint, and sort by those angles
    var angles = {};
    for(i = 0; i < points.length; i++) {
        angles[stringify_point(points[i])] = Math.atan(points[i].x - center.x, points[i].y - center.y);
    }
    points.sort(function(p1, p2) {
        return angles[stringify_point(p1)] - angles[stringify_point(p2)];
    });
    return points;
}

It sorts points (an array of objects like {x: 1, y: 1}) counter-clockwise.

Solution 3

For those arriving here having a similar problem a year later:

I do not agree with the chosen answer's bound walk. There is no singular solution to the order even with a given clock direction. The convex hull of the given coordinates eliminates points e and f. These can then be attached anywhere along the path. Objectively, h,e,f,c can be improved to h,f,e,c keeping the direction of the x component consistent - in this case, negative.

The significance of this makes it impossible to guarantee the inclusion of any map location in the resulted area bounded by the chosen walk.

Comments

[Public Profile]

This is a brilliant solution, but maybe I'm confused by the original question. Assuming clockwise, it seems like Bremen is out of place. I would have expected a sort order of: Hamburg,Praha,Stuttgard,Paris,Calais,Rotterdam,Amsterdam,Bre‌​men

[Bart Kiers]

@Jon Stevens, no, it sorts as it's supposed to. Think of it like this: attach a piece of rope to the top-city (Hamburg) with a small weight at the end. Drag the weight all the way to the right so that the rope is straight and then let the weight go. The order in which the rope goes through the cities is the proper sort-order. HTH

[Public Profile]

Bart, I get that part, but from your ascii map, it seems the sort order would be more in line with that. Also, the original requirements were to be (counter)-clockwise. Not necessarily string/rope based.

[Bart Kiers]

@Jon Stevens, no, the ascii map is also sorted like that with d being the top, followed by g,h,e,f,c,b,a. The original requirements are not all that precise: the OP mentions "convex hulls" as well while that had nothing to do with it. All I know is that the OP said: "... and is exactly what I was looking to accomplish" and accepted this answer as the correct one.

[Lee Louviere]

The only way to guarantee inclusion is to create a polygon with the outer most points, select two points in that shape and indent the line to loop through the second-most outer points, and recursively repeat that process. The shape you end up with is similar to a labyrinth. But, the complexity of the OP is not that great. They suggested 4 points, and we know that 4 points is always a quadrilateral, no matter where the points are. Therefore the given solutions will work.

[Tran Quan]

Awesome 💰 so simple and nice

[cbn]

I don't know about your code, but the algorithm is the correct answer - certainly correct for four points in a QUAD, as requested by the original poster. This should be upvoted as the correct answer.