LoginSignup
0
0

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-02-04
動作環境
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
の実装よりも精度が高いのだろうか。。。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0