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
ListDensityPlot[surf, ColorFunction -> "Rainbow", PlotRange -> All]
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]]
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]
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}] &)
]
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]
(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]
(4) Now show the holed polygon on top of the contour plot.
Show[cplot, hole]
(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]] &}]]
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]}
Let's extract first layer of point in 3D. 218 is found by try and error.
pts = Most /@ TakeSmallestBy[xyz, Last, 218];
ListPlot[pts]
These points are not ordered.
ListLinePlot[pts]
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}]]
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]}];
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]
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]
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