我正在尝试在某些成像数据上进行2D表面。我附上了此类数据的示例，该数据基本上是一个具有大量噪声的1014 x 1014阵列。example_image。该数组的一些补丁是无效的数据，我将其掩盖并将其设置为NAN值，如示例图像中的黄色所示。如您在图像中所看到的那样，我试图删除的背景梯度从左侧（更明亮）到右（调光器）。多项式不能很好地拟合梯度，因此我的目标是进行2D表面双变量样条拟合，然后减去梯度。
我尝试了Scipy中的许多任务，但其中大多数没有返回理想的结果。
-
首先，我尝试了[直接性比例]双变量结构化插值的大型阵列或NAN值或蒙版），但是由于我的图像中有NAN，因此运行rectbivariatespline仅给出了NAN的输出。
-
我还尝试了SmoothBivariateSpline，这是任务的不规则网格版本。我省略了那些具有NAN值并将其余的像素转换为1D数组作为输入。但是由于阵列大小太大而失败了。然后，我试图切碎数组以尝试在较小的块上运行它，但是它会出现以下错误并以分割故障退出，我不知道这是什么意思。
fitpack2.py:1044：用户沃宁：进入错误，未返回近似值。以下条件必须保持：Xb＆lt; = x [i]＆lt; = Xe，Yb＆lt; = y [i]＆lt; = ye，w [i]> 0，i = 0..m-1如果iopt == -1，那就XB
-
i然后尝试首先使用带有线性插值的griddata使用值填充图像中的NAN斑块。由于补丁很大，因此插值不是理想的，但是至少它给了我一个没有NAN的数组。然后，我使用此数组再次运行Rectbivariatespline。但是输出数组仍然是nans。
-
我怀疑我的图像中的噪音正在搞砸这两个任务的行为，因此我还试图在图像上首先运行高斯内核，以使其平滑，然后用Griddata填充Nan Patch，然后运行直接偏移量或平滑型载裂，但它们仍然给我以NAN值作为输出的数组。
我不确定我是否正确理解这两个任务的手册，因此我附上以下脚本：
#!/usr/bin/python
import matplotlib
matplotlib.use('qt5agg')
#matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
#matplotlib.rc('text.latex', preamble=r'usepackage{cmbright}')
#matplotlib.rc('text.latex', preamble=r'usepackage[scaled]{helvet} renewcommandfamilydefault{sfdefault} usepackage[T1]{fontenc}')
#matplotlib.rc('text', usetex=True)
import matplotlib.pyplot as plt
import numpy as np
import astropy.io.fits as pyfits
import scipy.interpolate as sp
from astropy.convolution import convolve
from astropy.convolution import Gaussian2DKernel
#------------------------------------------------------------
#Read in the arrays
hdulistorg = pyfits.open('icmj01jrq_flt.fits')
hdulistorg.info()
errarrorg = np.swapaxes(hdulistorg[1].data, 0,1)
hdulist = pyfits.open('jrq_sci_nan_deep.fits')
hdulist.info()
dataarrorg = np.swapaxes(hdulist[0].data, 0,1) #image array
errarrorg = np.swapaxes(hdulistorg[1].data, 0,1) #error array
#Flag some of the problematic values, turn NaNs into 0 for easier handling
dataarr = np.copy(dataarrorg)
w=np.isnan(dataarr)
ww=np.where(dataarr == 0)
www=np.where(dataarr > 100)
wwww=np.where(dataarr < 0)
errarr = 1.0 / (np.copy(errarrorg)+1e-5) # Try to use 1/error as the estimate for weight below
errarr[w] = 0
errarr[ww] = 0
errarr[www] = 0
errarr[wwww]=0
dataarr[w]= 0
dataarr[ww]= 0
dataarr[www]=0
dataarr[wwww]=0
#Make a gaussian kernel smoothed data
maskarr = np.copy(errarr) #For masking the nan regions so they dun get smoothed
maskarr[:]=0
maskarr[w]=1
maskarr[ww]=1
maskarr[www]=1
maskarr[wwww]=1
gauss = Gaussian2DKernel(stddev=5)
condataarr = convolve(dataarr,gauss,normalize_kernel=True,boundary='extend',mask=maskarr)
condataarr[w]=0
conerrarr = np.copy(errarr)
#Setting x,y arrays for the Spline functions
nx, ny = (1014,1014)
x = np.linspace(0, 1013, nx)
y = np.linspace(0, 1013, ny)
xv, yv = np.meshgrid(x, y)
#Make an 1D version of these 2D arrays
dataarrflat = np.ravel(condataarr[0:200,0:200]) #Try only a small chunk!
xvflat = np.ravel(xv[0:200,0:200])
yvflat = np.ravel(yv[0:200,0:200])
errarrflat = np.ravel(conerrarr[0:200,0:200])
notnanloc = np.where(dataarrflat != 0) #Not NaNs
#SmoothBivariateSpline!
rect_S_spline = sp.SmoothBivariateSpline(xvflat[notnanloc], yvflat[notnanloc], dataarrflat[notnanloc],w=errarrflat[notnanloc], kx=3, ky=3)
#Also try using grid data to fix the grid?
gddataarr = np.copy(condataarr)
gddataarrflat = np.ravel(gddataarr)
gdloc = np.where(gddataarrflat != 0) #Not NaNs
gdxvflat = np.ravel(xv)
gdyvflat = np.ravel(yv)
xyarr = np.c_[gdxvflat[gdloc],gdyvflat[gdloc]]
x_grid, y_grid = np.mgrid[0:1013:1014j,0:1013:1014j]
grid_z2 = sp.griddata(xyarr, gddataarrflat[gdloc], (x_grid, y_grid), method='linear')
plt.imshow(grid_z2.T)
#plt.show()
#RectBivariatSpline
rect_B_spline = sp.RectBivariateSpline(x, y, grid_z2.T)
#Result grid (same as input for now)
xnew = np.arange(0, 1013, 1)
ynew = np.arange(0, 1013, 1)
znewS = rect_S_spline(xnew, ynew)
znewB = rect_B_spline(xnew, ynew)
print 'znewS', znewS
print 'znewB', znewB
#Write FITS files
condataarr = np.swapaxes(condataarr, 0, 1)
hdu2 = pyfits.PrimaryHDU(condataarr)
hdulist2 = pyfits.HDUList([hdu2])
hdulist2.writeto('contest.fits',overwrite=True)
hdulist2.close()
hdu3 = pyfits.PrimaryHDU(znewS)
hdulist3 = pyfits.HDUList([hdu3])
hdulist3.writeto('Stest.fits',overwrite=True)
hdulist3.close()
我无法完全解决您的问题，但是我有一些代码与Python插入了Fortran插值例程。您只需直接从Python调用例程，就不需要Fortran。
您可以在此GitHub页面上找到代码和描述https://github.com/haakoan/inter