Очень медленная интерполяция с использованием `scipy.interpolate.griddata`

Я испытываю мучительно медленную работу scipy.interpolate.griddata при попытке интерполировать «почти» регулярно привязанные к сетке данные для сопоставления координат, чтобы как карта, так и данные могли быть построены с помощью matplotlib.pyplot.imshow потому что matplotlib.pyplot.pcolormesh занимает слишком много времени и, между прочим, плохо себя ведет с alpha .

Лучше всего показать пример (входные файлы можно скачать здесь ):

 import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import griddata map_extent = (34.4, 36.2, 30.6, 33.4) # data corners: lon = np.array([[34.5, 34.83806236], [35.74547079, 36.1173923]]) lat = np.array([[30.8, 33.29936152], [30.67890411, 33.17826563]]) # load saved files topo = np.load('topo.npy') lons = np.load('lons.npy') lats = np.load('lats.npy') data = np.load('data.npy') # get max res of data dlon = abs(np.array(np.gradient(lons))).max() dlat = abs(np.array(np.gradient(lats))).max() # interpolate the data to the extent of the map loni,lati = np.meshgrid(np.arange(map_extent[0], map_extent[1]+dlon, dlon), np.arange(map_extent[2], map_extent[3]+dlat, dlat)) zi = griddata((lons.flatten(),lats.flatten()), data.flatten(), (loni,lati), method='linear') 

Заговор:

 fig, (ax1,ax2) = plt.subplots(1,2) ax1.axis(map_extent) ax1.imshow(topo,extent=extent,cmap='Greys') ax2.axis(map_extent) ax2.imshow(topo,extent=extent,cmap='Greys') ax1.imshow(zi, vmax=0.1, extent=extent, alpha=0.5, origin='lower') ax1.plot(lon[0],lat[0], '--k', lw=3, zorder=10) ax1.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10) ax1.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10) ax1.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10) ax2.pcolormesh(lons,lats,data, alpha=0.5) ax2.plot(lon[0],lat[0], '--k', lw=3, zorder=10) ax2.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10) ax2.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10) ax2.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10) 

Результат:

введите описание изображения здесь

Обратите внимание , что это невозможно сделать, просто повернув данные с аффинными преобразованиями.

griddata занимает более 80 секунд за звонок с моими реальными данными, а pcolormesh занимает еще больше времени (более 2 минут!). Я посмотрел на ответ Джами и здесь ответ Джо Кингтона, но я не могу понять, как заставить его работать на меня.

Все мои наборы данных имеют точно такие же lons , что и в основном, мне нужно сопоставить их один раз с координатами карты и применить одно и то же преобразование к самим данным. Вопрос в том, как мне это сделать?

После долгого времени scipy.interpolate.griddata с scipy.interpolate.griddata медленной работой scipy.interpolate.griddata я решил griddata в пользу преобразования изображения с помощью OpenCV . В частности, Перспективная трансформация .

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

 import cv2 new_data = data.T[::-1] # calculate the pixel coordinates of the # computational domain corners in the data array w,e,s,n = map_extent dx = float(ew)/new_data.shape[1] dy = float(ns)/new_data.shape[0] x = (lon.ravel()-w)/dx y = (n-lat.ravel())/dy computational_domain_corners = np.float32(zip(x,y)) data_array_corners = np.float32([[0,new_data.shape[0]], [0,0], [new_data.shape[1],new_data.shape[0]], [new_data.shape[1],0]]) # Compute the transformation matrix which places # the corners of the data array at the corners of # the computational domain in data array pixel coordinates tranformation_matrix = cv2.getPerspectiveTransform(data_array_corners, computational_domain_corners) # Make the transformation making the final array the same shape # as the data array, cubic interpolate the data placing NaN's # outside the new array geometry mapped_data = cv2.warpPerspective(new_data,tranformation_matrix, (new_data.shape[1],new_data.shape[0]), flags=2, borderMode=0, borderValue=np.nan) 

Единственный недостаток, который я вижу в этом решении, – это небольшое смещение в данных, как показано неперекрывающимися контурами в прикрепленном изображении. (возможно, более точные) контуры данных в черном и warpПереводные контуры данных в «jet» colorscale.

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

Кто-то (а не я …) должен найти способ улучшить производительность griddata 🙂 Наслаждайтесь!

введите описание изображения здесь

Я использовал numpy ndimage.map_coordinates . Это здорово!

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html

скопировано из приведенной выше ссылки:

scipy.ndimage.interpolation.map_coordinates(input, coordinates, output=None, order=3, mode='constant', cval=0.0, prefilter=True)

Сопоставьте входной массив с новыми координатами путем интерполяции.

Массив координат используется для поиска для каждой точки вывода соответствующих координат на входе. Величина входа в этих координатах определяется сплайновой интерполяцией запрошенного порядка.

Форма вывода выводится из формы массива координат, отбрасывая первую ось. Значениями массива вдоль первой оси являются координаты во входном массиве, на котором найдено выходное значение.

  from scipy import ndimage a = np.arange(12.).reshape((4, 3)) a array([[ 0., 1., 2.], [ 3., 4., 5.], [ 6., 7., 8.], [ 9., 10., 11.]]) ndimage.map_coordinates(a, [[0.5, 2], [0.5, 1]], order=1) [ 2. 7.]