Matplotlib outline / outline ** concave ** without data binding

I would like to make a plot of matplotlib contour

or contourf

3D data without anchor (x, y, z) that is somehow C-shaped in x and y (see sketch) - so the part of the enclosing shell around the data is concave in x and y.

I usually do graphs of 3D data without reference, by first interpolating with

from matplotlib.mlab import griddata
griddata...

      

but this creates artifacts in the concave portion of the data so that the concave portion is filled with interpolation.

Is it possible to do interpolating or contour / contour plots such that the concave portion of the data is respected?

enter image description here

+2


source to share


2 answers


Below is an example of use tricontourf

with masking to obtain a concave shape without interpolated parts outside the data. It relies on the ability to mask data based on state.

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

      



enter image description here

+3


source


The answer provided by TheImportanceOfBeingErnest gave me a perfect starting point and I applied the above code to provide a general concave hull / alpha contour plot solution. I am using Python package to create a concave body polygon. This is not a built-in function that needs to be added, which I took from this post . When we have a polygon, we just check if the mean of the triangulated points is inside the polygon and use that as our condition to form the mask.

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 formula
      area = math.sqrt(s*(s-a)*(s-b)*(s-c))
      circum_r = a*b*c/(4.0*area)
      # Here 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()

      



tricontour with alpha mask shape

0


source







All Articles