【问题标题】:surface plot on irregular grid in python environmentpython环境中不规则网格上的曲面图
【发布时间】:2018-08-24 15:58:04
【问题描述】:

是否可以创建函数的曲面图和等高线图,例如u(x,y) = x^2 + y^2

在以下由方程界定的域上

r(t) = 1+(cos(4*t))^2, x = r(t)*cos(t),y = r(t)*sin(t), 0 < t < 2*pi

我想要以下散点图的曲面图变体。 我也使用了 scipy griddata 如下

from matplotlib.pyplot  import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import griddata

data = np.array([[ 2.00000000e+00,  0.00000000e+00],
   [ 1.89614525e+00,  1.51126793e-01],
   [ 1.62613327e+00,  2.60869688e-01],
   [ 1.29639187e+00,  3.15328472e-01],
   [ 1.03183519e+00,  3.39805963e-01],
   [ 9.22330309e-01,  3.87424263e-01],
   [ 9.85968247e-01,  5.09806830e-01],
   [ 1.16496791e+00,  7.25092951e-01],
   [ 1.35481653e+00,  1.00089121e+00],
   [ 1.45215333e+00,  1.26290100e+00],
   [ 1.39897768e+00,  1.42707450e+00],
   [ 1.20351374e+00,  1.44074409e+00],
   [ 9.29854857e-01,  1.31248166e+00],
   [ 6.63453200e-01,  1.11474417e+00],
   [ 4.70396020e-01,  9.55864325e-01],
   [ 3.70174853e-01,  9.32786757e-01],
   [ 3.33831717e-01,  1.08590518e+00],
   [ 3.06157327e-01,  1.37738581e+00],
   [ 2.38732323e-01,  1.70413196e+00],
   [ 1.15900227e-01,  1.94068353e+00],
   [-3.96388724e-02,  1.99329358e+00],
   [-1.83621616e-01,  1.84088591e+00],
   [-2.79344160e-01,  1.54430632e+00],
   [-3.22477316e-01,  1.21965478e+00],
   [-3.47231849e-01,  9.87829012e-01],
   [-4.09520279e-01,  9.23088740e-01],
   [-5.55288737e-01,  1.02352335e+00],
   [-7.90735832e-01,  1.21586141e+00],
   [-1.07111976e+00,  1.39113164e+00],
   [-1.31601256e+00,  1.45381050e+00],
   [-1.44544698e+00,  1.36169834e+00],
   [-1.42013085e+00,  1.13913085e+00],
   [-1.26561597e+00,  8.59376173e-01],
   [-1.06700598e+00,  6.06642746e-01],
   [-9.34540362e-01,  4.37040778e-01],
   [-9.54661396e-01,  3.57053472e-01],
   [-1.14896223e+00,  3.28361061e-01],
   [-1.46070082e+00,  2.94327457e-01],
   [-1.77632959e+00,  2.12929020e-01],
   [-1.97334870e+00,  7.85155418e-02],
   [-1.97334870e+00, -7.85155418e-02],
   [-1.77632959e+00, -2.12929020e-01],
   [-1.46070082e+00, -2.94327457e-01],
   [-1.14896223e+00, -3.28361061e-01],
   [-9.54661396e-01, -3.57053472e-01],
   [-9.34540362e-01, -4.37040778e-01],
   [-1.06700598e+00, -6.06642746e-01],
   [-1.26561597e+00, -8.59376173e-01],
   [-1.42013085e+00, -1.13913085e+00],
   [-1.44544698e+00, -1.36169834e+00],
   [-1.31601256e+00, -1.45381050e+00],
   [-1.07111976e+00, -1.39113164e+00],
   [-7.90735832e-01, -1.21586141e+00],
   [-5.55288737e-01, -1.02352335e+00],
   [-4.09520279e-01, -9.23088740e-01],
   [-3.47231849e-01, -9.87829012e-01],
   [-3.22477316e-01, -1.21965478e+00],
   [-2.79344160e-01, -1.54430632e+00],
   [-1.83621616e-01, -1.84088591e+00],
   [-3.96388724e-02, -1.99329358e+00],
   [ 1.15900227e-01, -1.94068353e+00],
   [ 2.38732323e-01, -1.70413196e+00],
   [ 3.06157327e-01, -1.37738581e+00],
   [ 3.33831717e-01, -1.08590518e+00],
   [ 3.70174853e-01, -9.32786757e-01],
   [ 4.70396020e-01, -9.55864325e-01],
   [ 6.63453200e-01, -1.11474417e+00],
   [ 9.29854857e-01, -1.31248166e+00],
   [ 1.20351374e+00, -1.44074409e+00],
   [ 1.39897768e+00, -1.42707450e+00],
   [ 1.45215333e+00, -1.26290100e+00],
   [ 1.35481653e+00, -1.00089121e+00],
   [ 1.16496791e+00, -7.25092951e-01],
   [ 9.85968247e-01, -5.09806830e-01],
   [ 9.22330309e-01, -3.87424263e-01],
   [ 1.03183519e+00, -3.39805963e-01],
   [ 1.29639187e+00, -3.15328472e-01],
   [ 1.62613327e+00, -2.60869688e-01],
   [ 1.89614525e+00, -1.51126793e-01],
   [ 2.00000000e+00, -4.89858720e-16],
   [ 0.00000000e+00, -1.50000000e+00],
   [-1.00000000e+00, -1.00000000e+00],
   [ 0.00000000e+00, -1.00000000e+00],
   [ 1.00000000e+00, -1.00000000e+00],
   [-5.00000000e-01, -5.00000000e-01],
   [ 0.00000000e+00, -5.00000000e-01],
   [ 5.00000000e-01, -5.00000000e-01],
   [-1.50000000e+00,  0.00000000e+00],
   [-1.00000000e+00,  0.00000000e+00],
   [-5.00000000e-01,  0.00000000e+00],
   [ 0.00000000e+00,  0.00000000e+00],
   [ 5.00000000e-01,  0.00000000e+00],
   [ 1.00000000e+00,  0.00000000e+00],
   [ 1.50000000e+00,  0.00000000e+00],
   [ 2.00000000e+00,  0.00000000e+00],
   [-5.00000000e-01,  5.00000000e-01],
   [ 0.00000000e+00,  5.00000000e-01],
   [ 5.00000000e-01,  5.00000000e-01],
   [-1.00000000e+00,  1.00000000e+00],
   [ 0.00000000e+00,  1.00000000e+00],
   [ 1.00000000e+00,  1.00000000e+00],
   [ 0.00000000e+00,  1.50000000e+00]])   
