# Plot geologic profile of layers

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: 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,
ninterf}];
densities = Flatten[Import["den.dat", "Table"]];

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

(*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[], rangex[], step}, {j,
rangez[], rangez[], -step}];

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

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

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


Here is the final picture-density profile 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? The data for this complex model can be found here.

• 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. – Carl Lange Jan 17 '19 at 14:17
• @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). – Moonwalk Jan 17 '19 at 14:36

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

ListLinePlot[
interfaces,
PlotStyle -> ({CapForm["Butt"], #} & /@ Map[ColorData["Rainbow"], Rescale[densities]]),
Filling -> Bottom,
FillingStyle -> Opacity
] Filling will extrude your curve in the vertical direction, and so your other dataset comes out more jagged with this approach: 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:

Show[
MeshRegion[lines, PlotTheme -> "Lines"],
Graphics[{Red, MeshPrimitives[lines, {0, bds}]}],
AspectRatio -> 1/GoldenRatio
] 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"},
MeshCoordinates[lines][[bds]],
{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]]][];
envelope = MeshRegion[MeshCoordinates[lines][[bds]], Line[boundary]];
Show[
MeshRegion[lines, PlotTheme -> "Lines"],
MeshRegion[envelope, MeshCellStyle -> {1 -> Red}, PlotTheme -> "Lines"],
AspectRatio -> 1/GoldenRatio
] Let's join these up and find faces using IGraphM:

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


Now color each face:

colors = Rescale[Sort[densities]];
coords2 = MeshCoordinates[tofill];

Graphics[
GraphicsComplex[
coords2,
Map[
{
EdgeForm[{Thin, GrayLevel[.2]}],
ColorData["Rainbow"][Median[Flatten[cnf[coords2[[#]]]]]],
Polygon[#]
}&,
faces
]
],
AspectRatio -> 1/GoldenRatio
] • 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. – Moonwalk Jan 21 '19 at 8:04
• @Moonwalk it’s possible older versions couldn’t discretize the scene. What version are you using? – Chip Hurst Jan 21 '19 at 14:22
• I am using 11.2.0.0. I have not downloaded my update to 11.3 yet. – Moonwalk Jan 21 '19 at 15:33
• @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. – Chip Hurst Jan 21 '19 at 23:41
• Really nice answer, very graceful! – Carl Lange Jan 21 '19 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[], xrange[]},
Axes -> False, ImageSize -> Large,
AspectRatio -> 1,
PlotStyle -> (Append[Directive[CapForm["Butt"], Thick, #] & /@
(ColorData["Rainbow"] /@Rescale[densities]), None]),
Filling -> Thread[Range -> List /@ Range[2, 16]],
FillingStyle -> Opacity]
` 