Bilineare Interpolation von Punktdaten auf einem Raster in Python?

12

Ich habe ein Raster, mit dem ich einige Punktinterpolationen durchführen möchte. Hier bin ich:

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

Bisher habe ich die interp2d- Funktion von SciPy ausprobiert :

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

Auf meinem 32-Bit-Windows-System mit einem Raster von 317 × 301 wird jedoch ein Speicherfehler angezeigt:

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

Ich gebe zu , ich begrenzte Vertrauen in diese SciPy Funktion haben, wie die bounds_erroroder fill_valueParameter funktionieren nicht dokumentiert. Ich verstehe nicht, warum ich einen Speicherfehler haben sollte, da mein Raster 317 × 301 ist und der bilineare Algorithmus nicht schwierig sein sollte.

Hat jemand einen guten bilinearen Interpolationsalgorithmus gefunden, vorzugsweise in Python, der möglicherweise auf NumPy zugeschnitten ist? Irgendwelche Hinweise oder Ratschläge?


(Hinweis: Der Interpolationsalgorithmus für den nächsten Nachbarn ist einfach:

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]

... aber ich bevorzuge bilineare Interpolationsmethoden)

Mike T
quelle
1
Vielleicht bekommst du das, MemoryErrorweil NumPy versucht, über deine Grenzen hinaus zuzugreifen band_array? Sie sollten axund überprüfen ay.
Alt
1
Axt, ay könnte ein Problem haben, wenn das Gitter überhaupt gedreht wird. Es ist möglicherweise besser, die Interpolationspunkte in Pixel- oder Datenkoordinaten umzuwandeln. Auch wenn es ein Problem gibt, bei dem es nur um eins geht, überschreiten Sie möglicherweise die Größe der Band.
Dave X
Richtige, gedrehte Gitter erfordern eine Transformation in den Gitterraum und dann zurück in den Koordinatenraum. Dies erfordert die Inverse der affinen Transformationskoeffizienten in gt.
Mike T

Antworten:

7

Ich habe die folgende Formel (aus Wikipedia ) in Python-Sprache übersetzt, um den folgenden Algorithmus zu erhalten, der zu funktionieren scheint.

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)

Beachten Sie, dass das Ergebnis mit einer offensichtlich höheren Genauigkeit als die Quelldaten zurückgegeben wird, da es dem dtype('float64')Datentyp von NumPy zugeordnet ist . Sie können den Rückgabewert mit verwenden .astype(band_array.dtype), um den Ausgabedatentyp mit dem Eingabearray gleichzusetzen.

bilineare Interpolationsformel

Mike T
quelle
3

Ich habe es lokal mit ähnlichen Ergebnissen versucht, aber ich arbeite auf einer 64-Bit-Plattform, sodass das Speicherlimit nicht überschritten wurde. Versuchen Sie stattdessen, kleine Teile des Arrays gleichzeitig zu interpolieren, wie in diesem Beispiel .

Sie könnten dies auch einfach mit GDAL tun, von der Kommandozeile aus wäre es:

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

Verwenden Sie ReprojectImage () , um die entsprechende Operation in Python auszuführen :

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
quelle
Meine Punktdaten, die ich interpolieren möchte, sind nicht regelmäßig verteilt, daher kann ich die integrierte ReprojectImageTechnik von GDAL nicht verwenden .
Mike T
1

Ich hatte das genaue Problem in der Vergangenheit und habe es nie mit interpolate.interp2d gelöst. Ich hatte Erfolg mit scipy.ndimage.map_coordinates . Versuche Folgendes:

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

Dies scheint die gleiche Ausgabe zu liefern wie bilinear.

Matthew Snape
quelle
Ich war ein bisschen verärgert, da ich nicht sicher bin, wie die Quell-Rasterkoordinaten verwendet werden (anstatt Pixelkoordinaten zu verwenden). Ich sehe, es ist "vektorisiert", um viele Punkte zu lösen.
Mike T
Einverstanden, ich verstehe Scipy nicht wirklich. Ihre numpy Lösung viel besser.
Matthew Snape
0

scipy.interpolate.interp2d () funktioniert gut mit moderneren scipy. Ich denke, ältere Versionen gehen von unregelmäßigen Rastern aus und nutzen die regulären Raster nicht aus. Ich bekomme den gleichen Fehler wie bei scipy. version = 0.11.0, aber auf scipy. version = 0.14.0, funktioniert problemlos mit einigen 1600x1600-Modellausgaben.

Vielen Dank für die Hinweise in Ihrer Frage.

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