下面是答案中的Swift实现

func degreesToRadians(degrees: Double) -> Double {
    return degrees * Double.pi / 180
}

func distanceInKmBetweenEarthCoordinates(lat1: Double, lon1: Double, lat2: Double, lon2: Double) -> Double {

    let earthRadiusKm: Double = 6371

    let dLat = degreesToRadians(degrees: lat2 - lat1)
    let dLon = degreesToRadians(degrees: lon2 - lon1)

    let lat1 = degreesToRadians(degrees: lat1)
    let lat2 = degreesToRadians(degrees: lat2)

    let a = sin(dLat/2) * sin(dLat/2) +
    sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2)
    let c = 2 * atan2(sqrt(a), sqrt(1 - a))
    return earthRadiusKm * c
}

下面是Kotlin的一个变种:

import kotlin.math.*

class HaversineAlgorithm {

    companion object {
        private const val MEAN_EARTH_RADIUS = 6371.008
        private const val D2R = Math.PI / 180.0
    }

    private fun haversineInKm(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double {
        val lonDiff = (lon2 - lon1) * D2R
        val latDiff = (lat2 - lat1) * D2R
        val latSin = sin(latDiff / 2.0)
        val lonSin = sin(lonDiff / 2.0)
        val a = latSin * latSin + (cos(lat1 * D2R) * cos(lat2 * D2R) * lonSin * lonSin)
        val c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a))
        return MEAN_EARTH_RADIUS * c
    }
}

对于java

public static double degreesToRadians(double degrees) {
    return degrees * Math.PI / 180;
}

public static double distanceInKmBetweenEarthCoordinates(Location location1, Location location2) {
    double earthRadiusKm = 6371;

    double dLat = degreesToRadians(location2.getLatitude()-location1.getLatitude());
    double dLon = degreesToRadians(location2.getLongitude()-location1.getLongitude());

    double lat1 = degreesToRadians(location1.getLatitude());
    double lat2 = degreesToRadians(location2.getLatitude());

    double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
            Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2);
    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    return earthRadiusKm * c;
}

一个T-SQL函数，我用来根据中心的距离选择记录

Create Function  [dbo].[DistanceInMiles] 
 (  @fromLatitude float ,
    @fromLongitude float ,
    @toLatitude float, 
    @toLongitude float
  )
   returns float
AS 
BEGIN
declare @distance float

select @distance = cast((3963 * ACOS(round(COS(RADIANS(90-@fromLatitude))*COS(RADIANS(90-@toLatitude))+ 
SIN(RADIANS(90-@fromLatitude))*SIN(RADIANS(90-@toLatitude))*COS(RADIANS(@fromLongitude-@toLongitude)),15)) 
)as float) 
  return  round(@distance,1)
END

这取决于你需要它有多准确。如果你需要精确到毫米的精度，最好看看使用椭球的算法，而不是球体，比如Vincenty的算法。

如果你需要更准确的数据，可以看看这个。

Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth. The first (direct) method computes the location of a point which is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (0.020″) on the Earth ellipsoid.