matplotlib.mlab.griddata очень медленный и возвращает массив nan при вводе правильных данных
Я пытаюсь сопоставить набор данных с нерегулярной сеткой (необработанные спутниковые данные) с соответствующими широтами и долготами с набором широт и долгот с регулярной сеткой, заданным basemap.makegrid(), я использую matplotlib.mlab.griddata с mpl_toolkits.natgrid установлены. Ниже приведен список переменных, используемых в качестве вывода whos в ipython и немного статистики по переменным:
Variable Type Data/Info
-------------------------------
datalat ndarray 666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
datalon ndarray 666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
gridlat ndarray 1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
gridlon ndarray 1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
var ndarray 666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
In [11]: var.min()
Out[11]: -30.0
In [12]: var.max()
Out[12]: 30.0
In [13]: datalat.min()
Out[13]: 27.339874
In [14]: datalat.max()
Out[14]: 47.05302
In [15]: datalon.min()
Out[15]: -137.55658
In [16]: datalon.max()
Out[16]: -108.41629
In [17]: gridlat.min()
Out[17]: 30.394031556984299
In [18]: gridlat.max()
Out[18]: 44.237140350357713
In [19]: gridlon.min()
Out[19]: -136.17646180595321
In [20]: gridlon.max()
Out[20]: -113.82353819404671
datalat а также datalon координаты данных оригинала
gridlat а также gridlon координаты для интерполяции
var содержит фактические данные
Используя эти переменные, когда я звоню griddata(datalon, datalat, var, gridlon, gridlat) это заняло до 20 минут, и возвращает массив nan, При просмотре данных широта и долгота кажутся правильными: исходные координаты перекрывают часть новой области, а несколько точек данных лежат за пределами новой области. У кого-нибудь есть предложения? Значения Nan предполагают, что я делаю что-то глупое...
4 ответа
Похоже, mlab.griddata рутина может вводить дополнительные ограничения на ваши выходные данные, которые могут не быть необходимыми. Хотя входные местоположения могут быть любыми, выходные местоположения должны быть регулярной сеткой - так как ваш пример находится в широтном / долгом пространстве, ваш выбор проекции карты может нарушить это (т.е. обычная сетка в x/y не является регулярной сеткой в лат. / долгота).
Вы можете попробовать interpolate.griddata подпрограмма из SciPy в качестве альтернативы - вам нужно будет объединить ваши переменные местоположения в один массив, так как сигнатура вызова отличается: что-то вроде
import scipy.interpolate
data_locations = np.vstack(datalon.ravel(), datalat.ravel()).T
grid_locations = np.vstack(gridlon.ravel(), gridlat.ravel()).T
grid_data = scipy.interpolate.griddata(data_locations, val.ravel(),
grid_locations, method='nearest')
для интерполяции ближайшего соседа. Это получает местоположения в массиве с 2 столбцами, соответствующими вашим 2 измерениям. Вы также можете выполнить интерполяцию в преобразованном пространстве вашей картографической проекции.
Скорее всего, griddata слишком сложна. Он предназначен для работы со случайно выбранными данными. Ваши данные почти наверняка регулярно отбираются - но не в той же самой сетке, что и целевая выходная сетка.
Посмотрите на гораздо более простой подход, такой как аффинное преобразование или серия аффинных преобразований на маленьких чипах, если топология или кривизна Земли влияют на ваши результаты.
Есть некоторые готовые решения, которые могут помочь. GDAL является хорошим примером.
Кроме того, этот тип проблемы часто обсуждается в ГИС. Увидеть:
https://gis.stackexchange.com/questions/10430/changing-image-projection-using-python
Если ваши данные находятся в сетке, так что данные указывают на точку (datalon[i], datalat[j]) в data[i,j]тогда вы можете использовать scipy.interpolate.RectBivariateSpline вместо griddata, Однако некоторые специфические для географии библиотеки могут предлагать больше функциональных возможностей.
Если вы используете pclormesh, вам не нужно выполнять какую-либо интерполяцию. pcolormesh с удовольствием примет структуру данных, как вы указали здесь:
from mpl_toolkits.basemap import Basemap
m = Basemap(-----)
x,y = m(datalon, datalat)
m.pcolormesh(x,y,var)
plt.show()
Пожалуйста, используйте это и скажите мне, работает ли это или нет.
Однако в pcolormesh возникает некоторая проблема, когда данные орбиты перекрываются. Пожалуйста, обратитесь к моему вопросу, вы можете найти что-то полезное.