差异 scipy interpolate vs mpl griddata

Posted 2023-03-27

【中文标题】差异 scipy interpolate vs mpl griddata

【英文标题】：differences scipy interpolate vs mpl griddata

【发布时间】：2021-12-31 22:15:36

【问题描述】：

我差点把头发扯下来，真的希望有人能在这里帮助我。

我需要将一个脚本从 python 2 转换为 python 3，它是一个简单的脚本，它接受坐标点x、y，以及这些点的值z。

z 在网格上插值，然后绘制。 python 2脚本中的插值目前正在使用from matplotlib.mlab import griddata

但是，python 3 中不存在 mlab griddata。所以我认为我们应该使用 scipy.interpolate 作为替代品。

问题是，它给了我非常奇怪的结果。

最近邻比较

Mpl griddata With interp='nn' 给我：

scipy interpolate using method='nearest' 给了我

这是最重要的，因为这是我需要的结果，但结果却如此不同！

我尝试了其他参数看看是否仍然存在

线性比较

带有interp='linear'的mpl griddata给了我

scipy interpolate with method='linear' 给了我

线性显示相同的结果，但我不在乎线性。

三次比较

mpl griddata 没有三次，但是使用 scipy interpolate 的结果非常奇怪。

数据

import numpy as np
# points
x = np.array([31132.71118116, 10763.66383076, 29340.62119857,  4025.2544491 ,
       19158.87683925, 36584.96725821, 17335.28090098, 27279.80878719,
        5893.32833709, 36392.22873482, 22350.16973122, 40235.86525991,
       19844.57893466, 11245.50791534, 17103.62528354, 29944.77964486,
       34021.92540595, 44485.36307614,  8797.36817447, 24733.44616607,
       30628.62812172, 27227.91858716, 26221.60497488, 27987.74985745,
       14366.37775758, 32035.24328142, 29779.84095199, 42821.14608703,
       42229.95963477, 22592.94483609, 26362.45107434, 22161.13844494,
       21562.53974235, 40784.42626907, 14195.87932099, 25628.56490406,
       19212.62065784, 28354.69294259, 34299.91259166, 28276.85653059,
       13291.49360218, 64014.88912569,     0.        , 64014.88912569,
           0.        ])
    
y = np.array([1.61640041e+04, 1.18226221e+04, 1.66496799e+04, 1.60268291e+04,
       1.47533260e+04, 1.72743862e+04, 2.58461610e+04, 1.11997818e+04,
       1.94747705e+04, 2.81559903e+04, 2.85786579e+04, 2.16418047e+04,
       2.58684064e+04, 2.17474690e+04, 3.21398453e+04, 2.13563269e+04,
       2.12339790e+04, 2.42537633e+04, 1.89242092e+04, 2.45410984e+04,
       1.90984612e+04, 2.45058762e+04, 1.76859142e+04, 2.13655953e+04,
       2.62595562e+04, 1.44826845e+04, 2.55328718e+04, 2.97076285e+04,
       1.81975446e+04, 3.27182416e+04, 1.88407910e+04, 1.58599943e+04,
       2.02514897e+04, 2.42723003e+04, 1.76265941e+04, 1.60045845e+04,
       2.03274941e+04, 1.67609034e+04, 1.74986867e+04, 1.78898244e+04,
       3.21546756e+04, 9.55042196e-08, 9.55042196e-08, 4.04414440e+04,
       4.04414440e+04])
    
z = np.array([0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.6, 0. , 0. , 0. , 0. ,
       0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
       0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
       0. , 0. , 0. , 0. , 0. , 0. ])
    
min_grid_val_x = 0
max_grid_val_x = 64014.88912569069

min_grid_val_y = 9.550421964377165e-08
max_grid_val_y = 40441.443993242734
    
# set up a square grid with the same extents as our measured data
numcols, numrows = 1000, 1000
xi = np.linspace(min_grid_val_x, max_grid_val_x, numcols)
yi = np.linspace(min_grid_val_y, max_grid_val_y, numrows)
xi, yi = np.meshgrid(xi, yi)

其他代码

import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.colors as mc
import matplotlib.patches as mpatches
from matplotlib.colors import ListedColormap
from matplotlib.colors import Normalize
import math

MODE == 'mpl'
if MODE == 'mpl':
    from matplotlib.mlab import griddata
    zi = griddata(x, y, z, xi, yi, interp='nn')
else:
    from scipy.interpolate import griddata
    points = np.array([[x[i],y[i]] for i in range(len(x))])
    zi = griddata(points, z, (xi, yi), method='cubic',fill_value=0)
    
min_val = math.floor(z.min())
max_val = math.ceil(z.max())
resolutions = [0.1, 0.2, 0.5, 1, 1.5, 2, 3, 5, 10, 
               15, 20, 30, 40, 50, 60, 70, 80]
    
max_of_Z = z.max()
if (max_of_Z==0):
    levels = np.arange(min, min + (0.1 * 12), 0.1)
else:
    hit = False
    for res in resolutions:
        condition = res*12
        if condition > z.max():
            print("Using ".format(res))
            levels = np.arange(min_val, max_val + condition, res)
            hit = True
            break
    if not hit:
        res = 90
        print("Using ".format(res))
        condition = res*12
        levels = np.arange(min_val, max_val + condition, res)
    
norm = mc.BoundaryNorm(levels, len(levels)-1)

cmapV=ListedColormap([[0.73,0.90,0.97],
[0.66,0.87,0.97],
[0.57,0.82,0.95],
[0.49,0.76,0.91],
[0.40,0.69,0.87],
[0.32,0.62,0.83],
[0.26,0.53,0.77],
[0.20,0.44,0.71],
[0.14,0.34,0.64],
[0.13,0.29,0.56],
[0.12,0.23,0.47],
[0.11,0.17,0.38]])

# define map extent
lllat=1.1497
urlat=1.5133
lllon=103.5822
urlon=104.1579
# set up plot
plt.clf()
fig = plt.figure(figsize=(13, 6))
ax = fig.add_subplot(111, facecolor='w', frame_on=False)

# Set up Basemap instance
m = Basemap(
    projection = 'merc',
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat,
    resolution='h')

# contour plot
con = m.contourf(xi, yi, zi, zorder=1, alpha=None, cmap=cmapV, levels=levels, norm=norm)

【参考方案1】：

根据STOQS 的问题Github 页面上的this discussion，最接近旧Matplotlib 的griddata 默认算法，即自然邻居，是使用scipy.interpolate.griddata 获得的与method='cubic' 和rescale=True。

也就是说，如果您发布的数据代表了您的实际问题，那么您的图表都没有洞察力，因为您不应该从单个数据点 (x=5893.32833709, y=19474.7705, z=0.6) 进行推断，而所有其他数据点都有一个 z纵坐标全等于零。

【讨论】：

嘿，感谢您的回复！我以前看过那个帖子，它对我也不起作用。尝试griddata(points, z, (xi, yi), method='cubic',fill_value=0, rescale=True) 给了我与问题中最后一个情节几乎完全相同的情节。如果它只给出左上角的区域（即大多数点为 0），我会接受这个图。是的，这是真实数据，将是一种常见情况。这种情况下的预期结果是理想情况下的第一个图 (NN)

@Wboy ①如果你在Q中提到你已经看过那个问题，我可能不会浪费时间回答②当你的数据有结构时，插值结果是明智的，尽管有算法;相反……

抱歉浪费您的时间。我正在寻找一种在 python 3 中复制 mpl 的旧 griddata 自然邻居插值的方法。