I am trying to plot a set of layers with different values of densities (density profile). Geometry ({x,z} values) of 15 interfaces of layers look like this: enter image description here

while the first interface is a planar surface. {x, z} values and densities of individual interfaces can be found here.

I tried to use ArrayPlot. Here is my attempt:

(*read data*)
ninterf = 15;
interfaces = 
Table[Import["interface_" <> ToString[i] <> ".dat", "Table"], {i, 1,
densities = Flatten[Import["den.dat", "Table"]];

rangex = MinMax[#[[1]] & /@ interfaces[[1]]]
rangez = {Min[#[[2]] & /@ interfaces[[-1]]], 
Max[#[[2]] & /@ interfaces[[1]]]}

(*interpolate data*)
f = Table[Unique[f$], {i, 1, ninterf}];
Table[f[[i]] = Interpolation[interfaces[[i]], InterpolationOrder -> 1], {i, 
1, ninterf}];
F = Through@*f;

(*create a grid of {x,z} values for ArrayPlot*)
step = 1.;
xzvalues = 
Table[{N[i], N[j]}, {i, rangex[[1]], rangex[[2]], step}, {j, 
rangez[[2]], rangez[[1]], -step}];

(*determine density for every point of the grid*)
densityfunc[xval_, zval_] := 
Module[{numb = 
Last[Flatten[Position[F[N[xval]], _?(# >= N[zval] &)]]]}, 

values = Apply[densityfunc, #, {1}] & /@ xzvalues;

(*Show result-ArrayPlot*)
ArrayPlot[Transpose[values], ColorFunction -> "Rainbow"]

Here is the final picture-density profile

enter image description here

The time needed for calculation of values on my computer is 689.04 seconds (for bigger models, e.g. 5 times bigger, this might take a lot of time). In this case, I set the grid spacing step=1. since interfaces should look smooth. With bigger step the calculation time would be shorter but interfaces would not be smooth enough.

I would like to speed up the evaluation of values(if it is possible; parallelization of code would be also welcomed) or maybe someone knows faster way how to get similiar picture of density profile. I will be grateful for any suggestions.

UPDATE: Carl Lange and kglr proposed using Filling which works great for this simple model of 15 interfaces. Could be Filling used also for this complex model of layers? enter image description here

The data for this complex model can be found here.

  • 2
    $\begingroup$ I'm not quite clear on what's going on here - perhaps you can post a little more context and some of your data. It seems like what you're generating can be achieved just by using Filling, but I think I'm just misunderstanding your exact problem and intention. $\endgroup$
    – Carl Lange
    Jan 17, 2019 at 14:17
  • $\begingroup$ @CarlLange I forgot about Filling.Thank you. Yes, in this case it solves my problem. I just need to think about it little more. Maybe is not applicable for all configurations (more complex density profiles). $\endgroup$
    – Moonwalk
    Jan 17, 2019 at 14:36
  • 2
    $\begingroup$ see also StackedListPlot $\endgroup$
    – kglr
    Jan 17, 2019 at 20:55

2 Answers 2


As suggested in the comments we can use ListLinePlot with Filling like so:

  PlotStyle -> ({CapForm["Butt"], #} & /@ Map[ColorData["Rainbow"], Rescale[densities]]),
  Filling -> Bottom,
  FillingStyle -> Opacity[1]

enter image description here

Filling will extrude your curve in the vertical direction, and so your other dataset comes out more jagged with this approach: enter image description here

To fix this, I will clean the data, then isolate and color each face manually. Since these faces will have edges of different colors, I choose the median of the face vertex colors.

First I convert our scene into a mesh, then find all unconnected endpoints. We'll need to close these up somehow.

lines = DiscretizeGraphics[Line /@ interfaces];
coords = MeshCoordinates[lines];
bds = Flatten[Position[Differences[lines["ConnectivityMatrix"[0, 1]]["RowPointers"]], 1]];

Let's pause and look at these endpoints. It's clear there's some that really shouldn't be there:

  MeshRegion[lines, PlotTheme -> "Lines"], 
  Graphics[{Red, MeshPrimitives[lines, {0, bds}]}], 
  AspectRatio -> 1/GoldenRatio

enter image description here

We'll need to clean up the original data a bit. I notice the endpoints we'd like to clean are all within 5.2 of another vertex, and the endpoints we want to keep are not (except for one). I will group these nearby points and merge them by taking the mean:

neardata = Nearest[
  MeshCoordinates[lines] -> {"Index", "Element"}, 
  {All, 5.2}

(coords[[#1]] = ConstantArray[#2, Length[#1]]) & @@@ MapAt[Mean, Transpose /@ neardata, {All, 2}];
lines = MeshRegion[coords, MeshCells[lines, 1]];
bds = Flatten[Position[Differences[lines["ConnectivityMatrix"[0, 1]]["RowPointers"]], 1]];

We'll now connect up the endpoints. If the set was convex, we'd use ConvexHullMesh. Instead we could find an alpha-shape, but something like FindShortestTour seems to work just fine.

boundary = FindShortestTour[MeshCoordinates[lines][[bds]]][[2]];
envelope = MeshRegion[MeshCoordinates[lines][[bds]], Line[boundary]];
  MeshRegion[lines, PlotTheme -> "Lines"], 
  MeshRegion[envelope, MeshCellStyle -> {1 -> Red}, PlotTheme -> "Lines"], 
  AspectRatio -> 1/GoldenRatio

enter image description here

Let's join these up and find faces using IGraphM:

tofill = RegionUnion[lines, envelope];
faces = Rest[IGFaces[IGMeshGraph[tofill]]];

Now color each face:

colors = Rescale[Sort[densities]];
cnf = Nearest[Catenate[MapThread[Thread[#1 -> #2] &, {interfaces, colors}]]];
coords2 = MeshCoordinates[tofill];

        EdgeForm[{Thin, GrayLevel[.2]}], 
  AspectRatio -> 1/GoldenRatio

enter image description here

  • $\begingroup$ Thank you for taking the time to write an elaborate answer for me. I think that this answer fuly accomplished what I need. Right now, I am trying to reproduce your code, but this first line of code lines=DiscretizeGraphics[Line[interfaces]]; gives me this error: DiscretizeGraphics: Graphics primitives with coordinates are not supported. $\endgroup$
    – Moonwalk
    Jan 21, 2019 at 8:04
  • $\begingroup$ @Moonwalk it’s possible older versions couldn’t discretize the scene. What version are you using? $\endgroup$
    – Greg Hurst
    Jan 21, 2019 at 14:22
  • $\begingroup$ I am using I have not downloaded my update to 11.3 yet. $\endgroup$
    – Moonwalk
    Jan 21, 2019 at 15:33
  • $\begingroup$ @Moonwalk Looks like lines = DiscretizeGraphics[Line /@ interfaces]; works in 11.2. I went ahead and made this change in my post. The rest worked for me in 11.2. $\endgroup$
    – Greg Hurst
    Jan 21, 2019 at 23:41
  • 1
    $\begingroup$ Really nice answer, very graceful! $\endgroup$
    – Carl Lange
    Jan 21, 2019 at 23:43
ifs = Interpolation[#, InterpolationOrder -> 1] & /@ interfaces;
xrange = MinMax[interfaces[[All, All, 1]]]; 
Plot[Evaluate[Append[Through@ifs@x,  1.05 Min[interfaces[[All, All, -1]]]]],
 {x, xrange[[1]], xrange[[2]]}, 
 Axes -> False, ImageSize -> Large, 
 AspectRatio -> 1, 
 PlotStyle -> (Append[Directive[CapForm["Butt"], Thick, #] & /@ 
    (ColorData["Rainbow"] /@Rescale[densities]), None]), 
 Filling -> Thread[Range[15] -> List /@ Range[2, 16]], 
 FillingStyle -> Opacity[1]]

enter image description here


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.