6
$\begingroup$

I have imported some CFD results (in csv format) into Mathematica to do some further processing of the data.

One of my 2D geometries is,

enter image description here

Now, when I try to do the contour plot, I get the following:

enter image description here

Why is contour plot filling those areas above and below the narrow section with purplish color. Those regions are not even within my X and Y coordinates (see the image for geometry above)?

Any ideas?

Thanks

Also cross-posted here.

$\endgroup$
9
  • 1
    $\begingroup$ I suspect the issue is due to the concave nature of your shape. Finding concave boundaries is much trickier than finding convex shapes. If you could find the boundary of your shape as a region, you might be able to use the RegionFunction argument of ContourPlot to clip it to the region of interest only. $\endgroup$
    – MassDefect
    Dec 10, 2019 at 23:20
  • $\begingroup$ The ListPlot shows an uneven density of points from which the features emerge. The ListContourPlot has a z value that is not necessarily associated with the density of the points. Without being able to see your data, it's certainly plausible that the two figures should not look similar (other than maybe having the same outer boundary). What can you show to counter that inference? $\endgroup$
    – JimB
    Dec 11, 2019 at 0:04
  • $\begingroup$ @JimB The ListPlot is showing a csv file containing x and y coordinates only. Then I have another csv file which contains the same x and y cooridinates as the first and second column and a velocity as the third column. So ListContourPlot is used to show the velocity contour. I am not sure if I answered your question correclty. Please let me know. $\endgroup$ Dec 11, 2019 at 1:13
  • $\begingroup$ @MassDefect Unfortunately the shape is random and I cannot specify it mathematically. $\endgroup$ Dec 11, 2019 at 1:14
  • 1
    $\begingroup$ Yes, you did answer my concern. I mistakenly assumed your issue was about the difference in the "patterns" observed in the two figures rather than the first figure simply demonstrating the region for which you wanted a dense contour plot. $\endgroup$
    – JimB
    Dec 11, 2019 at 16:54

4 Answers 4

8
$\begingroup$

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

$\endgroup$
12
  • $\begingroup$ Thanks a lot for all the time you spent on this. I honestly did not think it would be such a difficult thing to do. $\endgroup$ Dec 11, 2019 at 2:21
  • $\begingroup$ My sense is that it's a non-trivial problem. There might be a better approach, but if there is I don't know it $\endgroup$
    – b3m2a1
    Dec 11, 2019 at 6:00
  • 1
    $\begingroup$ Just did it. Thanks $\endgroup$ Dec 11, 2019 at 17:53
  • 1
    $\begingroup$ Since you are effectively using the z-values to get the cutout, shouldn't you be able to just use ListDensityPlot[surf, ColorFunction -> "Rainbow", PlotRange -> {0.1, All}]? $\endgroup$
    – Lukas Lang
    Dec 12, 2019 at 9:14
  • 1
    $\begingroup$ Good point - but I guess in those cases you could just use RegionFunction -> (!zmin < #3 < zmax&) or similar, without the need to extract the region. (Please don't take this as criticism of your solution - in fact, the way you extract the region from the plot is really nice :) ) $\endgroup$
    – Lukas Lang
    Dec 12, 2019 at 17:11
3
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Could you add a link to the function repo page? That's potentially very useful function. Also would it be possible to provide some example code to detail your proposed method? $\endgroup$
    – b3m2a1
    Dec 12, 2019 at 2:18
  • 1
    $\begingroup$ I’m not at my computer today, so can’t access any example code. But here’s the link to Jon McLoon’s function - resources.wolframcloud.com/FunctionRepository/resources/…. And, for plotting polygons with holes see the documentation pages at reference.wolfram.com/language/ref/Polygon.html. Just draw the holed polygon on top of your contour plot and hey presto. The contours are still there, but shouldn’t be visible because they are blanked-out by the overlying polygon. $\endgroup$
    – Ian
    Dec 12, 2019 at 7:39
  • $\begingroup$ (when you get the chance) it’d be much more helpful if you could actually add some code and examples to the post. More useful for posterity, you know? $\endgroup$
    – b3m2a1
    Dec 12, 2019 at 7:40
  • $\begingroup$ Agreed, will provide an example when I have a mo. $\endgroup$
    – Ian
    Dec 12, 2019 at 7:45
3
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ you can use pts with ListCurvePathPlot : Show[ListCurvePathPlot[pts, PlotStyle -> Red], ListPlot[xyz[[All, ;; 2]]], Axes -> False] $\endgroup$
    – kglr
    Dec 14, 2019 at 2:20
  • $\begingroup$ Good point, but ordering is necessary for Polygon, no? $\endgroup$ Dec 14, 2019 at 3:15
  • 1
    $\begingroup$ But we can use FindCurvePath[pts] :) $\endgroup$ Dec 14, 2019 at 3:18
2
$\begingroup$

Here is a solution based on a similar approach as shown by @b3m2a1: We cut out regions with a certain range of $z$ values by using PlotRange or RegionFunction:

Generate a example data set (taken from @b3m2a1's answer):

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];

Do the actual plotting for three different settings: (default, using PlotRange and using RegionFunction)

ListDensityPlot[surf, #] & /@ {
  {}, (* default behaviour, fills the hole *)
  PlotRange -> {0.07, All}, (* use PlotRange to limit the range of values *)
  RegionFunction -> (#3 > 0.07 &) (* use RegionFunction to limit the range *)
}

enter image description here

The first image shows how the plot looks like without any tweaking, so that a convex region is generated. The second example shows how PlotRange can be used to limit the range of $z$ values plotted. In more flexibility is needed (e.g. if you want to cut values in an interval), RegionFunction can be used: The function is called with $x,y,z$, where $x,y$ are the coordinates of the point and $z$ is the value at that point. To cut out values between 0.2 and 0.5 for example, we'd use

ListDensityPlot[surf, RegionFunction -> (! 0.2 <= #3 <= 0.5 &)]

enter image description here

$\endgroup$
1
  • $\begingroup$ In my opinion this is the correct answer. It’s simple and will always work for this type of problem. Mathematica does a constant fill when it fills in concavities so there should almost always be a way to describe that with RegionFunction. $\endgroup$
    – b3m2a1
    Dec 13, 2019 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.