Interpolación scipy con datos enmascarados?

Estoy intentando interpolar una matriz 2D que contiene datos enmascarados. He utilizado algunos de los métodos disponibles del módulo SciPy, incluidos interp2d , bisplrep/bisplev , así como RectBivariateSpline . Como información adicional, mis datos son una matriz regular, lo que significa que las cuadrículas tienen la misma dimensión (en este caso, 1ºX1º).

Dicho esto, ¿hay alguna forma de interpolar evitando datos enmascarados en una matriz con Python? Todavía soy nuevo usando los módulos Python y NumPy / SciPy.

Puede usar cada función que acepte x, y, z (que es el caso de interp2d y probablemente también las otras) con sus datos enmascarados. Pero necesitas crear explícitamente un mgrid :

 z = ... # Your data x, y = np.mgrid[0:z.shape[0], 0:z.shape[1]] 

Entonces necesitas eliminar todos los valores enmascarados en todas estas coordenadas:

 x = x[~z.mask] y = y[~z.mask] z = z[~z.mask] 

Con estos x, y, z finales x, y, z puede llamar a cada una de sus funciones especificadas (que acepta cuadrículas incompletas, por lo que RectBivariateSpline no funcionará). Sin embargo, tenga en cuenta que algunos de estos usan cuadros de interpolación, por lo que si hay una región demasiado grande donde descartó los datos debido a su máscara, la interpolación fallará allí (lo que resultará en np.nan o 0). Pero podría ajustar los parámetros para compensar eso, si sucede.

Por ejemplo:

 data = np.random.randint(0, 10, (5,5)) mask = np.random.uniform(0,1,(5,5)) > 0.5 z = np.ma.array(data, mask=mask) x, y = np.mgrid[0:z.shape[0], 0:z.shape[1]] x1 = x[~z.mask] y1 = y[~z.mask] z1 = z[~z.mask] interp2d(x1, y1, z1)(np.arange(z.shape[0]), np.arange(z.shape[1])) array([[ 1.1356716 , 2.45313727, 3.77060294, 6.09790177, 9.31328935], [ 3.91917937, 4. , 4.08082063, 3.98508121, 3.73406764], [ 42.1933738 , 25.0966869 , 8. , 0. , 0. ], [ 1.55118338, 3. , 4.44881662, 4.73544593, 4. ], [ 5. , 8. , 11. , 9.34152525, 3.58619652]]) 

puede ver el área pequeña de 0 porque la máscara tenía muchos valores enmascarados:

 mask array([[False, True, True, True, False], [False, False, True, False, False], [ True, True, False, True, True], [False, True, False, True, True], [False, True, False, False, True]], dtype=bool) data array([[2, 4, 4, 5, 5], [1, 4, 1, 3, 8], [9, 1, 8, 0, 9], [7, 2, 0, 3, 4], [9, 6, 0, 4, 4]]) 

Por lo general, sigo el enfoque descrito por @mseifert, pero agrego el siguiente refinamiento si estoy cansado del error de interpolación a través de las áreas enmascaradas. Esa parece ser una de tus preocupaciones, @hurrdrought? la idea es propagar la máscara al resultado interpolado. Un ejemplo simple para datos 1D es:

 def ma_interp(newx,x,y,mask,propagate_mask=True): newy = np.interp(newx,x[~mask],y[~mask]) # interpolate data if propagate_mask: # interpolate mask & apply to interpolated data newmask = mask[:] newmask[mask] = 1; newmask[~mask] = 0 newmask = np.interp(newx,x,newmask) newy = np.ma.masked_array(newy, newmask>0.5) return newy 

El problema con los enfoques descritos por @MSeifert es que la estructura de cuadrícula regular se pierde, lo que resulta en una interpolación ineficiente. Solo se justifica rellenar los datos faltantes por interpolación, pero no para la interpolación típica de una cuadrícula a otra, donde no se deben completar los datos faltantes.

En este caso, llenar los valores perdidos con np.nan es el enfoque más simple. Estos se propagarán en los cálculos y la matriz resultante tendrá nans siempre que se haya utilizado un valor faltante para la interpolación.

 # fast interpolator that use the regular grid structure (x and y are 1D arrays) z = z_masked.filled(np.nan) zinterp = RegularGridInterpolator((x, y), zT) # new grid to interpolate on X2, Y2 = np.meshgrid(x2, y2) newpoints = np.array((X2, Y2)).T # actual interpolation z2 = zinterp(newpoints) z2_masked = np.ma.array(z2, mask=np.isnan(z2)) 

Para completar, otro enfoque es interpolar una segunda matriz de mask (rellena con 1 donde faltan datos) para completar los valores faltantes en la nueva cuadrícula.

 # fast interpolator that use the regular grid structure (x and y are 1D arrays) zinterp = RegularGridInterpolator((x, y), zT) minterp = RegularGridInterpolator((x, y), (mask+0.).T) # actual interpolation z2 = zinterp(newpoints) mask2 = minterp(newpoints) > 0 # apply threshold, eg 0.5 is considered contaminated and will be removed. z2[mask2] = np.nan # fill with nans or whatever missing data flag 

Tenga en cuenta que ambos enfoques también deberían funcionar con RectBivariateSpline, si se desea la interoplación spline. Y de cualquier manera, esto debería ser mucho más rápido que usar interp2d …