Minimal Bounding Circle Area of Lon/Lat Coordinates

How do you check if the GPS trajectory stays around an area for a certain amount of time? While working with GPS data in Sakay, this is the problem we’ve encountered when trying to detect if a vehicle is not moving. There are many ways to do this, but we decided to check if the bounding circle of the GPS traces within the past five minutes is less than a certain number of square meters.

Bounding Region Selection

So why a bounding circle, and not other bounding regions such as an axis-aligned bounding box, oriented bounding box, or a convex hull? This is to be able to filter out cases of vehicles traveling along a straight line, which, if we resorted to any of those other methods, will result in a bounding area too tight, and will detect such behaviour as stopped.

For our specific case bounding circles are better than other bounding regions, such as axis-aligned bounding boxes, oriented bounding boxes, or convex hulls, because they can better filter out cases of vehicles traveling along a straight line. If we used any of the other methods, the bounding area would be too tight and would detect vehicles traveling along a straight line as stopped.

Choosing a Map Projection

GPS data is sent in EPSG:4326, which are longitude and latitude coordinates describing the angles on a spheroid (very similar to polar coordinates). So given the GPS trace, we can’t just use a minimum bounding circle algorithm over the points because the point’s coordinates describe angles and not distances. If we try to compute the area of the bounding circle, we get square radians, a unit that does not make sense in our context (unless we are dealing with solid degrees for lighting calculations, or something to that effect). To get something usable, we need to project the longitude and latitude into a Cartesian system, then do our computations in that space.

An aside: PostGIS has implemented a Minimum Bounding Circle computation since 2.3. However, it only accepts geometry which means that the calculation is done in Cartesian space. So just passing geography data without choosing an appropriate projection may give you wrong results.

Unfortunately, there’s no perfect way to turn a sphere into a flat plane — there will be some stretching, rotation, or warping somewhere, making things look larger than they actually are, or creating some other distortions. There are multiple map projections that have been made with different trade-offs for various purposes. They can be categorized as follows:

• Aesthetic/Compromise: The map looks nice but no property of the shapes is formally preserved (though they may be, to a limited extent).
• Equidistant: Distances between two points are preserved. This also means that computing the distance between two points on the projected map using the Pythagorean distance formula will result in the distance between to points on the surface of the earth, or at least be off by a constant factor. Though distances may be preserved, the shapes may be sheared to preserve these properties. Note however, that this is only possible for one or two points per projection. That is, you select one or two points (depending on the projection) as anchors, where among all distances computed with those points, one of the endpoints will be accurate, but the rest will not.
• Equal-area: Areas enclosing a space is preserved. This is a “weaker” version of the equidistant projections but still works in most cases.
• Conformal: Angles are preserved within a local area. This usually means that shapes are not warped, but this comes at the expense of preserving distances or areas.
• Gnomonic: Straight lines between two points map to the shortest path between two points along the surface of the earth (i.e. “great circle” paths).

In our case, where we need to do calculations to filter points in a local area, what we use needs to be:

1. Fast — we need to process thousands of GPS data points per second, so this means that any method that requires the solving of an implicit equation, or carrying out integration, is out.
2. Projected units is in meters — When distance or area is computed on the projected points, the result should be in meters. This makes controlling the parameters of the filter easier, since it maps out to physical measures directly instead of some abstract measure. In some cases this is not explicitly mentioned, but one way to find out if the projection units is in meters, is if the projection formula contains a scaling by the radius of the earth.
3. Equal-area, optionally equidistant — Since we are computing the area, this is our priority even at the expense of angles not being maintained, since we don’t actually care about angles.

While the go-to choice for map projection is the Universal Transverse Mercator (UTM), it’s unnecessarily complicated for our case. So based on the criteria listed, we limit our choices to the following projections:

1. Gall-Peters
2. Collignon
3. Equirectangular

We ended up going with the Collignon projection since the formula is straightforward, and being able to choose a central meridian enables us to localize the projection. Once projected, we can now directly use it to compute the area of the minimal bounding circle. For this case, we went for the approximate solution of deriving the area by:

1. Computing the axis-aligned bounding box (AABB) of the trajectory,
2. Retrieving the length of the longer side of box, then,
3. Using half the length as the radius to compute the area of the circle.

This solution is inaccurate (but fast) and will overestimate the area, but it works for our use case. For a more accurate computation, Wezl’s algorithm may be used.

Implementation of the Algorithm and Demonstration

Here’s a straightforward C++ snippet (I am using OSRM and Boost for geographic types and constants) for the Collignon projection

const auto SQRT_PI = sqrt(boost::math::constants::pi<double>());
typedef boost::geometry::model::d2::point_xy<double> projected_pt_t;
projected_pt_t collignon_projection(const osrm::util::FloatCoordinate& c, const osrm::util::FloatLongitude base_longitude) {
	using namespace osrm::util::coordinate_calculation::detail;
	const auto s = sqrt(1 - sin(static_cast<double>(c.lat * DEGREE_TO_RAD)));
	return projected_pt_t(
		(2/SQRT_PI) * EARTH_RADIUS * DEGREE_TO_RAD * static_cast<double>(c.lon - base_longitude) * s,
		SQRT_PI * EARTH_RADIUS * (1 - s)
	);
}

The Collignon Projection needs a central meridian (base_longitude). While it’s tempting to set it to the prime meridian (that is, longitude = 0), we will get more accurate results by setting it in the vicinity of where we are doing our calculations. This logic applies to all projection methods that need you to supply a longitude and/or latitude for which to center their projections on. This is because calculations are accurate near the center, and gets worse as you move away from it. In our case, the points form a trajectory so we can assume that they are close to each other (unless we have a teleporting vehicle) and thus, we can either choose the longitude from any point in the set, or compute the average longitude (there’s more analysis on what’s more appropriate, but we can choose any method).

For a simpler demonstration of how it works, suppose we have the following coordinates, forming a triangle: (121.0204195 14.56064563), (121.0207272 14.56117396), (121.0201416 14.56122199). Computing the area using QGIS with Ellipsoidal formula is 1932.51m². Using our method, we compute the Collignon projection of all three points with the base meridian set to the longitude of the first point, 121.0204195, then use the resulting three points to compute the area of the triangle using the determinant formula. The error is less than half of a percent. Not bad!

Note however, for other applications, such as real estate, the 8 square meter difference may still matter a lot. So a more complex projection method such as UTM in cases where areas are intended to be subdivided.

Conclusion

We present a quick way to computing areas given longitude/latitude coordinates. The accuracy is not the best but it works for basic use cases like filtering data points.

Need development or consulting for Geo-spatial work? We at Bonito Tech are offering various services to fit your needs. Our team has over a decade of experience creating solutions and transforming businesses through technology.

Let’s craft something awesome together.