python - ndimage map_coordinates with masked arrays

Question

I am using ndimage to interpolate as follows:

ndimage.map_coordinates(input_data, coords, output, 2, prefilter=False)

Now, the problem is that I do not have valid measurements over my whole input data. So, I have a masked array which tells me which data points are valid. So, when doing the interpolation, I would like to only use pixels which are valid and also adjust the weightings accordingly (to ensure that the weights sum to 1).

However, I see that there is no easy way to do this. I was wondering if someone knows of a good way to do this or can point me to some library or code that I can use. I am coming from a C++ background, so still finding my way around python.

Answer 1

It sounds like you need to focus on the interpolation of data and then extract values from desired coordinates. For 1D splrep and 2D bisplrep are the interpolation functions you need to check out (A good overview). Both of these functions can be weighted and provide fine tune control over the spline function you interpolate with.

Once you have filtered the data with the desired weights you can then determine the value at specified coordinates using.

ndimage.map_coordinates(input_data, coords, output, prefilter=True)

note the prefilter key word argument is not needed as that is the default value

Answer 2

You can use scipy.interpolate.griddata to get interpolated values from arbitrary known data points. This example

import numpy as np
from scipy.interpolate import griddata

# data array
a = np.arange(9).reshape((3,3)).astype(float)
a[1, 1] = np.nan
print(a)
a = a.flatten()

# coordinate arrays
ii, jj = np.indices((3,3))
ij = np.stack((ii.flatten(), jj.flatten()), axis=1)

# filter out unknowns
mask = ~np.isnan(a)
a = a[mask]
ij = ij[mask]

# interpolate for the missing a[1, 1] element
res = griddata(ij, a, (1, 1), method='cubic')
print(res)

produces

[[ 0.  1.  2.]
 [ 3. nan  5.]
 [ 6.  7.  8.]]

4.000000157826586