Interpolação bilinear de dados pontuais em uma varredura em Python?

12

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_errorou fill_valuenã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)

Mike T
fonte
1
Talvez você entenda MemoryErrorporque o NumPy tenta acessar além do seu band_array? Você deve verificar axe ay.
OLT
1
eixo, pode ter um problema se a grade for girada. Talvez seja melhor transformar seus pontos de interpolação em coordenadas de pixel ou dados. Além disso, se houver um problema pontual com eles, você poderá estar indo além do tamanho da banda.
Dave X
As grades rotativas corretas requerem transformação no espaço da grade e depois no espaço das coordenadas. Isso requer o inverso dos coeficientes de transformação afins em gt.
Mike T

Respostas:

7

Traduzi a fórmula abaixo (da Wikipedia ) para Python-speak para produzir o seguinte algoritmo, que parece funcionar.

from numpy import floor, NAN

def bilinear(px, py, no_data=NAN):
    '''Bilinear interpolated point at (px, py) on band_array
    example: bilinear(2790501.920, 6338905.159)'''
    ny, nx = band_array.shape
    # Half raster cell widths
    hx = gt[1]/2.0
    hy = gt[5]/2.0
    # Calculate raster lower bound indices from point
    fx = (px - (gt[0] + hx))/gt[1]
    fy = (py - (gt[3] + hy))/gt[5]
    ix1 = int(floor(fx))
    iy1 = int(floor(fy))
    # Special case where point is on upper bounds
    if fx == float(nx - 1):
        ix1 -= 1
    if fy == float(ny - 1):
        iy1 -= 1
    # Upper bound indices on raster
    ix2 = ix1 + 1
    iy2 = iy1 + 1
    # Test array bounds to ensure point is within raster midpoints
    if (ix1 < 0) or (iy1 < 0) or (ix2 > nx - 1) or (iy2 > ny - 1):
        return no_data
    # Calculate differences from point to bounding raster midpoints
    dx1 = px - (gt[0] + ix1*gt[1] + hx)
    dy1 = py - (gt[3] + iy1*gt[5] + hy)
    dx2 = (gt[0] + ix2*gt[1] + hx) - px
    dy2 = (gt[3] + iy2*gt[5] + hy) - py
    # Use the differences to weigh the four raster values
    div = gt[1]*gt[5]
    return (band_array[iy1,ix1]*dx2*dy2/div +
            band_array[iy1,ix2]*dx1*dy2/div +
            band_array[iy2,ix1]*dx2*dy1/div +
            band_array[iy2,ix2]*dx1*dy1/div)

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.

fórmula de interpolação bilinear

Mike T
fonte
3

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:

gdalwarp -ts $XSIZE*2 0 -r bilinear input.tif interp.tif

Para fazer a operação equivalente em Python, use ReprojectImage () :

mem_drv = gdal.GetDriverByName('MEM')
dest = mem_drv.Create('', nx, ny, 1)

resample_by = 2
dt = (gt[0], gt[1] * resample_by, gt[2], gt[3], gt[4], gt[5] * resample_by)
dest.setGeoTransform(dt)

resampling_method = gdal.GRA_Bilinear    
res = gdal.ReprojectImage(source, dest, None, None, resampling_method)

# then, write the result to a file of your choice...    
scw
fonte
Meus dados pontuais que eu gostaria de interpolar não são regularmente espaçados, portanto não posso usar a ReprojectImagetécnica incorporada da GDAL .
Mike T
1

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.

Matthew Snape
fonte
Fiquei um pouco impressionado com isso, pois não tenho certeza de como as coordenadas de varredura de origem são usadas (em vez de usar coordenadas de pixel). Eu vejo que é "vetorizado" para resolver muitos pontos.
Mike T
Concordo, eu realmente não entendo scipy. Sua solução numpy muito melhor.
Matthew Snape
0

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.

#!/usr/bin/env python

from osgeo import gdal
from numpy import array
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("filename",help='raster file from which to interpolate a (1/3,1/3) point from from')
args = parser.parse_args()

# Read raster
source = gdal.Open(args.filename)
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)])

from scipy import interpolate
bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')

x1 = gt[0] + gt[1]*nx/3
y1 = gt[3] + gt[5]*ny/3.

print(nx, ny, x1,y1,bilinterp(x1,y1))

####################################

$ time ./interp2dTesting.py test.tif 
(1600, 1600, -76.322, 30.70889, array([-8609.27777778]))

real    0m4.086s
user    0m0.590s
sys 0m0.252s
Dave X
fonte