matplotlibの輪郭/輪郭**凹形の非グリッド化データ

Question

私は何らかの形でC字型の非グリッド化された3Dデータ（x、y、z）のプロットをmatplotlib contourまたはcontourfしたいと思いますとy（スケッチを参照） - したがって、データの周りの包囲する船体の一部はxとyで凹です。

通常私は最初

from matplotlib.mlab import griddata 
griddata... 

とそれを補間することにより、非グリッド3Dデータのプロットを行うが、これは凹部が補間によって充填されるようにデータの凹部にアーティファクトを生成します。

データの凹み部分が尊重されるように、補間またはcontourf/contourプロットを行うことはできますか？以下

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()