Эффективно рассчитать расстояния между тысячами пар координат

Question

У меня есть каталог, который я открыл в python, в котором содержится около 70000 строк данных (ra, dec dec и имя объекта) для различных объектов.У меня также есть другой список из примерно 15 000 объектов, представляющих интерес, который также присутствует в ранее упомянутом каталоге.Для каждого из этих 15 000 объектов я хотел бы видеть, имеют ли какие-либо другие объекты в большом списке из 70 000 объектов координаты ra, dec в течение 10 угловых секунд объекта.Если это окажется правдой, я просто хочу пометить объект и перейти к следующему.Однако этот процесс занимает много времени, так как расстояния вычисляются между текущим объектом интереса (из 15 000) 70 000 раз.Это займет дни!Как я могу выполнить ту же задачу более эффективно?Ниже приведен мой текущий код, где all_objects - это список всех 15 000 имен объектов, представляющих интерес, а catalog - ранее упомянутые данные таблицы для 70000 объектов.

from astropy.coordinates import SkyCoord
from astropy import units as u

for obj_name in all_objects:
    obj_ind = list(catalog['NAME']).index(obj_name)
    c1 = SkyCoord(ra=catalog['RA'][obj_ind]*u.deg, dec=catalog['DEC'][obj_ind]*u.deg, frame='fk5')
    for i in range(len(catalog['NAME'])):
        if i != obj_ind:
            # Compute distance between object and other source
            c2 = SkyCoord(ra=catalog['RA'][i]*u.deg, dec=catalog['DEC'][i]*u.deg, frame='fk5')
            sep = c1.separation(c2)
            contamination_flag = False
            if sep.arcsecond <= 10:
                contamination_flag = True
                print('CONTAMINATION FOUND')
                break

Answer 1

1 Создайте свою собственную функцию разделения

Этот шаг действительно прост, если вы посмотрите на реализацию и спросите себя: «Как я могу сделать это быстрее»

def separation(self, other):
    from . import Angle
    from .angle_utilities import angular_separation # I've put that in the code bellow so it is clearer

    if not self.is_equivalent_frame(other):
        try:
            other = other.transform_to(self, merge_attributes=False)
        except TypeError:
            raise TypeError('Can only get separation to another SkyCoord '
                            'or a coordinate frame with data')

    lon1 = self.spherical.lon
    lat1 = self.spherical.lat
    lon2 = other.spherical.lon
    lat2 = other.spherical.lat

    sdlon = np.sin(lon2 - lon1)
    cdlon = np.cos(lon2 - lon1)
    slat1 = np.sin(lat1)
    slat2 = np.sin(lat2)
    clat1 = np.cos(lat1)
    clat2 = np.cos(lat2)

    num1 = clat2 * sdlon
    num2 = clat1 * slat2 - slat1 * clat2 * cdlon
    denominator = slat1 * slat2 + clat1 * clat2 * cdlon

    return Angle(np.arctan2(np.hypot(num1, num2), denominator), unit=u.degree)

Он вычисляетмного косинусов и синусов, затем создает экземпляр Angle и конвертирует в градусы, затем вы конвертируете в угловые секунды.

Возможно, вы не захотите использовать Angle, а также тесты и преобразования в начале, а такжеделать импорт в функции, и не делать так много назначения переменных, если вам нужна производительность.

Функция разделения кажется мне немного тяжелой, она должна просто брать числа и возвращать число.

2 Использовать дерево квадратов (требуется полная перезапись вашего кода)

Тем не менее, давайте посмотрим на сложность вашего алгоритма, он сравнивает каждый элемент с каждым другим элементом, сложность равна O(n**2) (нотация Big O).Можем ли мы сделать лучше ...

ДА Вы можете использовать Quad-дерево, в наихудшем случае сложность Quad-дерева равна O (N).В основном это означает, что если вы не знакомы с Big O, это то, что для элемента 15 000 поиск будет в 1021 * раз больше, чем для элемента 1 вместо 225 000 000 раз (15 000 в квадрате).... довольно хорошее улучшение ... У Сципи есть отличная библиотека Quad Tree (я всегда использовал свою).