Implementing Topography

[anneschulz1991]

Hello to all, I need help urgently. I have a very simple model in which I want to implement a topography. I had already managed to do this, but suddenly problems occur. I have now greatly simplified it and chosen a random topography. The command is also carried out without problems, however, in my output file for Pflotran all cells are included and also in the output of the inactive cells, the array is zero. This tells me that the topogrphy was not taken into account. Can someone explain this to me and possibly point out my errors? At the moment I am in despair and do not know how to help myself. Under the following link you can see my used files. https://upload.uni-jena.de/data/639aecdf5b66e9.53222843/Modell.zip Thanks for your help Anne

[Japhiolite]

Hi Anne,

I'll try to find time soon to have a look at it. Mind though, that the topography in GemPy is a masking array. So if you export, e.g. the lith_block, the cells above the topography will also be there.

Best Jan

[anneschulz1991]

Hello Jan, that would be fantastic! I have your comment in mind. But I thought inactive cells (i.e. cells that are above the topography) should be negatively assigned in the lith_block and unfortunately this is not the case for me, despite implemented topography. I have already checked this on an example (1.3c Adding topography to geological models) and there it works. Best Anne

[Japhiolite]

Hi Anne,

sorry for the late reply. I had a look into your code and found a peculiar behaviour (likely a bug). So what we expect from including topography is some kind of masked lith_block array in the solution.

In your example, the original solution yields a lith_block of:

In [1]: geo_model.solutions.lith_block
Out [1]: Lithology ids 
  [10.         10.         10.         ...  8.          8.
  7.95637271] 

the rounding to integers aside, the negative values for topography, which you mentioned are not visible there. However, once you call

In [2]: gp.plot_3d(geo_model,plotter_type="background",show_surfaces=True,show_lith=True,show_data=False,show_topography=True)

the lith_block becomes:

In [3]: geo_model.solutions.lith_block
Out [3]: array([  10.,   10.,   10., ..., -100., -100., -100.])

So now, the mask from topography is applied to the solutions (as the plotting apparently applies the mask to the lith_block array).

You can assign any value to the masked cells, it does not have to be -100. I'm not familiar with how Pflotran handles "inactive" cells, but maybe this method to give the topography-masked cells any value helps:

def topomask(geo_model, lith_block):
    """
    Method to mask the GemPy model with topography and change masked IDs to a new one for air.
    """
    model_shape = geo_model._grid.regular_grid.mask_topo.shape
    topo_mask = geo_model._grid.regular_grid.mask_topo
    
    # round lithblock and change to integers
    ids = np.round(lith_block)
    ids = ids.astype(int)
    
    mask_lith = ids.reshape(model_shape)
    mask_lith[topo_mask] = np.unique(mask_lith)[-1] + 1 # <-- here I give it the value max. Lith ID + 1, but can be ANY value
    
    return mask_lith

to get this into the solution you could simply assign it

lith_topo = topomask(geo_model, geo_model.solutions.lith_block)
lith_topo_flattened = lith_topo.flatten()
geo_model.solutions.lith_block = lith_topo_flattened

import gempy.utils.export as export
export.export_pflotran_input(geo_model, path='', filename="test.ugi")

Hope this helps.

Best and happy new year, Jan

[AlexanderJuestel]

@anneschulz1991, is this still an open issue?