Displaying Fence Diagrams in PyVista

[AlexanderJuestel]

Description

Hello,

I would like to plot a geological cross section as a fence diagram (see image below). I got to the point of defining a plane and adding the raster as texture. However, it seems like I have not fully understood how a plane is defined.

This is what the profile looks like in map view. We have cropped the profile so that the beginning and the end of the profile should be mapped to the beginning and the end of the plane. The profile (or a section of it) is 7878 m long. The profile is about 2000 m hight. The angle between the horizontal and the profile trace on a map is 42.5 degrees striking from SW to NE.

We have defined the plane the following way. The center coordinates are the mid point of the raster, the direction was calculated using the formulas here: https://en.wikipedia.org/wiki/Cylindrical_coordinate_system#:~:text=%20%20%20%20v%20t%20e%20Orthogonal,%20%20Cartesian%20Cylindrical%20Spherical%20Paraboli%20...%20. I tried playing around with the i and j size but couldn't get to a proper solution.

import pyvista as pv
plane = pv.Plane(center=[32367511.857,5787584.251,1000], direction=[5808,5322,1000], i_size=0.5)

# Loading Raster
profile = rasterio.open('Profile.tif')

# Creating texture from array
tex = pv.numpy_to_texture(profile.read(1))

# Plotting plane with texture
plane.plot(texture=tex)

This is the current result:

Cheers Alex

[endarthur]

The docstring for pyvista.Plane look strange, as it asks for the Direction cylinder points to in [x,y,z] for the direction parameter. Maybe it could be normal vector, or something similar?

Using plane = pv.Plane(center= [32367511.857,5787584.251,-1000], direction=[np.sin(np.radians(135)), np.cos(np.radians(135)), 0], i_size=-2000, j_size=7878) seems to work, in this case.

[AlexanderJuestel]

That works perfectly @endarthur. Thanks for the help! (We only plot a section of it since the entire profile is not completely straight)

[banesullivan]

@AlexanderJuestel, let me know if you have an example where the section is not straight. This is very common and something we can handle just as easily

[AlexanderJuestel]

Do you mean curved or with multiple segments? We usually have multiple segments which would be awesome to plot right away :)

[banesullivan]

Either case will work!

[banesullivan]

If you have an example, add it here, and I'll reopen and try to demo it

[banesullivan]

Basically, if you have those well locations, then I can help you generate a segmented surface between them on which we can map the texture

[AlexanderJuestel]

We would have the start and end points of each segmentand the depth of each section. Would that be enough?

[AlexanderJuestel]

This section goes from +200m to -2000m in vertical extent Coordinates of the start point: 32362837 5769796 Coordinates of the bend: 32368424 5765456 (the bending point is indicated with a "K" at around half the section) Coordinates of the end point: 32374114 5763507

[banesullivan]

Here you go! Hopefully, you can follow the logic of how I generate the mesh and the texture coordinates and for more complicated cross-sections (N-many bends), you could implement an automated way to do this quite easily.

import pyvista as pv
import numpy as np

# Parameters for cross section
zrng = [-2000, 200]
start = [32362837, 5769796]
bend = [32368424, 5765456] # (the bending point is indicated with a "K" at around half the section)
end = [32374114, 5763507]

# Make a surface mesh representing that coverage
# - this mesh consists of 6 points. Generate them
#   a-----b-----e
#   |     |     |
#   |     |     |
#   d-----c-----f
a = start + [zrng[1],]
b = bend + [zrng[1],]
c = bend + [zrng[0],]
d = start + [zrng[0],]
e = end + [zrng[1],]
f = end + [zrng[0],]
# Now make a poly data mesh of these points
points = np.array([a, b, c, d, e, f])
faces = np.array([4, 0, 1, 2, 3, 4, 1, 4, 5, 2])
surface = pv.PolyData(points, faces)

# Map the texture to the mesh.
# - We know the tcoords of a, d, e, & f, but not necessarily b & c
# - to find them, scale by cell sizes:
# - Get the width of the two cells to find those coords
w = surface.compute_cell_sizes()["Area"] / np.ptp(zrng)
tw = (w / np.sum(w))[0]

