5

我在 matplotlib 中生成地下水高程等值线和流图

等高线表明许多地区的海拔正在下降,但地下水流(流图)指向上坡。我圈出了似乎指向错误方向的箭头。

指向地图底部的箭头似乎指向了正确的方向。有谁知道这可能是为什么?

在此处输入图像描述

这是生成此图的大部分代码:

#create empty arrays to fill up!
x_values = []
y_values = []
z_values = []

#iterate over wells and fill the arrays with well data
for well in well_arr:
    x_values.append(well['xpos'])
    y_values.append(well['ypos'])
    z_values.append(well['value'])

#initialize numpy array as required for interpolation functions
x = np.array(x_values, dtype=np.float)
y = np.array(y_values, dtype=np.float)
z = np.array(z_values, dtype=np.float)

#create a list of x, y coordinate tuples
points = zip(x, y)

#create a grid on which to interpolate data
xi, yi = np.linspace(0, image['width'], image['width']),
         np.linspace(0, image['height'], image['height'])
xi, yi = np.meshgrid(xi, yi)

#interpolate the data with the matlab griddata function
zi = griddata(x, y, z, xi, yi, interp='nn')

#create a matplotlib figure and adjust the width and heights
fig = plt.figure(figsize=(image['width']/72, image['height']/72))

#create a single subplot, just takes over the whole figure if only one is specified
ax = fig.add_subplot(111, frameon=False, xticks=[], yticks=[])

#create the contours
kwargs = {}
if groundwater_contours:
    kwargs['colors'] = 'b'

CS = plt.contour(xi, yi, zi, linewidths=linewidth, **kwargs)

#add a streamplot
dx, dy = np.gradient(zi)
plt.streamplot(xi, yi, dx, dy, color='c', density=1, arrowsize=3)
4

2 回答 2

19

概括

我猜,但你的问题可能是因为你有一个固有的转置正在进行。2D numpy 数组被索引为行、列。"x, y" 索引是列、行。在这种情况下,numpy.gradient基本上是要返回 dy, dx 而不是 dx, dy。

尝试换行:

dx, dy = np.gradient(zi)

到:

dy, dx = np.gradient(zi)

此外,如果您的深度被定义为向上,它应该是:

dy, dx = np.gradient(-zi)

但是,我假设您有正向下的深度约定,所以我将在下面的示例中保留这部分内容。(因此在下面的示例数据中假设较高的值更深/更低,水将流向高值。)

重现问题

例如,如果我们修改您提供的代码以使用随机数据并填写一些来自代码示例范围之外的变量(因此它是一个独立的示例):

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata

# Generate some reproducible but random data
np.random.seed(1981)
width, height = 200, 300
x, y, z = np.random.random((3,10))
x *= width
y *= height

#create a list of x, y coordinate tuples
points = zip(x, y)

#create a grid on which to interpolate data
xi, yi = np.linspace(0, width, width), np.linspace(0, height, height)
xi, yi = np.meshgrid(xi, yi)

#interpolate the data with the matlab griddata function
zi = griddata(x, y, z, xi, yi, interp='nn')

#create a matplotlib figure and adjust the width and heights
fig = plt.figure()

#create a single subplot, just takes over the whole figure if only one is specified
ax = fig.add_subplot(111, frameon=False, xticks=[], yticks=[])

#create the contours
CS = plt.contour(xi, yi, zi, linewidths=1, colors='b')

#add a streamplot
dx, dy = np.gradient(zi)
plt.streamplot(xi, yi, dx, dy, color='c', density=1, arrowsize=3)

plt.show()

结果将如下所示: 在此处输入图像描述

请注意,有很多地方的流线不垂直于轮廓。这比箭头的错误方向更容易指示出现问题。(尽管“垂直”假设绘图的纵横比为 1,除非您设置它,否则这些绘图并不完全正确。)

解决问题

如果我们只是换行

dx, dy = np.gradient(zi)

到:

dy, dx = np.gradient(zi)

我们会得到正确的结果:

在此处输入图像描述


插值建议

另一方面,griddata在这种情况下是一个糟糕的选择。

首先,它不是一种“平滑”的插值方法。它使用 delaunay 三角测量,在三角形边界处形成“尖锐”的脊。这会导致这些位置出现异常梯度。

其次,它将插值限制为数据点的凸包,这可能是也可能不是一个好的选择。

径向基函数(或任何其他平滑插值)是插值的更好选择。

例如,如果我们修改您的代码片段以使用 RBF:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf

# Generate data
np.random.seed(1981)
width, height = 200, 300
x, y, z = np.random.random((3,10))
x *= width
y *= height

#create a grid on which to interpolate data
xi, yi = np.mgrid[0:width:1j*width, 0:height:1j*height]

#interpolate the data with the matlab griddata function
interp = Rbf(x, y, z, function='linear')
zi = interp(xi, yi)

#create a matplotlib figure and adjust the width and heights
fig, ax = plt.subplots(subplot_kw=dict(frameon=False, xticks=[], yticks=[]))

#create the contours and streamplot
CS = plt.contour(xi, yi, zi, linewidths=1, colors='b')
dy, dx = np.gradient(zi.T)
plt.streamplot(xi[:,0], yi[0,:], dx, dy, color='c', density=1, arrowsize=3)

plt.show()

在此处输入图像描述

(您会注意到,由于绘图的纵横比不相等,交叉点并不完全垂直。但是,如果我们将绘图的纵横比设置为 1,它们都是 90 度。)

作为两种方法的并排比较:

在此处输入图像描述

于 2013-06-03T17:57:00.113 回答
1

您可以使用 指定箭头样式arrowstyle='->'。尝试这两种方法,看看这是否适合您:

plt.streamplot(xi, yi, dx, dy, color='c', density=1, arrowsize=3, 
               arrowstyle='<-')

plt.streamplot(xi, yi, dx, dy, color='c', density=1, arrowsize=3, 
               arrowstyle='->')
于 2013-06-03T14:36:47.943 回答