Numpy > MATLABのsph2cart()の実装 > sph2cart_180204.py > 直交座標から球座標への変換 v0.1..v0.3


動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6

MATLABのsph2cart()

球座標から直交座標への変換sph2cart()

geometry > x,y,zベクトルから天頂角と方位角を計算する > 式の整理 | Numpy実装 polarAzimuthCalc_171209.py > v0.1, v0.2 | sys.float_info_epsilonのpitfalls
に記載の式(1)から(3)を使う。

ただし、天頂角(zenith angle)とelevationはpi/2で反転するためその部分を考慮した実装にする。

code v0.1 > 要素ずつの変換

test_sph2cart_180204.py
import numpy as np

azs = [0.7854, 0.7854, -0.7854, -0.7854, 2.3562, 2.3562, -2.3562, -2.3562]
els = [0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155]
rs = [1.7321, 1.7321, 1.7321, 1.7321, 1.7321, 1.7321, 1.7321, 1.7321]


def sph2cart(azimuth, elevation, radius):
    ax = np.cos(azimuth) * np.sin(np.pi/2 - elevation) * radius
    ay = np.sin(azimuth) * np.sin(np.pi/2 - elevation) * radius
    az = np.cos(np.pi/2 - elevation) * radius
    return ax, ay, az


for idx, (azi, elv, rad) in enumerate(zip(azs, els, rs)):
    ax, ay, az = sph2cart(azi, elv, rad)
    print("%d %.3f %.3f %.3f" % (idx, ax, ay, az))

run
$ python3 test_sph2cart_180204.py 
0 1.000 1.000 1.000
1 1.000 1.000 -1.000
2 1.000 -1.000 1.000
3 1.000 -1.000 -1.000
4 -1.000 1.000 1.000
5 -1.000 1.000 -1.000
6 -1.000 -1.000 1.000
7 -1.000 -1.000 -1.000

MATLABの上記ドキュメントと同じ結果になった。

code v0.2 > 全要素の一括変換

実際に使いたいのは全要素を一括変換する処理。

test_sph2cart_180204.py
import numpy as np

'''
v0.2 Feb.04, 2018
    - convert all elements at one
        + sph2cart() processes a set
v0.1 Feb.04, 2018
    - convert each element
        + add sph2cart()
'''

azs = [0.7854, 0.7854, -0.7854, -0.7854, 2.3562, 2.3562, -2.3562, -2.3562]
els = [0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155]
rs = [1.7321] * len(azs)
azs = np.array(azs)
els = np.array(els)
rs = np.array(rs)


def sph2cart(azimuth, elevation, radius):
    axs = np.cos(azimuth) * np.sin(np.pi/2 - elevation) * radius
    ays = np.sin(azimuth) * np.sin(np.pi/2 - elevation) * radius
    azs = np.cos(np.pi/2 - elevation) * radius
    return axs, ays, azs


axs, ays, azs = sph2cart(azs, els, rs)
np.set_printoptions(precision=5)  # for debug
print(axs)
print(ays)
print(azs)

run
$ python3 test_sph2cart_180204.py 
[ 1.00001  1.00001  1.00001  1.00001 -1.00002 -1.00002 -1.00002 -1.00002]
[ 1.00002  1.00002 -1.00002 -1.00002  1.00001  1.00001 -1.00001 -1.00001]
[ 1.00006 -1.00006  1.00006 -1.00006  1.00006 -1.00006  1.00006 -1.00006]

code v0.3 > import対応

sph2cart_180204.py
import numpy as np

'''
v0.3 Feb.04, 2018
    - can be used by "import"
    - renamed to [sph2cart_180204]
v0.2 Feb.04, 2018
    - convert all elements at one
        + sph2cart() processes a set
v0.1 Feb.04, 2018
    - convert each element
        + add sph2cart()
'''

# testd on Python 3.5.2
# coding rule:PEP8


def sph2cart(azimuth, elevation, radius):
    axs = np.cos(azimuth) * np.sin(np.pi/2 - elevation) * radius
    ays = np.sin(azimuth) * np.sin(np.pi/2 - elevation) * radius
    azs = np.cos(np.pi/2 - elevation) * radius
    return axs, ays, azs


def test_cubic_8vertices():
    AZS = [0.7854, 0.7854, -0.7854, -0.7854, 2.3562, 2.3562, -2.3562, -2.3562]
    ELS = [0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155, 0.6155, -0.6155]
    RS = [1.7321] * len(AZS)

    azs = np.array(AZS)
    els = np.array(ELS)
    rs = np.array(RS)

    axs, ays, azs = sph2cart(azs, els, rs)
    np.set_printoptions(precision=5)  # for debug
    print(axs)
    print(ays)
    print(azs)


if __name__ == '__main__':
    test_cubic_8vertices()

関連

直交座標から球座標への変換の高速版は下記で議論されている。

geometry > x,y,zベクトルから天頂角と方位角を計算する > 式の整理 | Numpy実装 polarAzimuthCalc_171209.py > v0.1, v0.2 | sys.float_info_epsilonのpitfalls
の実装よりも精度が高いのだろうか。。。