# Bilinear interpolation of point data on a raster in Python?

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

source = gdal.Open('my_raster.tif')
nx, ny = source.RasterXSize, source.RasterYSize
gt = source.GetGeoTransform()
# 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)])
``````

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 `bounds_error` or `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[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]
``````

... but I much prefer bilinear interpolation methods)

-
Maybe you get the `MemoryError` because NumPy tries to access beyond your `band_array`? You should check `ax` and `ay`. – olt Mar 23 '11 at 9:36
ax,ay might have a problem if the grid is rotated at all. It might be better to transform your interpolate-to points into pixel or data coordinates. Also, if there is an off-by-one problem with them, then you might be going beyond the size of the band. – Dave X Aug 4 '14 at 18:29
Correct, rotated grids require transformation to grid-space then back to coordinate-space. This requires the inverse of the affine transform coefficients in `gt`. – Mike T Aug 4 '14 at 21:07

I've translated the formula below (from Wikipedia) into Python-speak to yield the following algorithm, which appears to work.

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

Note that the result will be returned with an apparent higher precision than the source data, since it is up-classed to NumPy's `dtype('float64')` data type. You can use the return value with `.astype(band_array.dtype)` to make the output data type the same as the input array.

-
I've been comparing this with ArcGIS's Surface Spot results, and all I see are identical results. – Mike T Mar 24 '11 at 2:34

I tried it locally to similar results, but I'm on a 64-bit platform so it didn't hit the memory limit. Perhaps instead try interpolating small pieces of the array at a time, like in this example.

You could also just do this with GDAL, from the command line it'd be:

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

To do the equivalent operation in 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...
``````
-
My point data that I'd like to interpolate are not regularly spaced, so I can't used GDAL's built-in `ReprojectImage` technique. – Mike T Mar 23 '11 at 20:55

I have had the exact issue in the past, and never resolved it using interpolate.interp2d. I have had success using scipy.ndimage.map_coordinates. Try the following:

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

This seems to give the same output as bilinear.

-
I was a bit thrown off by this one, as I'm unsure how the source raster coordinates are used to be used (rather than using pixel coordinates). I see it is "vectorized" to solve many points. – Mike T Mar 24 '11 at 1:29
Agreed, I don't really understand scipy. Your numpy solution much better. – Matthew Snape Mar 24 '11 at 11:24

scipy.interpolate.interp2d() works fine with more modern scipy. I think older versions assume irregular grids and don't take advantage of the regular grids. I get the same error as you do with scipy.version = 0.11.0, but on scipy.version = 0.14.0, it happily works on some 1600x1600 model output.

Thank you for the hints in your 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()

source = gdal.Open(args.filename)
nx, ny = source.RasterXSize, source.RasterYSize
gt = source.GetGeoTransform()
# 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
``````
-