*** 以下の式導出は正しくないかもしれない ***

動作環境

Xeon E5-2620 v4 (8コア) x 2
32GB RAM
CentOS 6.8 (64bit)
openmpi-1.8.x86_64 とその-devel
mpich.x86_64 3.1-5.el6とその-devel
gcc version 4.4.7 (とgfortran)
NCAR Command Language Version 6.3.0
WRF v3.7.1を使用。
Python 2.6.6 (r266:84292, Aug 18 2016, 15:13:37) 
Python 3.6.0 on virtualenv

緯度緯度の２セットのデータから距離を求める。

自分が馴染みのある式は以下 (式1..3)。
参考 https://en.wikipedia.org/wiki/Spherical_coordinate_system
のCartesian coordinatesの2つ目の式。

x = rsin\theta cos\phi

y = rsin\theta sin\phi

z = rcos\theta

ここで、$r$はradius (地球の半径[km])。$\theta$は天頂角 (zenith angle)のradianを取ったもの。$\phi$は方位角 (polar angle)のradianを取ったもの。
参考 https://en.wikipedia.org/wiki/Spherical_coordinate_system
の右上の図。

上記をlatitude, longitudeにどう関連付けるか。

$\theta$の定義は[0, pi]。latitudeの定義は[pi/2, -pi/2]とする(pi/2は北極点とする)。
関係式は以下と置ける (式4)。

\theta = lat - pi/2

$\phi$の定義は[0, 2pi)として、どこを0に取るかは任意とすると、関係式を以下と置く (式5)。

\phi = lon

式1..3に式4,5を代入して、加法定理を使うと以下を得る。

x = -r cos(lat)cos(lon)

y = -r sin(lat)sin(lon)

z = rcos(lon)

http://qiita.com/s-wakaba/items/e12f2a575b6885579df7
のlatlng_to_xyzの式と正負の符号が異なるのは未消化。

2点の距離を計算するときはzip()した値同士でかけて負号は消えるので、同じ結果にはなる。
( Osaka, Tokyo ): 403.63km
( London, Tokyo ): 9571.22km

上記リンクでは「何の式か」はコード自体からは分からない。

自分が理解しやすいhaversine法を使う方が良さそう。