[FEA] Add Vincenty Distance to cuSpatial

[MurrayData]

Vincenty Distance is an iterative formula to calculate the shortest distance between two points on a spheroid. It is more accurate than simple great circle calculations, such as Haversine, which uses an average radius for the Earth.

Vincenty takes into the account the flattening of the Earth as an oblate spheroid to calcualte to an accuracy of 0.5 mm on the WGS84 ellipsoid.

While differences between Haversine and Vincenty are small over short distances, noticeable errors arise over longer distances, particularly if calculating aircraft trajectories.

I suggest adding Vincenty Distance to cuSpatial.

This is my implementation in CUDA JIT

@cuda.jit(device=True)
def vincenty_distance(lon_1, lat_1, lon_2, lat_2,a,f):
    lon_1 = lon_1 * math.pi / 180.0 
    lon_2 = lon_2 * math.pi / 180.0 
    lat_1 = lat_1 * math.pi / 180.0 
    lat_2 = lat_2 * math.pi / 180.0
    b = (1 - f) * a                  # semi-minor axis
    u_1 = math.atan((1 - f) * math.tan(lat_1))
    u_2 = math.atan((1 - f) * math.tan(lat_2))
    L = lon_2 - lon_1
    Lambda = L                       # set initial value of lambda to L
    Lambda_2 = L + 1.0
    sin_u1 = math.sin(u_1)
    cos_u1 = math.cos(u_1)
    sin_u2 = math.sin(u_2)
    cos_u2 = math.cos(u_2)
    iters = 0
    while iters == 0 or (iters < 200 and abs(Lambda_2 - Lambda) > 1E-13):
        iters += 1
        cos_lambda = math.cos(Lambda)
        sin_lambda = math.sin(Lambda)
        sin_sigma = math.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 = math.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_2 = Lambda
        Lambda = L + (1 - C) * f * sin_alpha * (sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m**2)))
    u_sq = cos_sq_alpha * ((a**2 - b**2) / b**2)
    A = 1 + (u_sq / 16384.0)*(4096.0 + u_sq * (-768.0 + u_sq * (320.0 - 175.0 * u_sq)))
    B = (u_sq / 1024.0) * (256.0 + u_sq * (-128.0 + u_sq * (74.0 - 47.0 * u_sq)))
    delta_sig = B * sin_sigma * (cos2_sigma_m + 0.25 * B * (cos_sigma * (-1 + 2 * cos2_sigma_m**2) - (1.0 / 6.0) * B * cos2_sigma_m * (-3.0 + 4.0 * sin_sigma**2) * (-3.0 + 4.0 * cos2_sigma_m**2)))
    return b * A * (sigma - delta_sig)                 # output distance in meters

where a = semi-major axis in metres and f = flattening

For the WGS84 ellipsoid, these values are:

a = 6378137.0
f = 0.0033528106647474805

I have tested the above against online converters provided by mapping agencies.

The sample below shows the difference between Vincenty (v_dist) against Haversine (h_dist) in km on the NY Taxi dataset:

Start_Lon	Start_Lat	End_Lon	End_Lat	h_dist	v_dist	diff
-73.991957	40.721567	-73.993803	40.695922	2.855836	2.852103	-0.003733
-73.982102	40.736290	-73.955850	40.768030	4.164867	4.163947	-0.000920
-74.002587	40.739748	-73.869983	40.770225	11.672168	11.698154	0.025987
-73.974267	40.790955	-73.996558	40.731849	6.835177	6.828220	-0.006957
-74.001580	40.719382	-74.008378	40.720350	0.582929	0.584338	0.001409
-73.989806	40.735006	-73.985021	40.724494	1.236468	1.235350	-0.001117
-73.984050	40.743544	-73.980260	40.748926	0.678293	0.677985	-0.000309
-73.992635	40.748362	-73.995585	40.728307	2.243822	2.240981	-0.002841
-73.969690	40.749244	-73.990413	40.751082	1.757570	1.761962	0.004392
-73.955173	40.783044	-73.958598	40.774822	0.958650	0.957733	-0.000917
-73.986824	40.750893	-73.984118	40.751437	0.235832	0.236374	0.000542
-74.006100	40.748432	-73.978437	40.762481	2.805283	2.809085	0.003802
-73.983339	40.744782	-73.981160	40.720835	2.669107	2.665647	-0.003460

[thomcom]

Hi @MurrayData ! We were wondering if you could provide a sort of ETA on this feature, if you're still interested? Thanks!

[github-actions]

This issue has been marked stale due to no recent activity in the past 30d. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed. This issue will be marked rotten if there is no activity in the next 60d.

[github-actions]

This issue has been marked rotten due to no recent activity in the past 90d. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed.

[MurrayData]

I struggled with the C++ element of the coding as I'm very rusty in it. I have it coded in Python.

[github-actions]

This issue has been labeled inactive-30d due to no recent activity in the past 30 days. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed. This issue will be labeled inactive-90d if there is no activity in the next 60 days.

[github-actions]

This issue has been labeled inactive-90d due to no recent activity in the past 90 days. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed.

[harrism]

This is still worthwhile. Reopening.

[MurrayData]

This is still worthwhile. Reopening.

Happy to help @harrism @thomcom but my C++ isn't up to the standard. I have published the Python version above.

I have also written some additional associated functions including bearing in degrees between 2 lat/lon points, Haversine Manhattan and Vincenty Manhattan distances if you are interesting in adding these.

[harrism]

We are of course interested, we are just a small team with a long list of requests, plans and priorities. :) If you think those APIs would be useful to have, it would be great if you can add an issue for each one. Links to similar functions in existing (Python or otherwise) libraries would be helpful. Thanks!