gg.utils.show_number_of_data_points(geo_model=geo_model) adds columns with make problems with GemPy

[flohorovicic]

Describe the bug Executing the command gg.utils.show_number_of_data_points(geo_model=geo_model) adds columns to GemPy dataframe, which then leads to problems in GemPy plotting.

To Reproduce Steps to reproduce the behavior:

1. Execute gg.utils.show_number_of_data_points(geo_model=geo_model)
2. Compute model, etc.
3. Try 2D plot of top:

gp.plot_2d(geo_model, show_topography=True, section_names=['topography'], show_lith=True,
           show_boundaries=False, figsize=(12,10),
           kwargs_topography={'cmap': 'gray', 'norm': None, 'azdeg': 235})

Leads to error: ValueError: cannot reshape array of size 82944 into shape (216,128)

=> GemPy array seems to be three times the expected size. Removing the line above gg.utils... with the same code otherwise works.

[Japhiolite]

Indeed, it's the addition of columns to the surfaces dataframe of a gempy model, e.g. from your Uttendorf example we looked at: .

I ran through it with the Tutorial Chapter 1 Model, where the original dataframe is:

After computing the solution, the geological map provided is:

array([array([[3., 4., 4., ..., 2., 2., 2.]]),
       array([[-0.15531003, -0.15973748, -0.16175674, ..., -1.15124612,
               -1.15291493, -1.15482889],
              [ 0.73627075,  0.72629464,  0.7215417 , ...,  0.82527073,
                0.82295558,  0.8198365 ]])                             ],
      dtype=object)

Important is the first array withing this array, so 3., 4., 4., etc. These are the IDs of the units provided in the geo_model.surfaces dataframe. Now using gg.utils.show_number_of_data_points(geo_model=geo_model) yields:

Rerunning the model, the geological map of the solution is:

array([array([[ 3.,  4.,  4., ...,  2.,  2.,  2.],
              [13., 16., 16., ...,  8.,  8.,  8.],
              [ 0.,  1.,  1., ...,  1.,  1.,  1.]]),
       array([[-0.15531003, -0.15973748, -0.16175674, ..., -1.15124612,
               -1.15291493, -1.15482889],
              [ 0.73627075,  0.72629464,  0.7215417 , ...,  0.82527073,
                0.82295558,  0.8198365 ]])                             ],
      dtype=object)

So now, there are three "maps", where gempy sees the entries as surface IDs. For instance, the Siltstone with ID 3 has 13 Interface points and 0 orientations, which are the first entries in the two extra "maps" 13., and 0.,. ID 4 has 16 interface points and 1 orientation as visible from the dataframe, and so on.

Possible Solutions:

• Only work on a copy of the dataframe in GemGIS, and not on the original one. Is the information of the Number or datapoints somehow relevant to gempy? Otherwise, it should maybe only use a copy.
• Specify which column should be used in GemPy for creating (geological) maps. Currently, there's never been the thought of adding more columns to geo_model.surfaces.

In principle, I'm more in favor of the first solution, as a software should not change some vital object of another software its using...on the other hand, the construction of the geological map seems to need some work, if it can be broken that easily by simply adding columns with geo_model.add_surface_values()...

[Japhiolite]

On further inspection, I think that the method can be removed alltogether? Information about the surfaces, i.e. the number of data points per surface, can be called via geo_model.additional_data:

We could add just a line len surfaces orientation_points to also get the number of orientations per surface/unit.

I think that's the cleanest way of solving this issue and gemgis can still use the information, but simply gets them from the additional_data wrapper.

[AlexanderJuestel]

👍 For your thoughts@Japhiolite :)

[AlexanderJuestel]

Instead of appending to the existing GemPy DataFrame, a new one is created and returned for now. This should be fixed with #229