I have a raster that I'd like to do some point interpolations with. Here is where I'm at:
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 + ix*gt + gt/2.0 for ix in range(nx)]) ay = array([gt + iy*gt + gt/2.0 for iy in range(ny)])
Up to now, I've tried SciPy's interp2d function:
from scipy import interpolate bilinterp = interpolate.interp2d(ax, ay, band_array, kind='linear')
however I get a memory error on my 32-bit Windows system with a 317×301 raster:
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
I'll admit, I have limited confidence in this SciPy function, as the
fill_value parameters don't work as documented. I don't see why I should have a memory error, since my raster is 317×301, and the bilinear algorithm should not be difficult.
Has anyone come across a good bilinear interpolation algorithm, preferably in Python, possibly tailored with NumPy? Any hints or advice?
(Note: the nearest neighbor interpolation algorithm is easy cake:
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 + gt/2.0))/gt)) iy = int(round((py - (gt + gt/2.0))/gt)) if (ix < 0) or (iy < 0) or (ix > nx - 1) or (iy > ny - 1): return no_data else: return band_array[iy, ix]
... but I much prefer bilinear interpolation methods)