Volume calculation between models with different topography

[dmartinezdi]

What is your question?
Hi, I am trying to compare two different models that have the same data but different topography and I am trying to calculate the volume difference between them, but the values that I got from 'raw_arrays' are the same I think because it took the same data to build the model, is there anyway to get the values with the topography cut to be able to measure the volume?

[javoha]

Hi,
this seems to be a little hidden in gempy and as far as I know this is not an implemented function. But you can access the vertices of the topogrpahy grid using geo_model.grid.topography.values. This might allow you to mask the solution from raw_arrays to get the clipped result.

Probably more straightforward: geo_model.grid.topography.topography_mask returns the mask that can directly be applied to the block solution (e.g. the solution in shape of your resolution).
Let me know if this is what you were looking for or if I got something wrong.
Cheers
Jan

Small edit: I will convert this issue to a discussion for now as it does not strike me as an issue or bug.

[javoha]

Hi, glad that it worked. In general it should be young to old, so a value of lith_block = 1 is the top element (youngest), lith_block = 3 in the given example would be the basement. I am not entirely sure if faults would change this but I think they should not. Here also another quick and dirty way to check:

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

p= ax.scatter(model_1.grid.regular_grid.values[:,0], 
           model_1.grid.regular_grid.values[:,1], 
           model_1.grid.regular_grid.values[:,2], 
           c=model_1.solutions.raw_arrays.lith_block,
      )

fig.colorbar(p)

Cheers,
Jan

[dmartinezdi]

Thanks Jan for all your help with this!

I have one more question but it is not related to this, I am doing a model where I have some Quaternary deposits, thus they are not regular through all the area, do you know any way to model this easily? Right now I am trying to add some points that are above the topography to correct the surfaces but I am not getting everything right.

Cheers,
David

[javoha]

Hi
sorry I somehow missed your new question: Yes, this might be tricky to be honest. If you are able (allowed) to, feel free to send me your model/code and I could take a look. Adding points manually above topography is the best gempy option I can think of. Alternatively you could try adding these structures manually after computing the gempy model. But it depends very much on the exact shapes we are looking at, so I would need to take a look.
Maybe also making it an onlap relationship might help, but again - hard to say without a more detailed look.

[dmartinezdi]

Hi @javoha no worries, thanks for your reply.
Yeah sure I can send it to you it but I would need an email address to be able to send it, I am going to try the onlap relationship but It will be amazing if you could take a look on it.
Thanks

[dmartinezdi]

Hi @javoha Jan I have sent you the email I do not know if you have seen it and had the time to check it, I am just wondering if it get through or if it might have gone to the spam or something.
Once again thank you,
David

[javoha]

Hi, sorry for the silence - I got your model and took a quick look but I didn't have the time to play around with some options. I will let you now as soon as I have tried things.

[dmartinezdi]

Hi Jan, no worries, thank you for taking time to look at it, I was just checking if the mail have gone through.

[javoha]

Hi,
so I was able to take a look: First, regarding the holes that sometimes appear in the meshes - These are artefacts of the dual contouring. We are aware of them and are working on a solution. It does not change your model though, as the underlying scalar field is the same, it is a bug purely in the mesh extraction step for visualization.
Regarding the bigger picture : I actually cannot run the model with your huge topography tif. What machine are you working on and how long didi it run? From the precomputed solution you send me I can see the issue. Have you tried to calculate the model without topography in higher resolution to identify areas where it does not fit your expectation?
Cheers
Jan

[dmartinezdi]

Hi Jan,

Oh okay okay I got it, I was trying to extract the mesh surface to make some calculations individually but I have not been able to do it successfully as I was trying to extract them with gemgis but I got something without the correct xyz coordinates, do you happen to have any idea how I can extract them with them?

Oh god, I can send you a less detailed topography if you like, I am running the model in a Thinkpad e14 with a i7 1165g7 16 gb ram and it takes like 3 min and a half to run with the topography. well the thing to calculate it without the topograpy is that I am not being able to check how the Quaternary behaves as it is not cut by the topography itself. Regarding the volume thing that you helped me through I was wondering if you might know a way to calculate it only inside a polygon that I have in a shp file?

Cheers
David

[javoha]

Ah yes, so extracting the meshes works, but you need to back-transform the coordinates of the vertices to get the correct coordinates. Check #916 for example. That way you are able to get the meshes to use in other calculations or software.
Regarding the volume by shapefile - I guess you would need to create some kind of mask to extract only the values within. I am not sure what the smartest way to do that would be, but it should be rather straightfoward.

[dmartinezdi]

Hi Jan

Thanks for that I have been able to extract those and generate the surface, and finishing the volume calculation.

Also I have sent you two new raster files to see if those could work to see if you can check the model with those, one of them the pixel size is 5 m and the other one is 10 m. I have attached those to your email.

Thanks again

[dmartinezdi]

Hi @javoha Jan, i have been trying to make the onlap surface but I do not understand very well how to do it, and have not been successful in the matter, would you be able to guide me through that process please?

[javoha]

Hi,
yeah I will try. Again I am not 100% sure it solves your issues but let's see:

First we need to recall again how gempy interpolates. Based on the strucutral frame, elements (rock units) within a group can never crosscut as they are interpolated with a single scalar field. GEmpy uses a bottom to top approach, so younger units within a group will generally conformibly lie on top of older groups. When we have multiple structural groups, younger groups will unconformably cut of older groups. This is indicated by the default structural relation ERODE.
When we want to model an onlap erlationship we need a younger group to NOT erode the older groups - this is a behavior that needs a little trick in gempy.
I do think the example we have is really helpful: LINK
Looking at the cross section image:

The normal gempy behavior would lead to the structural group right_series containing rock1 and rock2 cutting of the older group left_series with rock3. In order to get right_series to onlap on left_series we need to change two things:

1. We need to set its strucutral relation to ONLAP, forcing it to not cut older groups
2. As gempy normally goes from bottom to top, but top (younger) now is forced to not cut older (bottom), we need to add the upper bound of our left_series. This forces us to add a new structural group with an element we called onlap_surface in this example. Note that itself it still has the ERODE stack relation as it limits the left_series. This element is what our younger group reight series will now onlap on.

Does this make sense? I really hope it does :)

Feel free to ask if anything is unclear.

Cheers,
Jan

[dmartinezdi]

Hi @javoha,

Thank you, that sounds good but I do not understand very well how to setup the onlap relation, so I have this structural frame:

and I was trying to put it between Kd3 and Tba but I do not get pretty well if I have to had the data as the onlap and then data as an erode version of it? like I would put Tba as the onlap and then again as erode or how should I do that?

Also I have another little question if you happen to know, I was trying to get out the model plots in 3d and 2d with a scale bar is that possible? I saw like some kind of model in gpyviewer but I am not certain how to use it. Is it possible? And, is it possible to extract like the 2d plots as other format as a shp or a csv with the data of the cell or something like that?

Thank you

Cheers,
David

[javoha]

Hi,
so the onlap relation mostly makes sense if you have multiple elements in a structural group that onlap on a surface. And yes - then you need to define it as an additional surface. If the onlapping group is just a single element, using erode as of usual kind of works the same way. Maybe this ugly sketch helps:

Case A is when you need an additional onlap surface (red), as multiple structural elements of the younger group (orange, green) onlap on this surface. In Case B it does not matter as the structure is just filled with a single element and you can just put it to erode.

Regarding your second question: What do you mean by scale bar. It is possible to plot the scalar field instead of the lithologies. And as you can extract the lithology block you basically have the data per cell as well as the surfaces meshes for 3D plotting wherever you want. Is that what you mean?

Cheers,
Jan