Plot geologic profile of layers

Question

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[#[[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] &)]]]}, 
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.

Answer 1

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

Answer 2

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[1]
]

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]]][[2]];
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]];
cnf = Nearest[Catenate[MapThread[Thread[#1 -> #2] &, {interfaces, colors}]]];
coords2 = MeshCoordinates[tofill];

Graphics[
  GraphicsComplex[
    coords2, 
    Map[
      {
        EdgeForm[{Thin, GrayLevel[.2]}], 
        ColorData["Rainbow"][Median[Flatten[cnf[coords2[[#]]]]]], 
        Polygon[#]
      }&, 
      faces
    ]
  ], 
  AspectRatio -> 1/GoldenRatio
]