Saving a mesh (depth map) as tiff

[lucianomarquetto]

Hi all!

First I would like to thank the developers for this package, it helped me a lot in my research already.

Still, I do not have much experience in python, and thus I am having a hard time doing things that seem to be simple. Right now I want to transform a mesh created in Gemgis (depth map) to a tiff file I can load in QGIS.

I have computed my geological model in Gempy, then create depth maps of my units using:

#Creating` depth maps of units using Gemgis (based on the tutorial in 
#https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/18_creating_depth_maps_from_gempy_models.html).
dict_unit1 = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model,
                                                           surfaces='unit1')

mesh1 = dict_unit1['unit1'][0]
mesh1.save('mesh1.vtk')

mesh1

running 'mesh1' gives me the following result:

And then running print(type(mesh1)) shows it to be a <class 'pyvista.core.pointset.PolyData'>.

The problem is, I can't find a way of transforming this PolyData to an np.ndarray, which is needed to run 'raster.save_as_tiff' (as in https://gemgis.readthedocs.io/en/latest/getting_started/reference/raster_api/gemgis.raster.save_as_tiff.html?highlight=save%20as%20tif.

I tried using

mesh1_np = np.array(mesh1.points)

But this gives me an array with wrong max and min elevations, as it computes x,y and z, while it seems by other tutorials that only the z values are needed.

Then I tried something that is not in the tutorials, to extract only the z values from my array:

mesh1_np = np.array(mesh1.points) z_values = mesh1_np[:, 2] z_array = np.reshape(z_values, (mesh1.GetNumberOfPoints(), 1))

Running print(z_array.shape) gives me (4001, 1).

Then running:

gg.raster.save_as_tiff(raster=z_array, path=data_path + '/mesh1.tif', extent=[5.713e+05, 6.029e+05, 6.659e+06, 6.773e+06], crs='EPSG:31982', overwrite_file=True)

Generates a raster file but when I open this file in QGIS it shows like this:

I searched in other tutorials that had the "raster.save_as_tiff" expression, but could not figure out how to transform the PolyData to an np.ndarray. There are a lot of entries for this but I feel I am missing something, so decided to come here and ask if anyone would be willing to explain!

Thanks in advance!

Luciano.

[AlexanderJuestel]

Hey @lucianomarquetto,

I think I just have the thing you need :)

Have a look at https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/54_converting_pyvista_mesh_to_zmap_grid.html.

Use the function gg.utils.extract_zmap_data on your PyVista Meshes. The cell size will be the cell size of your raster later.

Avoid using the gg.utils.create_zmap_grid function. It is buggy once in a while.

Instead, use gg.raster.save_as_tiff on the z_values variable and provide all the information that are needed. Hope that helps!

[lucianomarquetto]

Hi @AlexanderJuestel - thanks to your help I am getting closer to what I need!

I followed your instructions and they worked. At first the raster was coming flipped (see images below), but then I switched the ymin and ymax in 'extent' and this solved the issue:

Instead of

gg.raster.save_as_tiff(raster=z_values, path=data_path + '/mesh1.tif', extent=[5.713e+05, 6.030e+05, 6.659e+06, 6.773e+06], crs='EPSG:31982', overwrite_file=True)

I fliped the y extents as:

extent=[5.713e+05, 6.030e+05, 6.773e+06, 6.659e+06],

Wrong raster orientation:

Right raster orientation:

But anyway, thanks for the help!