У меня есть растр, с которым я бы хотел провести точечную интерполяцию. Вот где я нахожусь:
from osgeo import gdal
from numpy import array
# Read raster
source = gdal.Open('my_raster.tif')
nx, ny = source.RasterXSize, source.RasterYSize
gt = source.GetGeoTransform()
band_array = source.GetRasterBand(1).ReadAsArray()
# Close raster
source = None
# Compute mid-point grid spacings
ax = array([gt[0] + ix*gt[1] + gt[1]/2.0 for ix in range(nx)])
ay = array([gt[3] + iy*gt[5] + gt[5]/2.0 for iy in range(ny)])
До сих пор я пробовал функцию SciPy interp2d :
from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
однако я получаю ошибку памяти в моей 32-битной системе Windows с растром 317 × 301:
Traceback (most recent call last):
File "<interactive input>", line 1, in <module>
File "C:\Python25\Lib\site-packages\scipy\interpolate\interpolate.py", line 125, in __init__
self.tck = fitpack.bisplrep(self.x, self.y, self.z, kx=kx, ky=ky, s=0.)
File "C:\Python25\Lib\site-packages\scipy\interpolate\fitpack.py", line 873, in bisplrep
tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
MemoryError
Я признаю, я не очень уверен в этой функции SciPy, так как параметры bounds_error
или fill_value
не работают так, как задокументировано. Я не понимаю, почему у меня должна быть ошибка памяти, поскольку мой растр 317 × 301, и билинейный алгоритм не должен быть сложным.
Кто-нибудь сталкивался с хорошим алгоритмом билинейной интерполяции, предпочтительно на Python, возможно, с учетом NumPy? Любые советы или советы?
(Примечание: алгоритм интерполяции ближайшего соседа прост:
from numpy import argmin, NAN
def nearest_neighbor(px, py, no_data=NAN):
'''Nearest Neighbor point at (px, py) on band_array
example: nearest_neighbor(2790501.920, 6338905.159)'''
ix = int(round((px - (gt[0] + gt[1]/2.0))/gt[1]))
iy = int(round((py - (gt[3] + gt[5]/2.0))/gt[5]))
if (ix < 0) or (iy < 0) or (ix > nx - 1) or (iy > ny - 1):
return no_data
else:
return band_array[iy, ix]
... но я предпочитаю билинейные методы интерполяции)
python
algorithm
interpolation
Майк Т
источник
источник
MemoryError
потому что NumPy пытается получить доступ за пределами вашегоband_array
? Вы должны проверитьax
иay
.gt
.Ответы:
Я перевел приведенную ниже формулу (из Википедии ) на язык Python, чтобы получить следующий алгоритм, который, похоже, работает.
Обратите внимание, что результат будет возвращен с явно более высокой точностью, чем исходные данные, поскольку он классифицируется с повышением до
dtype('float64')
типа данных NumPy . Вы можете использовать возвращаемое значение с,.astype(band_array.dtype)
чтобы сделать тип выходных данных таким же, как у входного массива.источник
Я попробовал это локально с похожими результатами, но я на 64-битной платформе, поэтому он не достиг предела памяти. Возможно, вместо этого попробуйте интерполировать небольшие фрагменты массива за раз, как в этом примере .
Вы также можете просто сделать это с GDAL, из командной строки это будет:
Чтобы сделать эквивалентную операцию в Python, используйте ReprojectImage () :
источник
ReprojectImage
технику GDAL .У меня была точная проблема в прошлом, и я никогда не решал ее с помощью interpolate.interp2d. Я имел успех, используя scipy.ndimage.map_coordinates . Попробуйте следующее:
scipy.ndimage.map_coordinates (band_array, [ax, ay]], order = 1)
Кажется, это дает тот же результат, что и билинейный.
источник
scipy.interpolate.interp2d () отлично работает с более современным scipy. Я думаю, что старые версии предполагают нерегулярные сетки и не используют преимущества обычных сеток. Я получаю ту же ошибку, что и вы со Сципи. версия = 0.11.0, но на scipy. Версия = 0.14.0, она успешно работает на некоторых моделях 1600x1600.
Спасибо за подсказки в вашем вопросе.
источник