Compute the angles between lon / lat points p1 and p2 given in radians.
    On the unit sphere, this also corresponds to the great circle distance.
    p1 and p2 can be numpy arrays of the same length.

def angular_separation(lonp1, latp1, lonp2, latp2):
    """
    Compute the angles between lon / lat points p1 and p2 given in radians.
    On the unit sphere, this also corresponds to the great circle distance.
    p1 and p2 can be numpy arrays of the same length.
    """

    xp1, yp1, zp1 = lonlat2xyz(lonp1, latp1)
    xp2, yp2, zp2 = lonlat2xyz(lonp2, latp2)

    ## dot products to obtain angles

    angles = np.arccos((xp1 * xp2 + yp1 * yp2 + zp1 * zp2))

    ## As this is a unit sphere, angle = length

    return angles