Numpy's arctan2(y, x) will compute the counterclockwise angle (a value in radians between -π and π) between the origin and the point (x, y).

You could do this for your points A and B, then subtract the second angle from the first to get the signed clockwise angular difference. This difference will be between -2π and 2π, so in order to get a positive angle between 0 and 2π you could then take the modulo against 2π. Finally you can convert radians to degrees using np.rad2deg.

import numpy as np

def angle_between(p1, p2):
    ang1 = np.arctan2(*p1[::-1])
    ang2 = np.arctan2(*p2[::-1])
    return np.rad2deg((ang1 - ang2) % (2 * np.pi))

For example:

A = (1, 0)
B = (1, -1)

print(angle_between(A, B))
# 45.

print(angle_between(B, A))
# 315.

If you don't want to use numpy, you could use math.atan2 in place of np.arctan2, and use math.degrees (or just multiply by 180 / math.pi) in order to convert from radians to degrees. One advantage of the numpy version is that you can also pass two (2, ...) arrays for p1 and p2 in order to compute the angles between multiple pairs of points in a vectorized way.

Note: all of the other answers here will fail if the two vectors have either the same direction (ex, (1, 0, 0), (1, 0, 0)) or opposite directions (ex, (-1, 0, 0), (1, 0, 0)).

Here is a function which will correctly handle these cases:

import numpy as np

def unit_vector(vector):
    """ Returns the unit vector of the vector.  """
    return vector / np.linalg.norm(vector)

def angle_between(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'::

            >>> angle_between((1, 0, 0), (0, 1, 0))
            1.5707963267948966
            >>> angle_between((1, 0, 0), (1, 0, 0))
            0.0
            >>> angle_between((1, 0, 0), (-1, 0, 0))
            3.141592653589793
    """
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

Your original code is pretty close. Adomas.m's answer is not very idiomatic numpy:

import numpy as np

a = np.array([32.49, -39.96,-3.86])
b = np.array([31.39, -39.28, -4.66])
c = np.array([31.14, -38.09,-4.49])

ba = a - b
bc = c - b

cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
angle = np.arccos(cosine_angle)

print np.degrees(angle)

Since the norm is just the magnitude or modulus or whatever you call it, this code should do the trick:

from numpy import arccos, array
from numpy.linalg import norm

# Note: returns angle in radians
def theta(v, w): return arccos(v.dot(w)/(norm(v)*norm(w)))

sulfur = array([0.0, 0.0, 0.102249])
hydrogen_1 = array([0.0, 0.968059, -0.817992])
hydrogen_2 = array([0.0, -0.968059, -0.817992])
print(theta(hydrogen_1-sulfur, hydrogen_2-sulfur))

The numpy function is documented here.

Here's a way of packing the calculation into on call to einsum.

def findAnglesBetweenTwoVectors1(v1s, v2s):
    dot = np.einsum('ijk,ijk->ij',[v1s,v1s,v2s],[v2s,v1s,v2s])
    return np.arccos(dot[0,:]/(np.sqrt(dot[1,:])*np.sqrt(dot[2,:])))

for random vectors of length 10, the speed for this function and yours is basically the same. But with 10000, (corrected times)

In [62]: timeit findAnglesBetweenTwoVectors0(v1s,v2s)
100 loops, best of 3: 2.26 ms per loop

In [63]: timeit findAnglesBetweenTwoVectors1(v1s,v2s)
100 loops, best of 3: 4.24 ms per loop

Yours is clearly faster.

We could try to figure out why the more complex einsum is slower (assuming the issue is there rather than the indexing in the last line). But this result is consistent with other einsum testing that I've seen and done. einsum on simple 2D cases is as fast as np.dot, and in some cases faster. But in more complex cases (3 or more indexes) it is slower than the equivalent constructed from multiple calls to np.dot or einsum.

I tried tried factoring out the array construction, and got a speed improvement.

In [12]: V1=np.array(v1s)  # shape (3,10000,3)
In [13]: V2=np.array(v2s)
In [22]: timeit findAnglesBetweenTwoVectors2(V1,V2)
100 loops, best of 3: 2.63 ms per loop

I've corrected the times, using v1s.shape==(10000,3)

Seems to me that it works fine. Due to the fact that I dont see the character of you numpy/pandas coordinates array I cant give you exact solution

3 versions:

arctan

>>> direction = np.rad2deg(np.arctan((2-1)/(2-1)))
>>> direction
45.0

math

>>> direction = math.degrees(math.atan((2 - 1) / (2 - 1)))
>>> direction
45.0

arctan2

>>> direction = np.rad2deg(np.arctan2((2-1),(2-1)))
>>> direction
45.0

With fake data (as I dont see your data), using your function,

reference point:

>>> Xr = 10
>>> Yr = 10

Check direction towards reference point with these coordinates:

>>> Xs = np.array([1,2,3,4,5,6])
>>> Ys = np.array([1,2,3,4,5,6])

I am expecting all of them to have direction of 45 degrees

>>> direction = np.rad2deg(np.arctan2(Yr-Ys,Xr-Xs))
>>> direction
array([45., 45., 45., 45., 45., 45.])

EDIT

Based on your provided data:

It seems you were correct with the notion that there is different result.

More in here: numpy arctan2 bug or usage issues?

using np.arctan should be OK

numpy module

>>> direction = np.rad2deg(np.arctan(data.y/data.x))
>>> direction
0     45.819271
1     45.805556
2     45.795956
3     45.790496
4     45.785684
5     45.782246
6     45.771168
7     45.759343
8     45.744059
9     45.724642
10    45.712089
dtype: float64

math MODULE

>>> for i in range(10):
...     math.degrees(math.atan((data.y[i]) / (data.x[i])))
... 
45.81927143732942
45.805555571378584
45.79595585089821
45.79049607412479
45.78568399098587
45.782246280023934
45.771167675155795
45.75934340359498
45.744059202888714
45.72464168864543
>>> 

NOTE

Be careful about sign +-. Y2-Y1/X2-X1 vs Y1-Y2/X1-X2 can give different result, but in really both correct. In our case, these results are 45 deg or -135 deg. They both are correct, just one of them is clockwise