边界插值蒙版

Question

bounty到期在6天内。这个问题的答案有资格获得+ 50声望赏金。 Tomas Jacobsen想要引起更多关注对此问题。

下图显示了具有X-Y坐标值和Z值的数据点。我想生成Z值之间的梯度的热图，并使用Plotly显示。

我遇到的问题是能够有效地掩盖在没有数据点的'凹形'区域中通过插值形成的有害数据。

这些是我到目前为止所做的一些尝试：

1. 我可以通过在Delaunay三角网格上执行插值来使用Matplotlib达到理想的结果，但是我不确定如何在Plotly中复制它。

2. 我已经尝试了一种KDTree方法，尝试将遮罩应用于特定距离内的点，但是，这并未在所需的边界“线”处应用遮罩。

下面是我的代码，可以在here中找到示例输入数据：

from matplotlib.tri import (Triangulation, UniformTriRefiner, CubicTriInterpolator)
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import plotly.graph_objs as go
from plotly.subplots import make_subplots
from scipy.spatial import KDTree


#-----------------------------------------------------------------------------
# Import data
#-----------------------------------------------------------------------------
df = pd.read_excel(
    r'C:\Users\morga\PycharmProjects\Sample_data'
    r'.xlsx')


x = df['X'].to_numpy()
y = df['Y'].to_numpy()
Z_value = df['Delay'].to_numpy()

#-----------------------------------------------------------------------------
# Creating a Triangulation
#-----------------------------------------------------------------------------
triang = Triangulation(x, y)

# Masking condition for unwanted triangles which have formed between concave areas.
max_radius = 10 #I need to make this dynamic
triangles = triang.triangles

#Mask off unwanted triangles by max_radius
xtri = x[triangles] - np.roll(x[triangles], 1, axis=1)
ytri = y[triangles] - np.roll(y[triangles], 1, axis=1)
maxi = np.max(np.sqrt(xtri**2 + ytri**2), axis=1)
triang.set_mask(maxi > max_radius)
masked_triang = triang
print(masked_triang)
print(triang.x)

#-----------------------------------------------------------------------------
# Plot masked triangulation mesh
#-----------------------------------------------------------------------------
fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')
ax1.triplot(masked_triang, 'bo-', lw=0.2)
plt.show()

#-----------------------------------------------------------------------------
# Calculate the gradient dx,dy at each mesh node point
#-----------------------------------------------------------------------------
tci = CubicTriInterpolator(masked_triang, Z_value)
# Gradient requested here at the mesh nodes:
(dx, dy) = tci.gradient(triang.x, triang.y)
Resultant_BR = np.sqrt(dx**2 + dy**2)


print("BR max:",np.nanmax(Resultant_BR),"BR min:",np.nanmin(Resultant_BR),"BR mean:",np.nanmean(Resultant_BR))
Resultant_BR= np.array(Resultant_BR,dtype=float)

#-----------------------------------------------------------------------------
# Plot the Resultant_BR vectors (Gradient) at each of the mesh nodepoints
#-----------------------------------------------------------------------------
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Enforce the margins, and enlarge them to give room for the vectors.
ax.use_sticky_edges = False
ax.margins(0.07)
ax.triplot(masked_triang, color='0.8')
ax.scatter(x,y,s=5)

# Plots direction and magnitude of the gradients at the mesh nodepoints
ax.quiver(x, y, dx, dy, color='blue')
ax.set_title('Gradient vectors at meshnode points')
plt.show()


#------------------ Plot Gradient on Masked Triangle mesh ---------------

fig1, ax1 = plt.subplots()
ax1.tricontourf(masked_triang, Resultant_BR,  cmap="Oranges")
ax1.scatter(x,y, s=3, color="k")
ax1.set(xlim=(min(x)-10,max(x)+10), ylim=(min(y)-10,max(y)+10), aspect="equal")
plt.show()



subset = df[['X', 'Y']]
List1 = [tuple(x) for x in subset.to_numpy()]  
points = np.array(List1)  # change list to array x-y
values = Resultant_BR  # z values as a numpy array size(n,)