# Generate Tcoords now!
t_coords = np.array([
    [0,1],
    [tw, 1],
    [tw, 0],
    [0, 0],
    [1, 1],
    [1, 0],
])
surface.t_coords = t_coords

# Load the texture image
texture = pv.Texture('./cross-section.png')

# Plot it up!
surface.plot(texture=texture, notebook=0, show_edges=True)

And look at the "K" in the correct spot 🚀

[banesullivan]

If you have a cross-section that has many more bends, post it here and I'll try to convert that snippet to a little function that takes a series of ordered control points and does all this automatically.

It'd be great to make that into an example in the docs

[AlexanderJuestel]

That is absolutely amazing! We have some more of these bigger (deep) sections that we could combine and that would intersect. They are part of a sedimentary basin in Northern Germany (the green part on the section). The data is publically available so we could create a nice example for that for the PyVista gallery :)

[akaszynski]

Keeping this open to remind us to add this as an example.

[AlexanderJuestel]

Thanks @akaszynski!

Here the example as promised, @banesullivan. The bends are indicated again with a "K" to keep it consistent. The vertical extent is +500 m to -6000 m. The coordinates are in order, so first coordinate is the point furthest to the left, last point is located furthest to the right. I also added a map view with the trace on the map for reference. The cross section goes from north to south in that case.

x values:

 [32383838.620400865,
 32379972.412131574,
 32378712.800021242,
 32376819.741358314,
 32375953.30297028,
 32378224.97336579,
 32381747.154628053,
 32383314.3887711,
 32385450.8929167,
 32386456.138052084]

y values:

[5806756.304725191,
 5804288.047468527,
 5798900.111274041,
 5797487.598271703,
 5786926.5152310245,
 5778557.011642429,
 5772623.000833637,
 5768227.100188477,
 5764516.343753345,
 5762795.47802485]

[akaszynski]

Vielen Dank!

[endarthur]

Hi! Sorry for barging in again, but I adapted Bane's code with some logic I used in another problem I had

import pyvista as pv
import numpy as np

X = np.array([
     [32383838.620400865,
 32379972.412131574,
 32378712.800021242,
 32376819.741358314,
 32375953.30297028,
 32378224.97336579,
 32381747.154628053,
 32383314.3887711,
 32385450.8929167,
 32386456.138052084],
    [5806756.304725191,
 5804288.047468527,
 5798900.111274041,
 5797487.598271703,
 5786926.5152310245,
 5778557.011642429,
 5772623.000833637,
 5768227.100188477,
 5764516.343753345,
 5762795.47802485]
]).T

n = len(X)
z0, z1 = -6000, 500

#duplicating the line, once with z=lower and another with z=upper values
vertices = np.zeros((2*n, 3))
vertices[:n, :2] = X
vertices[:n, 2] = z0
vertices[n:, :2] = X
vertices[n:, 2] = z1

# i+n --- i+n+1
# |\      |
# | \     |
# |  \    |
# |   \   |
# i  --- i+1
#
faces = np.array([[3, i, i + 1, i + n] for i in range(n-1)] + [[3, i+n+1, i+n, i+1] for i in range(n-1)])

# L should be the normalized to 1 cumulative sum of the segment lengths
L = np.linalg.norm(X[1:] - X[:-1], axis=1).cumsum()
L /= L[-1]
UV = np.zeros((2*n, 2))
UV[1:n, 0] = L
UV[n+1:, 0] = L
UV[:, 1] = np.repeat([0, 1], n)

surface = pv.PolyData(vertices, faces)

# Generate Tcoords now!
surface.t_coords = UV

# Load the texture image
texture = pv.Texture('./Gronau_GK_100_PyVista_Test.png')

# Plot it up!
surface.plot(texture=texture, notebook=0, show_edges=True)

This is the result:

[AlexanderJuestel]

That is what we will use it for later. We are going to display the sections together with a geological mode/a surface and borehole data and compare how they differ :)

I will add more cross sections tomorrow. It works perfectly!

[RichardScottOZ]

Nice! That was fast.

Probably worth a presentation this trick, not just an example writeup. Very impressive.

[banesullivan]

They are part of a sedimentary basin in Northern Germany (the green part on the section). The data is publically available so we could create a nice example for that for the PyVista gallery :)

