Eu tenho uma varredura com a qual gostaria de fazer algumas interpolações de pontos. Aqui é onde estou:
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)])
Até agora, eu tentei a função interp2d do SciPy :
from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
no entanto, recebo um erro de memória no meu sistema Windows de 32 bits com uma varredura de 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
Admito que tenho confiança limitada nessa função SciPy, pois os parâmetros bounds_error
ou fill_value
não funcionam como documentado. Não vejo por que devo ter um erro de memória, pois minha varredura é 317 × 301 e o algoritmo bilinear não deve ser difícil.
Alguém encontrou um bom algoritmo de interpolação bilinear, de preferência em Python, possivelmente adaptado ao NumPy? Alguma dica ou conselho?
(Nota: o algoritmo de interpolação do vizinho mais próximo é fácil:
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]
... mas eu prefiro métodos de interpolação bilineares)
fonte
MemoryError
porque o NumPy tenta acessar além do seuband_array
? Você deve verificarax
eay
.gt
.Respostas:
Traduzi a fórmula abaixo (da Wikipedia ) para Python-speak para produzir o seguinte algoritmo, que parece funcionar.
Observe que o resultado será retornado com uma precisão aparentemente mais alta que os dados de origem, pois é classificado de acordo com o
dtype('float64')
tipo de dados de NumPy . Você pode usar o valor de retorno com.astype(band_array.dtype)
para fazer com que o tipo de dados de saída seja o mesmo da matriz de entrada.fonte
Eu tentei localmente para obter resultados semelhantes, mas estou em uma plataforma de 64 bits para não atingir o limite de memória. Talvez, em vez disso, tente interpolar pequenos pedaços da matriz por vez, como neste exemplo .
Você também pode fazer isso com o GDAL, a partir da linha de comando:
Para fazer a operação equivalente em Python, use ReprojectImage () :
fonte
ReprojectImage
técnica incorporada da GDAL .Eu tive o problema exato no passado e nunca o resolvi usando o interpolate.interp2d. Eu tive sucesso usando scipy.ndimage.map_coordinates . Tente o seguinte:
scipy.ndimage.map_coordinates (matriz_de_ banda, [ax, ay]], ordem = 1)
Isso parece dar a mesma saída que o bilinear.
fonte
scipy.interpolate.interp2d () funciona bem com o scipy mais moderno. Acho que versões mais antigas assumem grades irregulares e não tiram vantagem das grades regulares. Eu recebo o mesmo erro que você faz com o scipy. versão = 0.11.0, mas no scipy. versão = 0.14.0, funciona felizmente em algumas saídas do modelo de 1600x1600.
Obrigado pelas dicas em sua pergunta.
fonte