Interpolation bilinéaire de données ponctuelles sur un raster en Python?

J'ai un raster avec lequel j'aimerais faire quelques interpolations ponctuelles. Voici où j'en suis:

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)])

Jusqu'à présent, j'ai essayé la fonction interp2d de SciPy :

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

cependant, j'obtiens une erreur de mémoire sur mon système Windows 32 bits avec un raster 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

Je l'admets, j'ai une confiance limitée dans cette fonction SciPy, car les paramètres bounds_errorou fill_valuene fonctionnent pas comme indiqué. Je ne vois pas pourquoi je devrais avoir une erreur de mémoire, car mon raster est 317 × 301, et l' algorithme bilinéaire ne devrait pas être difficile.

Quelqu'un a-t-il rencontré un bon algorithme d'interpolation bilinéaire, de préférence en Python, éventuellement adapté avec NumPy? Des astuces ou des conseils?

(Remarque: l' algorithme d'interpolation du plus proche voisin est un gâteau facile:

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]

... mais je préfère de loin les méthodes d'interpolation bilinéaire)

J'ai traduit la formule ci-dessous (de Wikipedia ) en langage Python pour produire l'algorithme suivant, qui semble fonctionner.

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)

Notez que le résultat sera retourné avec une précision apparemment plus élevée que les données source, car il est classé dans le dtype('float64')type de données NumPy . Vous pouvez utiliser la valeur de retour avec .astype(band_array.dtype)pour que le type de données de sortie soit identique au tableau d'entrée.

Je l'ai essayé localement pour des résultats similaires, mais je suis sur une plate-forme 64 bits, donc il n'a pas atteint la limite de mémoire. Essayez plutôt d'interpoler de petits morceaux du tableau à la fois, comme dans cet exemple .

Vous pouvez également le faire avec GDAL, à partir de la ligne de commande, ce serait:

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

Pour effectuer l'opération équivalente en Python, utilisez 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...    

J'ai eu le problème exact dans le passé et je ne l'ai jamais résolu en utilisant interpolate.interp2d. J'ai réussi à utiliser scipy.ndimage.map_coordinates . Essayez ce qui suit:

scipy.ndimage.map_coordinates (band_array, [ax, ay]], order = 1)

Cela semble donner la même sortie que bilinéaire.

scipy.interpolate.interp2d () fonctionne très bien avec scipy plus moderne. Je pense que les anciennes versions supposent des grilles irrégulières et ne profitent pas des grilles régulières. J'obtiens la même erreur que vous avec scipy. version = 0.11.0, mais sur scipy. version = 0.14.0, il fonctionne heureusement sur une sortie de modèle 1600x1600.

Merci pour les indices dans votre question.

#!/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