**凹面**非網格數據的matplotlib輪廓/ contourf

Question

我想做一個matplotlib contour或contourf非x格式的3D數據（x，y，z）和y（見草圖） - 因此，圍繞數據的封閉船體的一部分在x和y中是凹的。

一般我通過首先用

from matplotlib.mlab import griddata 
griddata... 

插值它做非網格化的三維數據的曲線圖，但是這產生在數據，使得所述凹部是由內插填充的凹部的僞影。

是否有可能進行插值或contourf /等高線圖，使數據的凹陷部分得到尊重？

Answer 1

下面是關於如何使用tricontourf與掩蔽，以獲得凹形狀而不數據外插部的例子。它依賴於根據條件掩蓋數據的能力。

import matplotlib.pyplot as plt 
import matplotlib.tri as tri 
import numpy as np 

# create some data 
rawx = np.random.rand(500) 
rawy = np.random.rand(len(rawx)) 
cond01 = (rawx-1)**2 + rawy**2 <=1 
cond02 = (rawx-0.7)**2 + rawy**2 >0.3 
x = rawx[cond01 & cond02] 
y = rawy[cond01 & cond02] 
f = lambda x,y: np.sin(x*4)+np.cos(y) 
z = f(x,y) 
# now, x,y are points within a partially concave shape 

triang0 = tri.Triangulation(x, y) 
triang = tri.Triangulation(x, y) 
x2 = x[triang.triangles].mean(axis=1) 
y2 = y[triang.triangles].mean(axis=1) 
#note the very obscure mean command, which, if not present causes an error. 
#now we need some masking condition. 
# this is easy in this case where we generated the data according to the same condition 
cond1 = (x2-1)**2 + y2**2 <=1 
cond2 = (x2-0.7)**2 + (y2)**2 >0.3 
mask = np.where(cond1 & cond2,0,1) 
# apply masking 
triang.set_mask(mask) 


fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(6,3)) 
ax.set_aspect("equal") 
ax2.set_aspect("equal") 

ax.tricontourf(triang0, z, cmap="Oranges") 
ax.scatter(x,y, s=3, color="k") 

ax2.tricontourf(triang, z, cmap="Oranges") 
ax2.scatter(x,y, s=3, color="k") 

ax.set_title("tricontourf without mask") 
ax2.set_title("tricontourf with mask") 
ax.set_xlim(0,1) 
ax.set_ylim(0,1) 
ax2.set_xlim(0,1) 
ax2.set_ylim(0,1) 

plt.show() 

Answer 2

由TheImportanceOfBeingErnest提供的答案給我的理想起點和我已經操縱上面的代碼，以提供到凹殼/α-形狀的等高線圖的一般解決方案。我使用Python軟件包來創建一個凹形的多邊形。這不是一個需要添加的內置函數，我從this post中獲得了這個函數。一旦我們有了多邊形，我們只需檢查三角點的平均值是否在多邊形內，並將其用作形成掩模的條件。

import matplotlib.pyplot as plt 
import matplotlib.tri as tri 
import numpy as np 
# Start Using SHAPELY 
import shapely.geometry as geometry 
from shapely.geometry import Polygon, MultiPoint, Point 
from shapely.ops import triangulate 
from shapely.ops import cascaded_union, polygonize 
from scipy.spatial import Delaunay 

from descartes.patch import PolygonPatch 
import math 

# create some data 
rawx = np.random.rand(500) 
rawy = np.random.rand(len(rawx)) 
cond01 = (rawx-1)**2 + rawy**2 <=1 
cond02 = (rawx-0.7)**2 + rawy**2 >0.3 
x = rawx[cond01 & cond02] 
y = rawy[cond01 & cond02] 
f = lambda x,y: np.sin(x*4)+np.cos(y) 
z = f(x,y) 
# now, x,y are points within a partially concave shape 

