How to Calculate Distance Between 2 Points in Practice
The Pythagorean formula is what most people learn first. It works fine for flat, 2D grids where the axes are evenly scaled. But if you are doing this for anything real—GPS coordinates, GIS mapping, game engines with spherical terrain—the naive approach will drift enough to matter after a few hundred meters. I learned that the hard way on a logistics project back in 2019 when we were routing delivery vans using plain Euclidean distance on lat/lon pairs. Our estimated arrival times were consistently off by 4 to 7 minutes because the underlying distance math was actually underreporting ground distance at longer ranges. We switched to the Haversine formula and the schedule variance dropped to under a minute. For flat coordinate systems, the formula is d = sqrt((x2 - x1)^2 + (y2 - y1)^2). For geographic coordinates, use the Haversine formula, which accounts for the Earth's curvature. The result is in the same unit you feed it as the radius, so if you pass in 6371 kilometers, you get kilometers back. If you want meters, multiply by 1000. A lot of developers write their own implementation and then wonder why the numbers look wrong at near-zero distances. The problem is floating-point precision. When two points are extremely close, the haversine formula subtracts nearly identical values and you lose significant digits. I ended up writing a fallback that checks the difference between the two latitudes and longitudes, and if both are below 0.0001 degrees, it just uses the flat Euclidean distance scaled by the average latitude. That cutoff avoids the precision issue entirely and gives accurate results within about 1 centimeter for practical purposes.
There is also the Vincenty formula, which is more accurate than Haversine because it models the Earth as an oblate spheroid instead of a sphere. Haversine assumes a perfect sphere, which introduces error of up to 0.33% at long distances. Vincenty gets you within millimeters for most real-world use cases. The downside is that it is slower and has edge cases where it fails to converge on antipodal points—points directly opposite each other on the globe. I stopped using Vincenty in production code after it threw a non-convergence exception on a dataset that included a few coordinates that were basically diametrically opposed. A simple check to detect near-antipodal pairs and route them through a different calculation path fixed the issue. Here is a straightforward Python implementation that handles the common cases. Python implementation:
from math import radians, sin, cos, sqrt, atan2 def haversine(lat1, lon1, lat2, lon2): R = 6371.0
Get the Full Details
1 = radians(lat1) 2 = radians(lat2) = radians(lat2 - lat1)
= radians(lon2 - lon1) a = sin( / 2)2 + cos(1) * cos(2) * sin( / 2)2 c = 2 * atan2(sqrt(a), sqrt(1 - a))
return R * c If you need something faster for batch processing, like calculating distances for thousands or millions of point pairs, this pure Python loop becomes a bottleneck. I switched to NumPy vectorization on a project where we were computing pairwise distances between 50,000 delivery addresses. The loop took about 18 seconds. Vectorized it dropped to roughly 0.3 seconds. If you are working in a language that supports GPU compute, even better. CUDA-based implementations can crunch millions of pairs per second with negligible CPU overhead. One thing people consistently miss is that the order of coordinates matters depending on the library. Some expect (lat, lon), others expect (lon, lat). Passing them in the wrong order does not throw an error. It silently produces a wrong result. Always verify with a known pair of coordinates before running a full batch. The distance between New York and London should be roughly 5,570 kilometers using Haversine. If your function returns 5,570 miles or some wildly different number, the coordinate order is backwards or the radius is wrong.
For indoor or short-range applications where the curvature of the Earth is irrelevant, stick with Euclidean distance. There is no reason to pull in spherical math for points that are 10 meters apart inside a warehouse. The extra trigonometric operations are wasted cycles. If you are working in a projected coordinate system like UTM, the distances are already in meters and you can use the straight-line formula directly. The projection handles the curvature transformation for you. I spent an afternoon debugging a discrepancy between two systems that were both "correct" but operating in different coordinate reference systems. One was WGS84 geographic and the other was a local UTM zone. Converting both to the same CRS before computing distance solved it immediately. If you want a ready-made solution rather than rolling your own, the geographiclib Python package is well-maintained and handles the Vincenty edge cases internally. It installs with pip install geographiclib and gives you geodesic distances that are accurate to within a millimeter over the entire Earth surface. The tradeoff is a slightly larger dependency footprint, which may not matter for most applications but could be relevant if you are shipping an embedded binary or a minimal container image. The main failure modes to watch for are: coordinate order mismatches, precision loss at very close ranges, non-convergence on antipodal points with Vincenty, and forgetting to convert degrees to radians. All of these produce silent wrong answers rather than exceptions, so unit tests with known values are not optional. Put a couple of hardcoded reference pairs in your test suite and run them on every deployment. It takes about ten seconds to add and has saved me from shipping incorrect distance calculations at least twice.
