Calculate distance between Latitude/Longitude points. All these formulas are for calculations on the basis of a spherical earth (ignoring ellipsoidal effects) – which is accurate enough for most purposes…

lat_long_distance.py

from math import radians, sin, cos, atan2, sqrt, tan, atan

def haversine_distance(long1, lat1, long2, lat2, degrees=False):
  '''
  The haversine formula determines the great-circle distance 
  between two points on a sphere given their longitudes and latitudes.
  '''
  #degrees vs radians
  if degrees == True:
    long1 = radians(long1)
    lat1 = radians(lat1)
    long2 = radians(long2)
    lat2 = radians(lat2)
  #implementing haversine
  a = sin((lat2-lat1) / 2)**2 + cos(lat1) * cos(lat2) *
  sin((long2-long1) / 2)**2
  c = 2*atan2(sqrt(a), sqrt(1-a))
  distance = 6371 * c #radius of earth in kilometers
  return(distance)

def vincenty(long1, lat1, long2, lat2, maxIter=200,tol=10**-12):
        # Modified from Nathan A. Rooy's code
        #--- CONSTANTS ------------------------------------+
        
        a=6378137.0                             # radius at equator in meters (WGS-84)
        f=1/298.257223563                       # flattening of the ellipsoid (WGS-84)
        b=(1-f)*a

        u_1=atan((1-f)*tan(radians(lat1)))
        u_2=atan((1-f)*tan(radians(lat2)))

        L=radians(long2-long1)

        Lambda=L                                # set initial value of lambda to L

        sin_u1=sin(u_1)
        cos_u1=cos(u_1)
        sin_u2=sin(u_2)
        cos_u2=cos(u_2)

        #--- BEGIN ITERATIONS -----------------------------+
        for i in range(0,maxIter):
            
            cos_lambda=cos(Lambda)
            sin_lambda=sin(Lambda)
            sin_sigma=sqrt((cos_u2*sin(Lambda))**2+(cos_u1*sin_u2-sin_u1*cos_u2*cos_lambda)**2)
            cos_sigma=sin_u1*sin_u2+cos_u1*cos_u2*cos_lambda
            sigma=atan2(sin_sigma,cos_sigma)
            sin_alpha=(cos_u1*cos_u2*sin_lambda)/sin_sigma
            cos_sq_alpha=1-sin_alpha**2
            cos2_sigma_m=cos_sigma-((2*sin_u1*sin_u2)/cos_sq_alpha)
            C=(f/16)*cos_sq_alpha*(4+f*(4-3*cos_sq_alpha))
            Lambda_prev=Lambda
            Lambda=L+(1-C)*f*sin_alpha*(sigma+C*sin_sigma*(cos2_sigma_m+C*cos_sigma*(-1+2*cos2_sigma_m**2)))

            # successful convergence
            diff=abs(Lambda_prev-Lambda)
            if diff<=tol:
                break
            
        u_sq=cos_sq_alpha*((a**2-b**2)/b**2)
        A=1+(u_sq/16384)*(4096+u_sq*(-768+u_sq*(320-175*u_sq)))
        B=(u_sq/1024)*(256+u_sq*(-128+u_sq*(74-47*u_sq)))
        delta_sig=B*sin_sigma*(cos2_sigma_m+0.25*B*(cos_sigma*(-1+2*cos2_sigma_m**2)-(1/6)*B*cos2_sigma_m*(-3+4*sin_sigma**2)*(-3+4*cos2_sigma_m**2)))

        return b*A*(sigma-delta_sig)                 # output distance in meters