Sometimes my job throws an interesting problem my way. This week I was presented with a very odd geometry problem 🙂

I needed to generate KML files from geographic data and one of my requirements was to represent certain geographic areas as polygons, the vertices of which would be supplied (along with the rest of the data) by another OSGi service running elsewhere. Seems fairly straightforward — In KML, the vertices of a polygon are usually specified as follows:

<Placemark> <name>LinearRing.kml</name> <Polygon> <outerBoundaryIs> <LinearRing> <coordinates> -122.365662,37.826988,0 -122.365202,37.826302,0 -122.364581,37.82655,0 -122.365038,37.827237,0 -122.365662,37.826988,0 </coordinates> </LinearRing> </outerBoundaryIs> </Polygon> </Placemark>

The coordinates tag requires at least four longitude/latitude/altitude triples to be specified, with the last coordinate being the same as the first. Here is where the problem comes in — The order in which these coordinates are specified matters (they must be specified in counter-clockwise order). To mix up the order of the coordinates would have unpredictable results (e.g. crazy geometric shapes) when the data is later displayed via Google Earth (or some other application that supports KML files). However, the area vertices are indeed fed to my KML generator in no particular order (and the services providing the data cannot be changed to guarantee a particular ordering).

So… how do I put the points in order? “Surely this is a solved problem.” I thought, turning to the all-knowing internet. A bit of searching turned up an algorithm called Graham’s Scan. Basically, this algorithm takes a bag of random coordinates and generates a convex hull with vertices defined in counter-clockwise order (Note: This may not be suitable if you’re trying to faithfully recreate complex geometries, fortunately I’m mostly concerned with rectangular areas). Roughly, the algorithm works as follows:

- Find the coordinate with the lowest y-value.
- Sort the remaining points by the polar angle between the line defined by the current point and the point found in step 1, and the x-axis.
- Go through the list, evaluating 3 points at a time (you’re concerned with the angle of the turn made by the two resulting line segments). Eliminate any non counter-clockwise turns (i.e., non-convex corners).

I found several example implementations for this algorithm in various languages: C++, Java, etc. Since I’m coding this up in Scala, I wasn’t too happy with any of those, and I couldn’t find an example in Scala to ~~rip off~~ draw inspiration from. However, I did manage to find an implementation in Haskell which I used as a rough guide. Anyway, here’s my attempt at Graham’s Scan in Scala:

type Coordinate = (Double, Double) //Longitude, Latitude protected def processArea(coords: List[Coordinate]): List[Coordinate] = { //returns > 0 if points form a counter clockwise turn, // < 0 if clockwise, and 0 if collinear def ccw(p1: Coordinate, p2: Coordinate, p3: Coordinate) = (p2._1 - p1._1)*(p3._2 - p1._2) - (p2._2 - p1._2)*(p3._1 - p1._1) //Scan the List of coordinates and find our vertices def scan(theCoords: List[Coordinate]): List[Coordinate] = theCoords match { case xs if xs.isEmpty => List() case xs if xs.size == 2 => theCoords case x::y::z::xs if ccw(x,y,z) > 0 => x::scan(y::z::xs) case x::y::z::xs => scan(x::z::xs) } //find the coordinate with the lowest latitude val origin = coords.minBy(_._2) //sort the rest of the points according to their polar angle (the angle between the line //defined by the origin and the current point, and the x-axis) val coordList = origin :: coords.filterNot(_ == origin). sortBy(point => atan2(point._2 - origin._2, point._1 - origin._1)) //do the graham scan scan(coordList) }

I think you’ll find that’s a bit shorter than some of the imperative language implementations out there 🙂

Hey man. I had to solve the same problem, your solution helped me to get some inspiration, but I must tell you, is wrong. After a couple of right turns, and then a left again, it “forgets” the last correct point, and degenerates. You should plot it in Octave/Matlab to check it, anyway here is my solution:

https://github.com/caente/convex-hull/blob/master/src/main/scala/com/miguel/GrahamScanScala.scala

BTW there is also a solution in Java 8 🙂

LikeLike

Interesting. Don’t believe I noted that in my testing, but then I’m not having to deal with a very large number of points for my particular problem. But I’m definitely curious now, so I’m looking at it.

In my final implementation, I did have to make some changes to this for geometries that straddled the anti-meridian (-180/180 longitude). This algorithm definitely did not deal well with that situation, but that’s a special case that people using a normal, cartesian coordinate system wouldn’t run into.

LikeLike

You mean points with the same X? That should be taken care of with the isLeft function, you only take the middle point if there is a left angle, ditch it otherwise…

LikeLike

No, the issue came up when displaying areas in something like google earth. An area that straddled the anti-meridian (roughly the international date line), wouldn’t display correctly, because the coordinate system abruptly went from -180.0 to 180.0 at that line. The solution was to split such polygons into two smaller ones, which was quite involved code-wise and involved calculating and adding new points.

LikeLike