Handling concavity in ListContourPlot

Here's an approach you could try, making use of ListContourPlot to build a RegionMemberFunction:

First I'll make some sample data. This should just be replaced by your actual stuff:

cutout = Region@Disk[{-2, 3}, 3];
concaveShape = RegionDifference[Disk[{0, 0}, 3], cutout];
distFunc = RegionDistance[cutout];
data = Flatten[
   Table[r*{Cos[q], Sin[q]}, {q, 0, 2 \[Pi], .1}, {r, .001, 3, .1}], 1];
data = Pick[data, RegionMember[cutout][data], False];
function = (distFunc@data);
surf = Join[data, ArrayReshape[function, {Length@function, 1}], 2];

We can see the ListDensityPlot has an issue:

concaveShape // RegionPlot

enter image description here

ListDensityPlot[surf, ColorFunction -> "Rainbow", PlotRange -> All]

enter image description here

Now here's the approach I'd advocate. First make ListContourPlot of your data to recover the shape of the cutout:

regions = ListContourPlot[surf, ColorFunction -> "Rainbow", 
  Contours -> Subdivide[0, .1, 5]]

enter image description here

Now, we'll use some post-processing to get all of the polygons in that plot. Then we'll find only the purple ones (i.e. those in the cutout). Then we'll take the RegionDifference of the RegionUnion of all of them and the RegionUnion of just the cutout:

poly =
  With[{points = regions[[1, 1]]},
   Association@
    Cases[
     regions[[1, 2]], 
     {
       ___,
       c_?ColorQ,
       ___,
       g : _GraphicsGroup | _Polygon
       } :> (c -> (g /. i_Integer :> points[[i]])), 
     Infinity
     ]
   ];

fullRegion =
  RegionUnion@Cases[Values[poly], p_Polygon :> Region[p], Infinity];