x=df['X'].to_numpy()
y=df['Y'].to_numpy()
z=values

# create meshgrid
grid_x, grid_y = np.meshgrid((np.arange(min(x), max(x),1 )), (np.arange(min(y), max(y),1)))
grid_z = griddata(points, values, (grid_x, grid_y), method='linear')
print("grid_z",grid_z)


#------------Plotly Heatmap ---------------------------------------------------------------

fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=x,
        y=y,mode="markers"
    ))

fig.add_trace(
    go.Heatmap(z=grid_z,x0=min(x),y0=min(y),showscale=True, zsmooth='best', connectgaps=False,
               colorscale='Hot'
               )
)

fig.update_layout(
    width = 1200,
    height = 1200,
    title = "Gradient Heatmap Plot",
    yaxis = dict(
      scaleanchor = "x",
      scaleratio = 1,
    ))

fig.show()

#-----------------------------------------------------------------------------
# Apply KDTree to mask by distance to neighbours
#-----------------------------------------------------------------------------
grid_x, grid_y = np.meshgrid((np.arange(min(x), max(x),1 )), (np.arange(min(y), max(y),1)))
grid_z = griddata(points, values, (grid_x, grid_y), method='linear')
tree = KDTree(np.c_[points])
dist, _ = tree.query(np.c_[grid_x.ravel(), grid_y.ravel()], k=1)
dist = dist.reshape(grid_x.shape)
grid_z[dist > 6] = np.nan #Insert NaN values at Distance greater than 6

#------------Plotly Heatmap ---------------------------------------------------------------
fig = go.Figure()

fig.add_trace(
    go.Scatter(
        x=x,
        y=y,mode="markers"
    ))

fig.add_trace(
    go.Heatmap(z=grid_z,x0=min(x),y0=min(y),showscale=True, zsmooth='best', connectgaps=False,
               colorscale='Hot'
               )
)

fig.update_layout(
    width = 1200,
    height = 1200,
    title = "Gradient Heatmap Plot",
    yaxis = dict(
      scaleanchor = "x",
      scaleratio = 1,
    ))

fig.show()

不幸的是，我已经花了很多时间试图解决这个问题，但是没有运气。您的帮助将不胜感激。

Answer 1

一种解决方法（也许不是最优雅的方法是找到点的边界（凹壳），然后将边界之外的任何内容设置为nan。

要找到边界，可以使用alphashape并确定grid_z点是否在边界内（或在边界上，可以使用shapely。

这是在第一个情节剧情开始之前接手的示例：

from shapely.geometry import Polygon, Point
import alphashape

mpoints = [Point(X, Y) for X, Y in zip(x, y)]

alpha=.125
hull = alphashape.alphashape(mpoints, alpha)
poly = Polygon(hull)

grid_gz = grid_z
gx = np.arange(min(x), max(x),1)
gy = np.arange(min(y), max(y),1)
for i, gxi in enumerate(gx):
    for j, gyi in enumerate(gy):
        p1 = Point(gxi, gyi)
        test = poly.contains(p1) | poly.touches(p1)
        if test==False:
            grid_gz[j,i]=np.nan


fig = go.Figure()

fig.add_trace(
    go.Heatmap(z=grid_gz,x0=min(x),y0=min(y),showscale=True, zsmooth='best', connectgaps=False,
               colorscale='Hot'
               )
)

fig.add_trace(
    go.Scatter(
        x=x,
        y=y,mode="markers",marker_size=2,marker_color="black",
    ))

fig.update_layout(
    width = 1200,
    height = 1200,
    title = "Gradient Heatmap Plot",
    yaxis = dict(
      scaleanchor = "x",
      scaleratio = 1,
    ))

fig.show()

进一步的评论：

• 我最终手动选择了alpha值。您可以使用optimizealpha查找Alpha值，但是该值最终还是有点偏，并且一条边在热图中有台阶。另请参见here。
• 可以在here上找到alphashape的更多信息。
• 可以在形状的here上找到contains，在here上找到touches更多信息，>