triang0 = tri.Triangulation(x, y) 
triang = tri.Triangulation(x, y) 
# Function for finding an alpha function 
def alpha_shape(points, alpha): 
    """ 
    Compute the alpha shape (concave hull) of a set 
    of points. 
    @param points: Iterable container of points. 
    @param alpha: alpha value to influence the 
     gooeyness of the border. Smaller numbers 
     don't fall inward as much as larger numbers. 
     Too large, and you lose everything! 
    """ 
    if len(points) < 4: 
    # When you have a triangle, there is no sense 
    # in computing an alpha shape. 
    return geometry.MultiPoint(list(points)).convex_hull 
    def add_edge(edges, edge_points, coords, i, j): 
    """ 
    Add a line between the i-th and j-th points, 
    if not in the list already 
    """ 
    if (i, j) in edges or (j, i) in edges: 
     # already added 
     return 
    edges.add((i, j)) 
    edge_points.append(coords[ [i, j] ]) 

    coords = np.array([point.coords[0] 
        for point in points]) 

tri = Delaunay(coords) 
    edges = set() 
    edge_points = [] 
    # loop over triangles: 
    # ia, ib, ic = indices of corner points of the 
    # triangle 
    for ia, ib, ic in tri.vertices: 
     pa = coords[ia] 
     pb = coords[ib] 
     pc = coords[ic] 
     # Lengths of sides of triangle 
     a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2) 
     b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2) 
     c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2) 
     # Semiperimeter of triangle 
     s = (a + b + c)/2.0 
     # Area of triangle by Heron's formula 
     area = math.sqrt(s*(s-a)*(s-b)*(s-c)) 
     circum_r = a*b*c/(4.0*area) 
     # Here's the radius filter. 
     #print circum_r 
     if circum_r < 1.0/alpha: 
      add_edge(edges, edge_points, coords, ia, ib) 
      add_edge(edges, edge_points, coords, ib, ic) 
      add_edge(edges, edge_points, coords, ic, ia) 
    m = geometry.MultiLineString(edge_points) 
    triangles = list(polygonize(m)) 
    #import ipdb; ipdb.set_trace() 
    return cascaded_union(triangles), edge_points 

# create array of points from reduced exp data to convert to Polygon 
crds=np.array([x,y]).transpose() 
# Adjust the length of acceptable sides by adjusting the alpha parameter 
concave_hull, edge_points = alpha_shape(MultiPoint(crds), alpha=2.3) 

# View the polygon and adjust alpha if needed 
def plot_polygon(polygon): 
    fig = plt.figure(figsize=(10,10)) 
    ax = fig.add_subplot(111) 
    margin = .3 
    x_min, y_min, x_max, y_max = polygon.bounds 
    ax.set_xlim([x_min-margin, x_max+margin]) 
    ax.set_ylim([y_min-margin, y_max+margin]) 
    patch = PolygonPatch(polygon, fc='#999999', 
         ec='#000000', fill=True, 
         zorder=-1) 
    ax.add_patch(patch) 
    return fig 
plot_polygon(concave_hull); plt.plot(x,y,'.'); #plt.show() 


# Use the mean distance between the triangulated x & y poitns 
x2 = x[triang.triangles].mean(axis=1) 
y2 = y[triang.triangles].mean(axis=1) 
##note the very obscure mean command, which, if not present causes an error. 
##now we need some masking condition. 

# Create an empty set to fill with zeros and ones 
cond = np.empty(len(x2)) 
# iterate through points checking if the point lies within the polygon 
for i in range(len(x2)): 
    cond[i] = concave_hull.contains(Point(x2[i],y2[i])) 

mask = np.where(cond,0,1) 
# apply masking 
triang.set_mask(mask) 

fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(6,3)) 
ax.set_aspect("equal") 
ax2.set_aspect("equal") 

ax.tricontourf(triang0, z, cmap="Oranges") 
ax.scatter(x,y, s=3, color="k") 

ax2.tricontourf(triang, z, cmap="Oranges") 
ax2.scatter(x,y, s=3, color="k") 

ax.set_title("tricontourf without mask") 
ax2.set_title("tricontourf with mask") 
ax.set_xlim(0,1) 
ax.set_ylim(0,1) 
ax2.set_xlim(0,1) 
ax2.set_ylim(0,1) 

plt.show()