@AlexanderJuestel, do you have a link to where these data are publically hosted?

[AlexanderJuestel]

Hey @banesullivan,

the data is available here (in German): https://www.opengeodata.nrw.de/produkte/geologie/geologie/GK/

We used either the 1:50000 or 1:100000 maps. What we did is we cropped the maps with the profiles and turned them a little bit around so that they are straight. We then made a final cropping to cut everything out that is not needed.

Let me know if you need anything else. Since this is working now, we may also provide you with some already cropped profiles until next week :)

[AlexanderJuestel]

Hey @banesullivan,

what information do you need to create an example or do you want me to create one?

Cheers Alex

[banesullivan]

that link doesn't load for me. Basically we just need the data files (or direct download links) and then the coordinates for making the various sections of them.

It would also be ideal to include a DEM and overlay a map like the one in https://github.com/pyvista/pyvista/issues/6196

[AlexanderJuestel]

Hey @banesullivan,

here is a link to a ZIP folder where all the data is in there. It includes an uncovered geological map, the cross sections, a shape file containing the coordinates of the bending points and a DEM (which needs to be cropped). Let me know if you need anything else.

https://drive.google.com/file/d/1dRZxOyeiXdGviNzSsohkC46_rlIGwsrP/view?usp=sharing

[banesullivan]

😍 That notebook is awesome!

@AlexanderJuestel, is it okay for me to pretty much turn that notebook you included into an example in the gallery?

[AlexanderJuestel]

Yes, go ahead :)

[RichardScottOZ]

@banesullivan not sure I have seen this one, is there a link? Thanks!

[banesullivan]

I hadn't gotten around to making a gallery example of this quite yet. Once I do, I'll post here and close the issue

[RichardScottOZ]

And I guess other than python, what can we export to for others to use these?

[RichardScottOZ]

Hi @endarthur in your bendy example above, what does the Gronau text png look like for that to work?

[RichardScottOZ]

e.g. I am just trying something out to see what happens (which isn't currently correct)

with here:-

and a texture

just end up with black, however

[RichardScottOZ]

e.g. this is a data version

but the texture created variety could be a small size approximation to locate things

[FCalderon3006]

That is what we will use it for later. We are going to display the sections together with a geological mode/a surface and borehole data and compare how they differ :)

I will add more cross sections tomorrow. It works perfectly!

Hi! Can you teach us how to show borehole with information like Lithology or % of ore grade ? Thank you in advance

[FCalderon3006]

Hi, I have a problem, I tried activate the path picking but the code only recognize the boundary of the section for picking, but i want to pick within the picture texture.

X = np.array(datos_organizados).T

n = len(X) z0, z1 = 3685, 4108

duplicating the line, once with z=lower and another with z=upper values

vertices = np.zeros((2*n, 3)) vertices[:n, :2] = X vertices[:n, 2] = z0 vertices[n:, :2] = X vertices[n:, 2] = z1

i+n --- i+n+1

|\ |

| \ |

| \ |

| \ |

i --- i+1

# faces = np.array([[3, i, i + 1, i + n] for i in range(n-1)] + [[3, i+n+1, i+n, i+1] for i in range(n-1)])

L = np.linalg.norm(X[1:] - X[:-1], axis=1).cumsum() L /= L[-1] UV = np.zeros((2*n, 2)) UV[1:n, 0] = L UV[n+1:, 0] = L UV[:, 1] = np.repeat([0, 1], n)

surface = pv.PolyData(vertices, faces) surface.t_coords = UV surface.active_t_coords = UV texture = pv.Texture('D:/2024/Modelo de Fundaciones/Otros/NE1.png')

Plot it up!

surface.plot(texture=texture, notebook=0, show_edges=False)

p = pv.Plotter() p.add_mesh(surface,texture = texture)

p.enable_path_picking(show_path=True)

p.enable_path_picking(pickable_window=True, show_path = True) p.show()

g = p.picked_path

p = pv.Plotter() p.add_mesh(surface, texture = texture) p.add_mesh(g) p.show()

horizonte = os.path.join(file_path, seccion.replace(".dxf","_horizonte.stl"))

g.save(horizonte)