I have a program that takes as input an array of lat/long points. I need to perform a check on that array to ensure that all of the points are within a certain radius. So, for example, the maximum radius I will allow is 100 miles. Given an array of lat/long (coming from a MySQL database, could be 10 points could be 10000) I need to figure out if they will all fit in a circle with radius of 100 miles.
Kinda stumped on how to approach this. Any help would be greatly appreciated.
Find the smallest circle containing all points, and compare its radius to 100.
It's easiest way for me to solve this is by converting the coordinates to (X,Y,Z), then finding the distance along a sphere.
Assuming Earth is a sphere (totally untrue) with radius R...
X = R * cos(long) * cos(lat)
Y = R * sin(long) * cos(lat)
Z = R * sin(lat)
At this point, you can approximate the distance between the points using the extension of the pythagorean theorem for threespace:
dist = sqrt((x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2)
But to find the actual distance along the surface, you're going to need to know the angle subtended by the two points from the origin (center of the Earth).
Representing your locations as vectors V1 = (X1, Y1, Z1) and V2 = (X2, Y2, Z2), the angle is:
angle = arcsin((V1 x V2) / (|V1||V2|)), where x is the cross-product.
The distance is then:
dist = (Earth's circumference) * angle / (2 * pi)
Of course, this doesn't take into account changes in elevation or the fact that the Earth is wider at the equator.
Apologies for not writing my math in LaTeX.
The answer below involves pretending that the earth is a perfect sphere, which should give a more accurate answer than treating the earth as a flat plane.
To figure out the radius of a set of lat/lon points, you must first ensure that your set of points is "hemispherical", ie. all the points can fit into some arbitrary half of your perfect sphere.
Check out Section 3 in the paper "Optimal algorithms for some proximity problems on the Gaussian sphere with applications" by Gupta and Saluja. I don't have a specific link, but I believe that you can find a copy online for free. This paper is not sufficient to implement a solution. You'll also need Appendix 1 in "Approximating Centroids for the Maximum Intersection of Spherical Polygons" by Ha and Yoo.
I wouldn't use Megiddo's algorithm for doing the linear programming part of the hemisphericity testing. Instead, use Seidel's algorithm for solving Linear Programming problems, described in "Small-Dimensional Linear Programming and Convex Hulls Made Easy" by Raimund Seidel. Also see "Seidel’s Randomized Linear Programming Algorithm" by Kurt Mehlhorn and Section 9.4 from "Real-Time Collision Detection" by Christer Ericson.
Once you have determined that your points are hemispherical, move on to Section 4 of the paper by Gupta and Saluja. This part shows how to actually get the "smallest enclosing circle" for the points.
To do the required quadratic programming, see the paper "A Randomized Algorithm for Solving Quadratic Programs" by N.D. Botkin. This tutorial is helpful, but the paper uses (1/2)x^T G x - g^T x and the web tutorial uses (1/2)x^T H x + c^T x. One adds the terms and the other subtracts, leading to sign-related problems. Also see this example 2D QP problem. A hint: if you're using C++, the Eigen library is very good.
This method is a little more complicated than some of the 2D methods above, but it should give you more accurate results than just ignoring the curvature of the earth completely. This method also has O(n) time complexity, which is likely asymptotically optimal.
Note: The method described above may not handle duplicate data well, so you may want to check for duplicate lat/lon points before finding the smallest enclosing circle.
Check out the answers to this question. It gives a way to measure the distance between any two (lat,long) points. Then use a smallest enclosing circle algorithm.
I suspect that finding a smallest enclosing circle may be difficult enough on a plane, so to eliminate the subtleties of working with latitude and longitude and spherical geometry, you should probably consider mapping your points to the XY plane. That will introduce some amount of distortion, but if your intended scale is 100 miles you can probably live with that. Once you have a circle and its center on the XY plane, you can always map back to the terrestial sphere and re-check your distances.
