Greg Troxel said: > Or go full Pythagoras > (geodesic_distance)^2 + (h_start-h_end)^2)^1/2 That's the approach taken in the proj_lpz_dist function in the PROJ 4D API: https://github.com/OSGeo/PROJ/blob/298f7f6f967bb62fa8abf6a243dce46ae05cd474/src/4D_api.cpp#L148-L153