Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
1  
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

3 Answers 3

up vote 3 down vote accepted

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.

bilinear interpolation formula

share|improve this answer
    
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...    
share|improve this answer
    
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.

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.