NaN-friendly 2D convolution using Weave. Also very useful for convolving with big kernels, where scipy's convolve routines run out of RAM.

nan_convolve.py

from scipy.weave import converters
import scipy.weave as weave
import numpy as np


def convolve(image, kernel):

    image = image.astype(float)
    kernel = kernel.astype(float)

    nx, ny = image.shape
    nkx, nky = kernel.shape

    if nkx % 2 == 0 or nky % 2 == 0:
        raise Exception("Kernel dimensions should be odd")

    smoothed = np.zeros(image.shape)
    isvalid = ~np.isnan(image)

    code = """
            double top, bot;
            int wkx, wky, iimin, iimax, jjmin, jjmax;
            wkx = (nkx-1)/2;
            wky = (nky-1)/2;
            for (int i=0; i<nx; ++i) {
                for (int j=0; j<ny; ++j) {
                    if(isvalid(i,j)) {
                        top = 0.;
                        bot = 0.;
                        if(i-wkx > 0) { iimin = i-wkx; } else { iimin = 0; };
                        if(i+wkx < nx-1) { iimax = i+wkx; } else { iimax = nx-1; };
                        if(j-wkx > 0) { jjmin = j-wky; } else { jjmin = 0; };
                        if(j+wkx < ny-1) { jjmax = j+wky; } else { jjmax = ny-1; };
                        for (int ii=iimin; ii <= iimax ; ++ii) {
                            for (int jj=jjmin; jj <= jjmax; ++jj) {
                                if(isvalid(ii,jj)) {
                                    top = top + kernel(wkx + ii-i,wky + jj-j) * image(ii,jj);
                                    bot = bot + kernel(wkx + ii-i,wky + jj-j);
                                }
                            }
                        }
                        smoothed(i,j) = top / bot;
                    } else {
                        smoothed(i,j) = image(i,j);
                    }
                }
            }
            return_val = 1;
            """

    weave.inline(code, ['image', 'isvalid', 'nx', 'ny', 'kernel', 'nkx', 'nky', 'smoothed'],
                 type_converters=converters.blitz, compiler = 'gcc')

    return smoothed