二维散布点的等值线

随机的散布点

下面我们首先生成二维平面上的随机散列点(x,y)，以及每个点对应的数值v。

mport numpy as np
np.random.seed(0)
x = np.random.normal(size=200)
y = np.random.normal(size=200)
v = np.sqrt(x**2+y**2)

我们可以用matplotlib的scatter()绘制散布点：

import matplotlib.pyplot as plt
plt.scatter(x, y, c=v)
plt.color()
plt.show()

griddata的用法

griddata()的调用方式如下：

griddata(points, values, xi, method='linear', fill_value=nan)

其中points为二维平面上坐标点数据，可以是形状为(N, 2)的数组，其中N为坐标点数；或者为一个包含各个轴数据的元组。values为每个坐标点所对应的数值，它是一个长度为N的一维数组。xi为需要进行插值的坐标点，我们将使用包含各个轴坐标的元组。method为插值方式，有’linear’、’nearest’和’cubic’三种方式。

griddata()也支持多维空间散列点的插值。

下面让我们通过一个例子说明各个参数的用法。

import numpy as np
from scipy.interpolate import griddata
from enthought.mayavi import mlab

x = np.array([-1, 0, 0, 1])
y = np.array([0, -2, 2, 0])
v = np.array([0, 1.5, 1.5, 0])

xg, yg = np.mgrid[-2:2:100j, -2:2:100j]
for method in ['nearest', 'linear', 'cubic']:
    vg = griddata((x, y), v, (xg, yg), method=method)
    mlab.figure()
    mlab.surf(xg, yg, vg)
    mlab.points3d(x, y, v, v, scale_mode="none", scale_factor=0.1)
    mlab.axes()
mlab.show()

在这个例子中，我们对二维平面上的(-1, 0), (0, -2), (0, 2), (1, 0)等4点的数据进行插值，这些点所对应的数值分别为0, 1.5, 1.5, 0。为了直观地观察插值结果，我们使用Mayavi将结果用三维空间中的曲面显示，曲面的高度就为平面上每点所对应的插值。

其运行结果如下图所示，图中使用蓝色和红色的圆点表示原始数据点，曲面为使用原始数据点进行插值的结果：

由图可知，’nearest’找到离插值点最近的原始数据点，并使用它所对应的数值作为插值点的值。’nearest’能够计算平面上的任意点的插值。而’linear’和’cubic’使用线性和三次方插值，它们只能对原始数据点的凸包所包围的区域进行插值。由于’linear’为线性插值，因此整个插值曲面是由许多小三角形面构成。

将散布点用栅格点插值

现在我们将前面的散布点用griddata进行栅格化：

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
np.random.seed(0)

x = np.random.normal(size=200)
y = np.random.normal(size=200)
v = np.sqrt(x**2+y**2)
xg, yg = np.mgrid[-2:2:100j, -2:2:100j]
vg = griddata((x, y), v, (xg, yg), method='cubic')
plt.contourf(xg, yg, vg)
plt.show()

我们用matplotlib的contourf()将栅格数据用等值线显示。其效果如下：

如果使用pyplot.imshow()显示栅格数据，需要对vg进行转置。

从栅格数据抽取等值线

matplotlib提供了contour()和contourf()绘制栅格数据的等值线，在《Python科学计算》中介绍过从所绘制的等值线提取坐标信息的方法，使用这种方法的缺点是需要先绘制等值图。实际上我们只需要稍微分析一下matplotlib的源码就可以找到计算等值线的程序，从而实现直接计算等值线。下面先看完整的源程序：

import matplotlib._cntr as _cntr ❶

c = _cntr.Cntr(xg, yg, vg) ❷
plt.figure()
for v in [0, 0.5, 1.0, 1.5]:
    t = c.trace(v) ❸
    for p in t[:len(t)//2]: ❹
        if np.all(np.isnan(p)): ❺
            continue
        plt.plot(p[:,0], p[:,1], lw=2) ❻
        r = np.sqrt(np.sum(p**2, 1)) ❼
        print r
plt.show()

❶在matplotlib中计算等值线由_cntr模块中的Cntr类实现，❷首先用栅格点坐标xg、yg和各点所对应的值vg，创建一个Cntr对象。❸调用其trace(v)方法可以获得平面中与值v对应的坐标点。

❹由于等值线不一定是完全连续的，因此trace()可能会返回多条等值线的信息。它返回的数据分为两部分，前半部分为每条等值线上各点的值，而后半部分为各点的绘制方法，这部分在matplotlib绘制等值线时使用，在我们的程序中可以将其忽略。

❺由于我们的栅格上并非所有的点都对应有效数据，因此trace()可能会返回许多全部为NaN的数据，在对结果进行处理时，需要将其忽略。❻等值线上各点的坐标采用形状为(N, 2)的数组保存。❼为了验证结果，我们用各点坐标计算按照公式直接计算其对应的值。

程序的输出如下：