Попробуйте Python, за исключением того, что не работает с np.arccos

Question

Это мой первый вопрос здесь, поэтому укажите мне, если я что-то не так делаю.

Я реализовывал уравнение для вычисления сферического угла (https://en.wikipedia.org/wiki/Solution_of_triangles#Solving_spherical_triangles). Как math.acos и np.arccos определены только для [-1,1], я должен закодировать оператор try / исключением. Он хорошо работал с math.acos, но я продолжаю получать

RuntimeWarning: недопустимое значение, встречающееся в arccos

когда я использую np.arccos. Поскольку я хочу вычислить сферический угол для очень длинного файла, мне нужен np.arccos, чтобы я мог воспользоваться преимуществами векторизации. Любые идеи о том, почему это происходит или как я могу улучшитькод?

Вот моя попытка:

import numpy as np
import math

def spherical_to_cartesian(latr, lonr):
    """Convert a point given latitude and longitude in radians to
    3D cartesian coordinates, assuming a sphere radius of one."""
    return np.array((
        np.cos(latr) * np.cos(lonr),
        np.cos(latr) * np.sin(lonr),
        np.sin(latr)
    ))

def angle_between_vectors(u, v):
    'Return the angle between two vectors in any dimension space in radians'
    return np.arccos(np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v)))

def angle_3pts(p1, p2, p3):
    'Convert the points to numpy latitude/longitude radians space'
    p1_rad = np.radians(np.array(p1))
    p2_rad = np.radians(np.array(p2))
    p3_rad = np.radians(np.array(p3))

    'The points in 3D cartesian coordinates'
    p1_cartesian = spherical_to_cartesian(*p1_rad)
    p2_cartesian = spherical_to_cartesian(*p2_rad)
    p3_cartesian = spherical_to_cartesian(*p3_rad)

    'The angle between vectors'
    ang_12 = angle_between_vectors(p1_cartesian, p2_cartesian)
    ang_23 = angle_between_vectors(p2_cartesian, p3_cartesian)
    ang_31 = angle_between_vectors(p3_cartesian, p1_cartesian)
    print(ang_12, ang_23, ang_31)

    '''Calculate spherical angle: 
       https://en.wikipedia.org/wiki/Solution_of_triangles#Solving_spherical_triangles'''
    print((np.cos(ang_31) - np.cos(ang_12)*np.cos(ang_23))/(np.sin(ang_12)*np.sin(ang_23)))

    try:
        angle_3d = np.degrees(np.arccos((np.cos(ang_31) - np.cos(ang_12)*np.cos(ang_23))/(np.sin(ang_12)*np.sin(ang_23))))
    except:
        angle_3d = 180.

    return angle_3d


point1 = (0, 0)
point2 = (0, 10)
point3 = (0, 20)

ang = angle_3pts(point1, point2, point3)

Когда я использую math.acos внутри уравнения try, он работает хорошо, но np.arccos выдает ошибку

Answer 1

arccos является ufunc и принимает параметры where и out:

In [226]: x = np.linspace(-2,2,20)
In [227]: np.arccos(x, where=(abs(x)<1), out=np.full_like(x,np.pi))
Out[227]: 
array([3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265,
       2.8157097 , 2.39918372, 2.12505816, 1.89208492, 1.67625485,
       1.4653378 , 1.24950774, 1.01653449, 0.74240893, 0.32588296,
       3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265])

Это позволяет указать, какие входные значения следует пропустить, а какие использовать вместо.

По сравнению с проверенной math версией:

def foo(x):
    try:
        return math.acos(x)
    except ValueError:
        return math.pi
In [239]: np.array([foo(i) for i in x])
Out[239]: 
array([3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265,
       2.8157097 , 2.39918372, 2.12505816, 1.89208492, 1.67625485,
       1.4653378 , 1.24950774, 1.01653449, 0.74240893, 0.32588296,
       3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265])