Error in geo_model.add_surface_points()

[TobiasGlaubach]

Describe the bug Hi, when I try to generate a model from scratch, as shown in notebook: https://github.com/cgre-aachen/gempy/blob/master/notebooks/tutorials/ch1-Fundamentals/ch1-6_CreatingModels.ipynb I get an error in the add_surface_points method.

The error arises in the following code passage in gempy/core/data:

        try:
            self.df.loc[idx, 'surface'] = surface
        # ToDO test this
        except ValueError as error:
            self.del_surface_points(idx)
            print('The surface passed does not exist in the pandas categories. This may imply that'
                  'does not exist in the surface object either.')
            raise ValueError(error)

I am confused, since in your tutorial notebook it seems to be working with only slightly different values for the input data.

To Reproduce

import gempy as gp


extend = [0.0, 1.0, 0.0, 1.0, 0.0, 1.1]
discretization = [5,20,20]

x, y, z, f = 0.0, 0.0, 0.5, 'surface_2'

#%%
geo_model = gp.create_model('test')
gp.init_data(geo_model, extend, discretization)

geo_model.set_default_surfaces()
geo_model.set_default_orientation()

strats = ['surface_1', 'surface_2', 'basement']

gp.map_series_to_surfaces(geo_model, {'Strat_Series':  strats})
                  
geo_model.add_surface_points(x, y, z, f)
geo_model.add_orientations(x, y, z, f, pole_vector=N)  

Expected behavior should add a point to the geo model at 'surface_2'

Error Message:

The surface passed does not exist in the pandas categories. This may imply thatdoes not exist in the surface object either.
Traceback (most recent call last):

  File "C:\Users\tobia\Anaconda3\lib\site-packages\gempy\core\data.py", line 1396, in add_surface_points
    self.df.loc[idx, 'surface'] = surface

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 670, in __setitem__
    self._setitem_with_indexer(indexer, value)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 1023, in _setitem_with_indexer
    setter(item, v)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 961, in setter
    s._data = s._data.setitem(indexer=pi, value=v)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\internals\managers.py", line 561, in setitem
    return self.apply("setitem", **kwargs)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\internals\managers.py", line 442, in apply
    applied = getattr(b, f)(**kwargs)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\internals\blocks.py", line 1797, in setitem
    self.values[indexer] = value

  File "C:\Users\tobia\Anaconda3\lib\site-packages\pandas\core\arrays\categorical.py", line 2047, in __setitem__
    "Cannot setitem on a Categorical with a new "

ValueError: Cannot setitem on a Categorical with a new category, set the categories first


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File "C:\Users\tobia\OneDrive\Desktop\Projects\GeologicalModelling\GeoModellingJupyter\M2_GeoModelling\untitled0.py", line 31, in <module>
    geo_model.add_surface_points(x, y, z, f)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\gempy\plot\decorators.py", line 9, in pasp
    surface_points, idx = func(*args, **kwargs)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\gempy\core\model.py", line 560, in add_surface_points
    self.surface_points.add_surface_points(X, Y, Z, surface, idx)

  File "C:\Users\tobia\Anaconda3\lib\site-packages\gempy\core\data.py", line 1402, in add_surface_points
    raise ValueError(error)

ValueError: Cannot setitem on a Categorical with a new category, set the categories first

Desktop (please complete the following information):

• OS: Windows10
• Anaconda: Version 2020.02 Build py37_0
• gempy==2.0.1
• Theano==1.0.4

[Leguark]

Thank you for the detailed issue. You were right that it seemed quite confusing and probably we need much better error messages here. The problem is with the name of the surfaces, when you call set_default_surfaces you create Surface1 and Surface2 and then you try to map and add a point to Surface_2 so it does not exist.

Do you want to look where it fits and add a couple of Attribute Errors? Thanks

[TobiasGlaubach]

Hi, thanks for clearing this up. I would maybe check in the beginning of add_surface_points, if all surfaces which are being tried to add actually exist in the models surfaces, and if not throw an exception with proper info about the problem. Downside would be that this would introduce computational overhead.

Agree? If so I will implement this and make a Pull request.

[Leguark]

Sorry for the veeery late response I missed this issue completely. Yes a good error message in case the surface is not present is always a good idea. GemPy needs many of those specially on the public api