Return xy 坐标的 z 值
Return z-value of xy coordinate
我有一组生成轮廓的 xy 坐标。对于下面的代码，这些坐标来自 df 中的组 A 和 B。我还创建了一个单独的 xy 坐标，从 C1_X 和 C1_Y 调用。然而，这并不用于生成轮廓本身。它是一个单独的xy坐标。
问题：是否可以returnC1_XC1_Y坐标处等高线的z值？
我发现了一个类似的单独问题：multivariate spline interpolation in python scipy?。该问题中的数字显示了我希望 return 但我只想要一个 xy 坐标的 z 值。
此问题中的 contour 已标准化，因此值介于 -1 和 1 之间。我希望 return C1_X 和 C1_Y 的 z 值，这是代码下方图中看到的白色散点。
我已尝试 return 此点的 z 值使用：
# Attempt at returning the z-value for C1
f = RectBivariateSpline(X, Y, normPDF)
z = f(d['C1_X'], d['C1_Y'])
print(z)
但是我 return 出错了：raise TypeError('x must be strictly increasing')
TypeError: x must be strictly increasing
我已经注释掉这个函数所以代码可以运行。
旁注：此代码是为动画编写的。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline
DATA_LIMITS = [0, 15]
def datalimits(*data):
return DATA_LIMITS
def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
XY = np.stack([X, Y], 2)
PDF = sts.multivariate_normal([x, y]).pdf(XY)
return X, Y, PDF
def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
PDFs = []
for i,(x,y) in enumerate(zip(xs,ys)):
X, Y, PDF = mvpdf(x, y, xlim, ylim)
PDFs.append(PDF)
return X, Y, np.sum(PDFs, axis=0)
fig, ax = plt.subplots(figsize = (10,6))
ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)
line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=5,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=5,zorder=3)
scat = ax.scatter([], [], s=5**2,marker='o', c='white', alpha = 1,zorder=3)
lines=[line_a,line_b]
scats=[scat]
cfs = None
def plotmvs(tdf, xlim=datalimits(df['X']), ylim=datalimits(df['Y']), fig=fig, ax=ax):
global cfs
if cfs:
for tp in cfs.collections:
tp.remove()
df = tdf[1]
PDFs = []
for (group, gdf), group_line in zip(df.groupby('group'), (line_a, line_b)):
group_line.set_data(*gdf[['X','Y']].values.T)
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, xlim, ylim)
PDFs.append(PDF)
for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
if group in ['A','B']:
group_line.set_data(*gdf[['X','Y']].values.T)
kwargs = {
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
PDFs.append(PDF)
#plot white scatter point from C1_X, C1_Y
elif group in ['C']:
gdf['X'].values, gdf['Y'].values
scat.set_offsets(gdf[['X','Y']].values)
# normalize PDF by shifting and scaling, so that the smallest value is -1 and the largest is 1
normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())
''' Attempt at returning z-value for C1_X, C1_Y '''
''' This is the function that I am trying to write that will '''
''' return the contour value '''
#f = RectBivariateSpline(X[::-1, :], Y[::-1, :], normPDF[::-1, :])
#z = f(d['C1_X'], d['C1_Y'])
#print(z)
cfs = ax.contourf(X, Y, normPDF, cmap='jet', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = fig.colorbar(cfs, ax=ax, cax=cax)
cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])
return cfs.collections + [scat] + [line_a,line_b]
''' Sample Dataframe '''
n = 1
time = range(n)
d = ({
'A1_X' : [3],
'A1_Y' : [6],
'A2_X' : [6],
'A2_Y' : [10],
'B1_X' : [12],
'B1_Y' : [2],
'B2_X' : [14],
'B2_Y' : [4],
'C1_X' : [4],
'C1_Y' : [6],
})
# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
for k,v in d.items() for i,t in enumerate(time) ]
df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']
#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)
plt.show()
我希望return白色散点的z值。预期输出将显示 C1_X、C1_Y 的规范化 z 值 (-1,1)。
根据目视检查，这将介于 0.6 和 0.8
之间
编辑 2：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline
import matplotlib.transforms as transforms
DATA_LIMITS = [-85, 85]
def datalimits(*data):
return DATA_LIMITS # dmin - spad, dmax + spad
def rot(theta):
theta = np.deg2rad(theta)
return np.array([
[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]
])
def getcov(radius=1, scale=1, theta=0):
cov = np.array([
[radius*(scale + 1), 0],
[0, radius/(scale + 1)]
])
r = rot(theta)
return r @ cov @ r.T
def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
XY = np.stack([X, Y], 2)
x,y = rot(theta) @ (velocity/2, 0) + (x, y)
cov = getcov(radius=radius, scale=scale, theta=theta)
PDF = sts.multivariate_normal([x, y], cov).pdf(XY)
return X, Y, PDF
def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
PDFs = []
for i,(x,y) in enumerate(zip(xs,ys)):
kwargs = {
'radius': radius[i] if radius is not None else 0.5,
'velocity': velocity[i] if velocity is not None else 0,
'scale': scale[i] if scale is not None else 0,
'theta': theta[i] if theta is not None else 0,
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdf(x, y,**kwargs)
PDFs.append(PDF)
return X, Y, np.sum(PDFs, axis=0)
fig, ax = plt.subplots(figsize = (10,6))
ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)
line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=3,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=3,zorder=3)
lines=[line_a,line_b] ## this is iterable!
offset = lambda p: transforms.ScaledTranslation(p/82.,0, plt.gcf().dpi_scale_trans)
trans = plt.gca().transData
scat = ax.scatter([], [], s=5,marker='o', c='white', alpha = 1,zorder=3,transform=trans+offset(+2) )
scats=[scat]
cfs = None
def plotmvs(tdf, xlim=None, ylim=None, fig=fig, ax=ax):
global cfs
if cfs:
for tp in cfs.collections:
tp.remove()
df = tdf[1]
if xlim is None: xlim = datalimits(df['X'])
if ylim is None: ylim = datalimits(df['Y'])
PDFs = []
for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
if group in ['A','B']:
group_line.set_data(*gdf[['X','Y']].values.T)
kwargs = {
'radius': gdf['Radius'].values if 'Radius' in gdf else None,
'velocity': gdf['Velocity'].values if 'Velocity' in gdf else None,
'scale': gdf['Scaling'].values if 'Scaling' in gdf else None,
'theta': gdf['Rotation'].values if 'Rotation' in gdf else None,
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
PDFs.append(PDF)
elif group in ['C']:
gdf['X'].values, gdf['Y'].values
scat.set_offsets(gdf[['X','Y']].values)
normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())
def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):
CS = ax.contour(X, Y, Z, 100)
containing_levels = []
for cc, lev in zip(CS.collections, CS.levels):
for pp in cc.get_paths():
if pp.contains_point((point_x, point_y)):
containing_levels.append(lev)
if max(containing_levels) == 0:
return 0
else:
if max(containing_levels) > 0:
lev = max(containing_levels)
adj = 1. / precision
elif max(containing_levels) < 0:
lev = min(containing_levels)
adj = -1. / precision
is_inside = True
while is_inside:
CS = ax.contour(X, Y, Z, [lev])
for pp in CS.collections[0].get_paths():
if not pp.contains_point((point_x, point_y)):
is_inside = False
if is_inside:
lev += adj
return lev - adj
print(get_contour_value_of_point(d['C1_X'], d['C1_Y'], X, Y, normPDF))
cfs = ax.contourf(X, Y, normPDF, cmap='viridis', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = fig.colorbar(cfs, ax=ax, cax=cax)
cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])
return cfs.collections + [scat] + [line_a,line_b]
''' Sample Dataframe '''
n = 10
time = range(n)
d = ({
'A1_X' : [3],
'A1_Y' : [6],
'A2_X' : [6],
'A2_Y' : [10],
'B1_X' : [12],
'B1_Y' : [2],
'B2_X' : [14],
'B2_Y' : [4],
'C1_X' : [4],
'C1_Y' : [6],
})
# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
for k,v in d.items() for i,t in enumerate(time) ]
df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']
#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)
plt.show()
如果您有一个任意的 (X, Y, Z) 点云，并且您想要对某些 (x, y) 点的 z-coordinate 进行插值，您有许多不同的选择。最简单的可能是只使用 scipy.interpolate.interp2d 来获得 z-value:
f = interp2d(X.T, Y.T, Z.T)
z = f(x, y)
由于您的网格看起来很规则，您最好使用 scipy.interpolate.RectBivariateSpline，它具有非常相似的界面，但专为规则网格而设计：
f = RectBivariateSpline(X.T, Y.T, Z.T)
z = f(x, y)
既然你有一个规则的网格，你也可以这样做
f = RectBivariateSpline(X[0, :], Y[:, 0], Z.T)
z = f(x, y)
请注意，尺寸在绘图数组和插值数组之间翻转。绘图将轴 0 视为行，即 Y，而插值函数将轴 0 视为 X。除了转置，您还可以切换 X 和 Y 输入，保持 Z 不变以获得类似的最终结果，例如：
f = RectBivariateSpline(Y, X, Z)
z = f(y, x)
或者，您也可以更改所有绘图代码以交换输入，但此时这样做工作量太大。无论你做什么，选择一种方法并坚持下去。只要你始终如一地这样做，它们应该都能奏效。
如果您使用其中一种 scipy 方法（推荐），请保持对象 f 周围以插入您可能需要的任何其他点。
如果你想要更手动的方法，你可以做一些事情，比如找到离 (x, y) 最近的三个 (X, Y, Z) 点，然后找到它们之间的平面在 (x, y）。例如：
def interp_point(x, y, X, Y, Z):
"""
x, y: scalar coordinates to interpolate at
X, Y, Z: arrays of coordinates corresponding to function
"""
X = X.ravel()
Y = Y.ravel()
Z = Z.ravel()
# distances from x, y to all X, Y points
dist = np.hypot(X - x, Y - y)
# indices of the nearest points
nearest3 = np.argpartition(dist, 2)[:3]
# extract the coordinates
points = np.stack((X[nearest3], Y[nearest3], Z[nearest3]))
# compute 2 vectors in the plane
vecs = np.diff(points, axis=0)
# compute normal to plane
plane = np.cross(vecs[0], vecs[1])
# rhs of plane equation
d = np.dot(plane, points [:, 0])
# The final result:
z = (d - np.dot(plane[:2], [x, y])) / plane[-1]
return z
print(interp_point(x, y, X.T, Y.T, Z.T))
由于您的数据位于规则网格上，因此在围绕 (x, y) 的四边形上执行类似双线性插值的操作可能更容易：
def interp_grid(x, y, X, Y, Z):
"""
x, y: scalar coordinates to interpolate at
X, Y, Z: arrays of coordinates corresponding to function
"""
X, Y = X[:, 0], Y[0, :]
# find matching element
r, c = np.searchsorted(Y, y), np.searchsorted(X, x)
if r == 0: r += 1
if c == 0: c += 1
# interpolate
z = (Z[r - 1, c - 1] * (X[c] - x) * (Y[r] - y) +
Z[r - 1, c] * (x - X[c - 1]) * (Y[r] - y) +
Z[r, c - 1] * (X[c] - x) * (y - Y[r - 1]) +
Z[r, c] * (x - X[c - 1]) * (y - Y[r - 1])
) / ((X[c] - X[c - 1]) * (Y[r] - Y[r - 1]))
return z
print(interpolate_grid(x, y, X.T, Y.T, Z.T))
这是一种不优雅的蛮力方法。*假设我们有 X、Y 和 Z 值，让我们定义一个函数，一遍又一遍地绘制自定义等高线，直到它们与 user-defined 处的点相交精度水平（在您的数据中，Z = normPDF）。
def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, 100)
containing_levels = []
for cc, lev in zip(CS.collections, CS.levels):
for pp in cc.get_paths():
if pp.contains_point((point_x, point_y)):
containing_levels.append(lev)
if max(containing_levels) == 0:
return 0
else:
if max(containing_levels) > 0:
lev = max(containing_levels)
adj = 1. / precision
elif max(containing_levels) < 0:
lev = min(containing_levels)
adj = -1. / precision
is_inside = True
while is_inside:
CS = ax.contour(X, Y, Z, [lev])
for pp in CS.collections[0].get_paths():
if not pp.contains_point((point_x, point_y)):
is_inside = False
if is_inside:
lev += adj
return lev - adj
更详细地说：这是在绘制一个具有 100 个级别的初始等高线图，然后找到其多边形包含相关点的等高线级别列表。然后我们找到最窄的级别（如果级别为正则为最高，如果级别为负则为最低）。从那里开始，我们通过小步骤（对应于您所需的精度级别）收紧级别，检查该点是否仍在多边形内。当该点不再位于轮廓多边形内时，我们知道我们找到了正确的级别（最后一个包含该点的级别）。
举个例子，我们可以使用Matplotlib库中的等高线：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2
使用此设置，get_contour_value_of_point(0, -0.6) returns 1.338399999999998 目视检查似乎匹配。 get_contour_value_of_point(0, -0.6) returns -1.48，这似乎也匹配。下图用于视觉验证。
*我不能保证这会涵盖所有用例。它涵盖了我尝试过的那些。在将其用于任何类型的生产环境之前，我会对其进行相当严格的测试。我希望有比这更优雅的解决方案（例如 ），但这是我想到的并且似乎可以直接工作，如果 brute-force，方式。
我有一组生成轮廓的 xy 坐标。对于下面的代码，这些坐标来自 df 中的组 A 和 B。我还创建了一个单独的 xy 坐标，从 C1_X 和 C1_Y 调用。然而，这并不用于生成轮廓本身。它是一个单独的xy坐标。
问题：是否可以returnC1_XC1_Y坐标处等高线的z值？
我发现了一个类似的单独问题：multivariate spline interpolation in python scipy?。该问题中的数字显示了我希望 return 但我只想要一个 xy 坐标的 z 值。
此问题中的 contour 已标准化，因此值介于 -1 和 1 之间。我希望 return C1_X 和 C1_Y 的 z 值，这是代码下方图中看到的白色散点。
我已尝试 return 此点的 z 值使用：
# Attempt at returning the z-value for C1
f = RectBivariateSpline(X, Y, normPDF)
z = f(d['C1_X'], d['C1_Y'])
print(z)
但是我 return 出错了：raise TypeError('x must be strictly increasing')
TypeError: x must be strictly increasing
我已经注释掉这个函数所以代码可以运行。
旁注：此代码是为动画编写的。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline
DATA_LIMITS = [0, 15]
def datalimits(*data):
return DATA_LIMITS
def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
XY = np.stack([X, Y], 2)
PDF = sts.multivariate_normal([x, y]).pdf(XY)
return X, Y, PDF
def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
PDFs = []
for i,(x,y) in enumerate(zip(xs,ys)):
X, Y, PDF = mvpdf(x, y, xlim, ylim)
PDFs.append(PDF)
return X, Y, np.sum(PDFs, axis=0)
fig, ax = plt.subplots(figsize = (10,6))
ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)
line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=5,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=5,zorder=3)
scat = ax.scatter([], [], s=5**2,marker='o', c='white', alpha = 1,zorder=3)
lines=[line_a,line_b]
scats=[scat]
cfs = None
def plotmvs(tdf, xlim=datalimits(df['X']), ylim=datalimits(df['Y']), fig=fig, ax=ax):
global cfs
if cfs:
for tp in cfs.collections:
tp.remove()
df = tdf[1]
PDFs = []
for (group, gdf), group_line in zip(df.groupby('group'), (line_a, line_b)):
group_line.set_data(*gdf[['X','Y']].values.T)
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, xlim, ylim)
PDFs.append(PDF)
for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
if group in ['A','B']:
group_line.set_data(*gdf[['X','Y']].values.T)
kwargs = {
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
PDFs.append(PDF)
#plot white scatter point from C1_X, C1_Y
elif group in ['C']:
gdf['X'].values, gdf['Y'].values
scat.set_offsets(gdf[['X','Y']].values)
# normalize PDF by shifting and scaling, so that the smallest value is -1 and the largest is 1
normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())
''' Attempt at returning z-value for C1_X, C1_Y '''
''' This is the function that I am trying to write that will '''
''' return the contour value '''
#f = RectBivariateSpline(X[::-1, :], Y[::-1, :], normPDF[::-1, :])
#z = f(d['C1_X'], d['C1_Y'])
#print(z)
cfs = ax.contourf(X, Y, normPDF, cmap='jet', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = fig.colorbar(cfs, ax=ax, cax=cax)
cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])
return cfs.collections + [scat] + [line_a,line_b]
''' Sample Dataframe '''
n = 1
time = range(n)
d = ({
'A1_X' : [3],
'A1_Y' : [6],
'A2_X' : [6],
'A2_Y' : [10],
'B1_X' : [12],
'B1_Y' : [2],
'B2_X' : [14],
'B2_Y' : [4],
'C1_X' : [4],
'C1_Y' : [6],
})
# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
for k,v in d.items() for i,t in enumerate(time) ]
df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']
#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)
plt.show()
我希望return白色散点的z值。预期输出将显示 C1_X、C1_Y 的规范化 z 值 (-1,1)。
根据目视检查，这将介于 0.6 和 0.8
编辑 2：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import RectBivariateSpline
import matplotlib.transforms as transforms
DATA_LIMITS = [-85, 85]
def datalimits(*data):
return DATA_LIMITS # dmin - spad, dmax + spad
def rot(theta):
theta = np.deg2rad(theta)
return np.array([
[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]
])
def getcov(radius=1, scale=1, theta=0):
cov = np.array([
[radius*(scale + 1), 0],
[0, radius/(scale + 1)]
])
r = rot(theta)
return r @ cov @ r.T
def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
XY = np.stack([X, Y], 2)
x,y = rot(theta) @ (velocity/2, 0) + (x, y)
cov = getcov(radius=radius, scale=scale, theta=theta)
PDF = sts.multivariate_normal([x, y], cov).pdf(XY)
return X, Y, PDF
def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
PDFs = []
for i,(x,y) in enumerate(zip(xs,ys)):
kwargs = {
'radius': radius[i] if radius is not None else 0.5,
'velocity': velocity[i] if velocity is not None else 0,
'scale': scale[i] if scale is not None else 0,
'theta': theta[i] if theta is not None else 0,
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdf(x, y,**kwargs)
PDFs.append(PDF)
return X, Y, np.sum(PDFs, axis=0)
fig, ax = plt.subplots(figsize = (10,6))
ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)
line_a, = ax.plot([], [], 'o', c='red', alpha = 0.5, markersize=3,zorder=3)
line_b, = ax.plot([], [], 'o', c='blue', alpha = 0.5, markersize=3,zorder=3)
lines=[line_a,line_b] ## this is iterable!
offset = lambda p: transforms.ScaledTranslation(p/82.,0, plt.gcf().dpi_scale_trans)
trans = plt.gca().transData
scat = ax.scatter([], [], s=5,marker='o', c='white', alpha = 1,zorder=3,transform=trans+offset(+2) )
scats=[scat]
cfs = None
def plotmvs(tdf, xlim=None, ylim=None, fig=fig, ax=ax):
global cfs
if cfs:
for tp in cfs.collections:
tp.remove()
df = tdf[1]
if xlim is None: xlim = datalimits(df['X'])
if ylim is None: ylim = datalimits(df['Y'])
PDFs = []
for (group, gdf), group_line in zip(df.groupby('group'), lines+scats):
if group in ['A','B']:
group_line.set_data(*gdf[['X','Y']].values.T)
kwargs = {
'radius': gdf['Radius'].values if 'Radius' in gdf else None,
'velocity': gdf['Velocity'].values if 'Velocity' in gdf else None,
'scale': gdf['Scaling'].values if 'Scaling' in gdf else None,
'theta': gdf['Rotation'].values if 'Rotation' in gdf else None,
'xlim': xlim,
'ylim': ylim
}
X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
PDFs.append(PDF)
elif group in ['C']:
gdf['X'].values, gdf['Y'].values
scat.set_offsets(gdf[['X','Y']].values)
normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())
def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):
CS = ax.contour(X, Y, Z, 100)
containing_levels = []
for cc, lev in zip(CS.collections, CS.levels):
for pp in cc.get_paths():
if pp.contains_point((point_x, point_y)):
containing_levels.append(lev)
if max(containing_levels) == 0:
return 0
else:
if max(containing_levels) > 0:
lev = max(containing_levels)
adj = 1. / precision
elif max(containing_levels) < 0:
lev = min(containing_levels)
adj = -1. / precision
is_inside = True
while is_inside:
CS = ax.contour(X, Y, Z, [lev])
for pp in CS.collections[0].get_paths():
if not pp.contains_point((point_x, point_y)):
is_inside = False
if is_inside:
lev += adj
return lev - adj
print(get_contour_value_of_point(d['C1_X'], d['C1_Y'], X, Y, normPDF))
cfs = ax.contourf(X, Y, normPDF, cmap='viridis', alpha = 1, levels=np.linspace(-1,1,10),zorder=1)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = fig.colorbar(cfs, ax=ax, cax=cax)
cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])
return cfs.collections + [scat] + [line_a,line_b]
''' Sample Dataframe '''
n = 10
time = range(n)
d = ({
'A1_X' : [3],
'A1_Y' : [6],
'A2_X' : [6],
'A2_Y' : [10],
'B1_X' : [12],
'B1_Y' : [2],
'B2_X' : [14],
'B2_Y' : [4],
'C1_X' : [4],
'C1_Y' : [6],
})
# a list of tuples of the form ((time, group_id, point_id, value_label), value)
tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
for k,v in d.items() for i,t in enumerate(time) ]
df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']
#Code will eventually operate with multiple frames
interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)
plt.show()
如果您有一个任意的 (X, Y, Z) 点云，并且您想要对某些 (x, y) 点的 z-coordinate 进行插值，您有许多不同的选择。最简单的可能是只使用 scipy.interpolate.interp2d 来获得 z-value:
f = interp2d(X.T, Y.T, Z.T)
z = f(x, y)
由于您的网格看起来很规则，您最好使用 scipy.interpolate.RectBivariateSpline，它具有非常相似的界面，但专为规则网格而设计：
f = RectBivariateSpline(X.T, Y.T, Z.T)
z = f(x, y)
既然你有一个规则的网格，你也可以这样做
f = RectBivariateSpline(X[0, :], Y[:, 0], Z.T)
z = f(x, y)
请注意，尺寸在绘图数组和插值数组之间翻转。绘图将轴 0 视为行，即 Y，而插值函数将轴 0 视为 X。除了转置，您还可以切换 X 和 Y 输入，保持 Z 不变以获得类似的最终结果，例如：
f = RectBivariateSpline(Y, X, Z)
z = f(y, x)
或者，您也可以更改所有绘图代码以交换输入，但此时这样做工作量太大。无论你做什么，选择一种方法并坚持下去。只要你始终如一地这样做，它们应该都能奏效。
如果您使用其中一种 scipy 方法（推荐），请保持对象 f 周围以插入您可能需要的任何其他点。
如果你想要更手动的方法，你可以做一些事情，比如找到离 (x, y) 最近的三个 (X, Y, Z) 点，然后找到它们之间的平面在 (x, y）。例如：
def interp_point(x, y, X, Y, Z):
"""
x, y: scalar coordinates to interpolate at
X, Y, Z: arrays of coordinates corresponding to function
"""
X = X.ravel()
Y = Y.ravel()
Z = Z.ravel()
# distances from x, y to all X, Y points
dist = np.hypot(X - x, Y - y)
# indices of the nearest points
nearest3 = np.argpartition(dist, 2)[:3]
# extract the coordinates
points = np.stack((X[nearest3], Y[nearest3], Z[nearest3]))
# compute 2 vectors in the plane
vecs = np.diff(points, axis=0)
# compute normal to plane
plane = np.cross(vecs[0], vecs[1])
# rhs of plane equation
d = np.dot(plane, points [:, 0])
# The final result:
z = (d - np.dot(plane[:2], [x, y])) / plane[-1]
return z
print(interp_point(x, y, X.T, Y.T, Z.T))
由于您的数据位于规则网格上，因此在围绕 (x, y) 的四边形上执行类似双线性插值的操作可能更容易：
def interp_grid(x, y, X, Y, Z):
"""
x, y: scalar coordinates to interpolate at
X, Y, Z: arrays of coordinates corresponding to function
"""
X, Y = X[:, 0], Y[0, :]
# find matching element
r, c = np.searchsorted(Y, y), np.searchsorted(X, x)
if r == 0: r += 1
if c == 0: c += 1
# interpolate
z = (Z[r - 1, c - 1] * (X[c] - x) * (Y[r] - y) +
Z[r - 1, c] * (x - X[c - 1]) * (Y[r] - y) +
Z[r, c - 1] * (X[c] - x) * (y - Y[r - 1]) +
Z[r, c] * (x - X[c - 1]) * (y - Y[r - 1])
) / ((X[c] - X[c - 1]) * (Y[r] - Y[r - 1]))
return z
print(interpolate_grid(x, y, X.T, Y.T, Z.T))
这是一种不优雅的蛮力方法。*假设我们有 X、Y 和 Z 值，让我们定义一个函数，一遍又一遍地绘制自定义等高线，直到它们与 user-defined 处的点相交精度水平（在您的数据中，Z = normPDF）。
def get_contour_value_of_point(point_x, point_y, X, Y, Z, precision=10000):
fig, ax = plt.subplots()
CS = ax.contour(X, Y, Z, 100)
containing_levels = []
for cc, lev in zip(CS.collections, CS.levels):
for pp in cc.get_paths():
if pp.contains_point((point_x, point_y)):
containing_levels.append(lev)
if max(containing_levels) == 0:
return 0
else:
if max(containing_levels) > 0:
lev = max(containing_levels)
adj = 1. / precision
elif max(containing_levels) < 0:
lev = min(containing_levels)
adj = -1. / precision
is_inside = True
while is_inside:
CS = ax.contour(X, Y, Z, [lev])
for pp in CS.collections[0].get_paths():
if not pp.contains_point((point_x, point_y)):
is_inside = False
if is_inside:
lev += adj
return lev - adj
更详细地说：这是在绘制一个具有 100 个级别的初始等高线图，然后找到其多边形包含相关点的等高线级别列表。然后我们找到最窄的级别（如果级别为正则为最高，如果级别为负则为最低）。从那里开始，我们通过小步骤（对应于您所需的精度级别）收紧级别，检查该点是否仍在多边形内。当该点不再位于轮廓多边形内时，我们知道我们找到了正确的级别（最后一个包含该点的级别）。
举个例子，我们可以使用Matplotlib库中的等高线：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2
使用此设置，get_contour_value_of_point(0, -0.6) returns 1.338399999999998 目视检查似乎匹配。 get_contour_value_of_point(0, -0.6) returns -1.48，这似乎也匹配。下图用于视觉验证。
*我不能保证这会涵盖所有用例。它涵盖了我尝试过的那些。在将其用于任何类型的生产环境之前，我会对其进行相当严格的测试。我希望有比这更优雅的解决方案（例如