Эффективно находите ближайшие точки для отслеживания в пространстве и времени на данных с координатной сеткой

Общая/упрощенная версия

Дан список точек трека, определяемый тремя одномерными массивами (lats, lons и dtime, все одинаковой длины) и трехмерным массивом с координатной сеткой rr (определяемым двумерными координатными массивами lat_radar, lon_radar и одномерным временным массивом time_radar) Я хочу извлечь все значения сетки в rr, где координаты (включая широту, долготу и время) ближе всего к трем одномерным массивам.

Мне удалось использовать cKDTree для выбора точек в пространстве, но я не знаю, как обобщить решение для пространства и времени вместе. Прямо сейчас я должен делать выборку по времени отдельно, и это делает код довольно громоздким и трудным для чтения.

подробнее об этой проблеме см. далее

Расширенная версия

Я пытаюсь разработать приложение, которое использует данные об осадках, полученные из композитов метеорологических радаров, для прогнозирования осадков вдоль трассы. Большинство приложений обычно предсказывают осадки в точке без учета движения точки во времени.

Идея состоит в том, чтобы, учитывая точки, идентифицирующие трек в пространстве и времени, найти ближайшие точки сетки по данным радара, чтобы получить оценку осадков по треку (см. График). Конечной целью будет смещение времени начала, чтобы определить лучшее время для выезда, чтобы избежать дождя.

Я только что оптимизировал свой предыдущий алгоритм, в котором использовались простые циклы, для использования cKDTree из scipy. Время выполнения сократилось с 30 до 380 мс :). Однако я думаю, что код все еще можно оптимизировать. Вот моя попытка.

В качестве входных данных мы имеем

• lons, lats: координаты трека в виде N-мерных массивов
• dtime: timedelta Т-мерный массив, содержащий время, прошедшее на треке
• lon_radar, lat_radar: матрицы M x P, содержащие координаты данных радара.
• dtime_radar: timedelta Q-мерный массив, содержащий радиолокационный прогноз
• rr: Массив M x P X Q, содержащий радиолокационный прогноз на каждом временном шаге

Сначала найдите точки сетки, ближайшие к траектории, используя cKDTree:

combined_x_y_arrays = np.dstack([lon_radar.ravel(), 
                                 lat_radar.ravel()])[0]
points_list = list(np.vstack([lons, lats]).T)

def do_kdtree(combined_x_y_arrays, points):
    mytree = cKDTree(combined_x_y_arrays)
    dist, indexes = mytree.query(points)
    return indexes

results = do_kdtree(combined_x_y_arrays, points_list)
# As we have many duplicates, since the itinerary has a much higher resolution than the radar,
# we only select the unique points 
inds_itinerary = np.unique(results)
lon_lat_itinerary = combined_x_y_arrays[inds_itinerary]

затем найдите ближайшие точки на дорожке, чтобы подмножить ее. Не имеет смысла иметь разрешение трека 10 м, если радар имеет точки сетки только через каждый км.

combined_x_y_arrays = np.vstack([lons, lats]).T
points_list = list(lon_lat_itinerary)

results = do_kdtree(combined_x_y_arrays, points_list)

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

dtime_itinerary = dtime[results]
# find indices of these dtimes in radar dtime 
inds_dtime_radar = np.abs(np.subtract.outer(dtime_radar, dtime_itinerary)).argmin(0)

Теперь у нас есть все, что нам нужно, чтобы найти осадки, поэтому нам нужен только один последний цикл. Я также зацикливаюсь на shifts, чтобы получить прогноз с разным временем начала.

shifts = (1, 3, 5, 7, 9)

rain = np.empty(shape=(len(shifts), len(inds_itinerary)))
for i, shift in enumerate(shifts):
    temp = []
    for i_time, i_space in zip(inds_dtime_radar, inds_itinerary):
        temp.append(rr[i_time+shift].ravel()[i_space])
    rain[i, :] = temp

В частности, я хотел бы найти способ объединить поиск по времени с поиском ближайших точек по широте и долготе.