cutoutPolygons =
  With[{test = ColorData["Rainbow"][0]},
   KeySelect[poly, ColorDistance[#, test] < .1 &]
   ];

cutoutRegion =
  RegionUnion@Cases[Values[cutoutPolygons], p_Polygon :> Region[p], Infinity];

diff = RegionDifference[fullRegion, cutoutRegion]

enter image description here

Finally we can use RegionMember to get a new RegionFunction to use in ListDensityPlot:

rf = RegionMember[diff];
ListDensityPlot[
 surf, ColorFunction -> "Rainbow", 
 PlotRange -> All,
 RegionFunction -> (rf[{#, #2}] &)
 ]

enter image description here


This is a problem that I often encounter and am working on a solution. In the meantime, a simple hack that I’ve found useful is to draw a polygon with a hole on top of the contour plot. If you define the hole such that it corresponds to the boundary of the area within which you want the contours To be shown then the contours that lie outside the hole will be blanked out by the overlying polygon. The main hassle with this approach is in defining the hole geometry, especially if it’s highly irregular. There’s a function in the function repository that aims to semi-automatically define non-convex hulls for irregularly spaced points. This function generally works ok, but not for all cases so some manual intervention may still be necessary.

Here's an example...

(1) Generate some {x,y,z} data.

xyz = Table[{x, y, x*y}, {x, 0, 10}, {y, 0, 10}] // Flatten[#, 1] &;

(2) Plot the data.

cplot=ListContourPlot[xyz]

enter image description here

(3) Define polygon with a hole and plot it.

hole = Graphics[{EdgeForm[{Gray, Thickness[Tiny]}], FaceForm[White], Polygon[{{-1, -1}, {11, -1}, {11, 11}, {-1, 11}} -> {{1, 1}, {5, 
   5}, {9, 1}, {5, 9}}]}, Frame -> True]

enter image description here

(4) Now show the holed polygon on top of the contour plot.

Show[cplot, hole]

enter image description here

(5) If you wanted to, you could also play some tunes with FaceForm Opacity to feintly show the blanked-out contours, e.g.

Show[ListContourPlot[xyz],Graphics[{EdgeForm[{Gray, Thickness[Tiny]}], FaceForm[{White, Opacity[0.8]}],Polygon[{{-1, -1}, {11, -1}, {11, 11}, {-1, 11}} -> {{1, 1}, {5,5}, {9, 1}, {5, 9}}] // Rotate[#, 270 \[Degree]] &}]]

enter image description here

Hope this helps.


Here is my solution based on @lan's idea.

Let xyz 1st, 2nd and 6th column of data.

Let's plot it in 2D and 3D.

  {ListPlot[Most /@ xyz, Frame -> True, Axes -> False], ListPointPlot3D[xyz]}

enter image description here

Let's extract first layer of point in 3D. 218 is found by try and error.

pts = Most /@ TakeSmallestBy[xyz, Last, 218];
ListPlot[pts]

enter image description here

These points are not ordered.

ListLinePlot[pts]

enter image description here

Let's order them counterclockwise.

q1 = ReverseSortBy[Select[pts, #[[1]] > 0 && #[[2]] > 0 &], ArcTan[#2/#1 &]];
q2 = ReverseSortBy[Select[pts, #[[1]] < 0 && #[[2]] > 0 &], ArcTan[#2/#1 &]];
q3 = SortBy[Select[pts, #[[1]] < 0 && #[[2]] < 0 &], ArcTan[#2/#1 &]];
q4 = SortBy[Select[pts, #[[1]] > 0 && #[[2]] < 0 &], ArcTan[#2/#1 &]];

ListLinePlot[Join[q1, q2, q3, q4, {First@q1}]]

enter image description here

Now we can use these points as a mask. I chose {{1, 0}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 0}} as frame coordinates.

        mask = Join[q1, q2, q3,  q4, {First@q1}, {{1, 0}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 0}}];
        hole = Graphics[{FaceForm[White], EdgeForm[{Gray, Thickness[Tiny]}], 
        Polygon[mask]}];

enter image description here

Replace EdgeForm[{Gray, Thickness[Tiny]}] with EdgeForm[ColorData["Rainbow"] /@ Subdivide[5] // First]

And all together, it gives

   pts = Most /@ TakeSmallestBy[xyz, Last, 218];
q1 = ReverseSortBy[Select[pts, #[[1]] > 0 && #[[2]] > 0 &], 
   ArcTan[#2/#1 &]];
q2 = ReverseSortBy[Select[pts, #[[1]] < 0 && #[[2]] > 0 &], 
   ArcTan[#2/#1 &]];
q3 = SortBy[Select[pts, #[[1]] < 0 && #[[2]] < 0 &], ArcTan[#2/#1 &]];
q4 = SortBy[Select[pts, #[[1]] > 0 && #[[2]] < 0 &], ArcTan[#2/#1 &]];
mask = Join[q1, q2, q3, 
   q4, {First@q1}, {{1, 0}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 0}}];
hole = Graphics[{FaceForm[White], 
    EdgeForm[ColorData["Rainbow"] /@ Subdivide[5] // First], 
    Polygon[mask]}];
Show[ListDensityPlot[xyz, ColorFunction -> "Rainbow", 
  AspectRatio -> Automatic, ImageSize -> 600], hole]

enter image description here

Alternatively we can use FindCurvePath to sort the points (Or FindShortestTour[pts] will also work, i.e., orderedPts =pts[[Last@FindShortestTour[pts]]]).

curve = First@FindCurvePath[pts];
orderedPts = pts[[curve]];
ListLinePlot[orderedPts]

enter image description here

mask = Join[orderedPts, {{-0.005, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 1}, {-0.005, 1}}];
hole = Graphics[{FaceForm[White], 
    EdgeForm[ColorData["Rainbow"] /@ Subdivide[5] // First], 
    Polygon[mask]}];
Show[ListDensityPlot[xyz, ColorFunction -> "Rainbow", 
  AspectRatio -> Automatic, ImageSize -> 600], hole]

Same picture