ua = data[:,0]**2+data[:,1]**2 # u=x^2+y^2


xx,yy = np.meshgrid(np.linspace(-2,2,100),np.linspace(-2,2,100))
Ua = griddata((data[:,0],data[:,1]),ua,(xx,yy),method='cubic') 

fig = figure(1)
plot (data[:,0], data[:,1], '*'); # 
fig = figure(2)
ax = fig.gca(projection='3d')

ax.plot_wireframe(xx,yy,Ua,rstride=1,cstride=1,linewidth=.5) 

但结果并不好,如下所示

【问题讨论】:

  • 你用什么代码来创建第二个情节?你希望它看起来像什么?
  • @Joooeey:我添加了代码。正如我所说,我希望给定域的表面图和等高线图。
  • 好的。结果有什么问题?它是否显示正确的域?看起来该域可能已关闭。此外,如果您提供一个 Minimal, Complete, and Verifiable example (MCVE),我可以将其复制到我的 IDE 中,这将更容易找出问题所在。即,定义data的代码在哪里?
  • 你展示的是一个曲面图。您可以使用contour(xx, yy, Ua) 获取等高线图。
  • 我知道这是一个曲面图,但它不正确。 (见散点图)。

标签: python matplotlib plot


【解决方案1】:

首先在矩形网格上计算函数!

scipy.interpolate.griddata 在凸包上进行插值。您可以在the docs 中阅读。这意味着它还会在您的叶之间进行插值。这就是为什么你没有得到正确的域。插值在本质上也不如仅计算跨网格的函数准确。

plot_wireframe 想要一个您已经创建的矩形网格。您需要做的就是计算矩形网格上的函数值。要仅在域边界上绘制值,请将其之外的所有内容设置为 np.NaN(不是数字)。

这是怎么做的:

import matplotlib.pyplot as plt
import numpy as np

# cartesian coordinates
xx,yy = np.meshgrid(np.linspace(-2,2,100),np.linspace(-2,2,100))

# function value across square domain
Ua = xx**2 + yy**2

# polar coordinates
tt = np.arctan2(yy, xx)
rr = np.sqrt(Ua) # re-using x^2 + y^2 -- only works for this function

# r coordinate of domain boundary
domain_boundary = 1 + (np.cos(4*tt))**2

# function value across actual domain, with rest set to NaN
Ua[rr > domain_boundary] = np.NaN

# plotting
fig = plt.figure(2)
ax = fig.gca(projection='3d')
ax.plot_wireframe(xx,yy,Ua,rstride=1,cstride=1,linewidth=.5)

该解决方案并不完美,因为您可以识别结果中的矩形网格。您可以尝试在极坐标中工作,如 this official matplotlib example 所示。

【讨论】:

  • 其实我的Ua一般没有x^2+y^2这样的公式。如果我的 Ua 不能明确说明怎么办。上述代码也不适用于 plot_surface 命令。
  • 1.如果 Ua 有不同的公式,您只需替换 rr = np.sqrt(xx**2 + yy**2) 即可获得半径坐标。 2. 只要 Ua 是 2D 域上的函数,就应该有一种方法可以明确地说明它(即从 x 和 y 值计算它)。如果它不是一个函数(即,如果一对 xy 有两个 z 值),事情就会变得复杂。最好为此提出一个新问题。
  • ax.plot_surface 对我来说不会失败。域边界看起来非常崎岖。您必须花很多心思来制作网格以使其看起来不错。您也可以尝试ax.plot_trisurf(如已删除的 Q 答案中所建议的那样)。 Trisurf 将 xx、yy 和 Ua 作为一维数组。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2011-06-22
  • 2011-12-03
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2018-08-14
相关资源
最近更新 更多