Faster alternative to SkyCoord.separation?

Question

I am using SkyCoord.separation to find the pixels around a given coordinate from an array. I have noticed the separation function takes too long. I have to repeat the calculation on a huge data set and using this function does not seem practical. Is there a faster alternative to SkyCoord.separation which handles arrays? The other similar functions do not seem to take array inputs and can only calculate separation between two sets of coordinates.

For example: I have longitude and latitude arrays (size 50,331,648). I need to find the separation of each row with the rest of the array. So I run a loop 50,331,648 times.

Any suggestions would be of great help! Thank you

update: Using cartesian coordinates and dot product for the angle now. It is around 7 times faster.

Answer 1

Here is a pretty fast implementation that uses numexpr for extra speed. You can also do this with plain numpy but it isn't as fast.

import numexpr

def sphere_dist(ra1, dec1, ra2, dec2):
    """
    Haversine formula for angular distance on a sphere: more stable at poles.
    This version uses arctan instead of arcsin and thus does better with sign
    conventions.  This uses numexpr to speed expression evaluation by a factor
    of 2 to 3.

    :param ra1: first RA (deg)
    :param dec1: first Dec (deg)
    :param ra2: second RA (deg)
    :param dec2: second Dec (deg)

    :returns: angular separation distance (deg)
    """

    ra1 = np.radians(ra1).astype(np.float64)
    ra2 = np.radians(ra2).astype(np.float64)
    dec1 = np.radians(dec1).astype(np.float64)
    dec2 = np.radians(dec2).astype(np.float64)

    numerator = numexpr.evaluate('sin((dec2 - dec1) / 2) ** 2 + '
                                 'cos(dec1) * cos(dec2) * sin((ra2 - ra1) / 2) ** 2')

    dists = numexpr.evaluate('2 * arctan2(numerator ** 0.5, (1 - numerator) ** 0.5)')
    return np.degrees(dists)