gempy
gempy copied to clipboard
KeyError when show_topography=True with layer not in lith_block (e.g., layer below resolution of blocks)
Describe the bug
Whenever a layer is too thin/substantially thinner than the resolution, its colour gets mapped but is not actually present in the lith_block data.
Plotting this in 2D/3D is fine as long as the plot has Topography disabled, even if topography is applied to the generated surfaces.
However, plotting in 3D with show_topography= True
results in the following KeyError:
Log dump/KeyError
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Input In [63], in <cell line: 1>()
----> 1 gp.plot_3d(geo_model, plotter_type="basic", off_screen=False,
2 show_topography=True,
3 show_scalar=False,
4 show_lith=True, show_data=False)
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\gempy\plot\plot_api.py:336, in plot_3d(model, plotter_type, show_data, show_results, show_surfaces, show_lith, show_scalar, show_boundaries, show_topography, scalar_field, ve, kwargs_plot_structured_grid, kwargs_plot_topography, kwargs_plot_data, image, off_screen, **kwargs)
333 gpv.plot_data(**kwargs_plot_data)
335 if show_topography and model._grid.topography is not None:
--> 336 gpv.plot_topography(**kwargs_plot_topography)
338 if ve is not None:
339 gpv.p.set_scale(zscale=ve)
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\gempy\plot\vista.py:557, in GemPyToVista.plot_topography(self, topography, scalars, contours, clear, **kwargs)
553 colors_rgb = pd.DataFrame(colors_rgb_.to_list(), index=colors_hex.index) * 255
555 sel = np.round(self.model.solutions.geological_map[0]).astype(int)[0]
--> 557 scalars_val = numpy_to_vtk(colors_rgb.loc[sel], array_type=3)
558 cm = mcolors.ListedColormap(list(self._get_color_lot(is_faults=True, is_basement=True)))
559 rgb = True
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\pandas\core\indexing.py:931, in _LocationIndexer.__getitem__(self, key)
928 axis = self.axis or 0
930 maybe_callable = com.apply_if_callable(key, self.obj)
--> 931 return self._getitem_axis(maybe_callable, axis=axis)
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\pandas\core\indexing.py:1153, in _LocIndexer._getitem_axis(self, key, axis)
1150 if hasattr(key, "ndim") and key.ndim > 1:
1151 raise ValueError("Cannot index with multidimensional key")
-> 1153 return self._getitem_iterable(key, axis=axis)
1155 # nested tuple slicing
1156 if is_nested_tuple(key, labels):
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\pandas\core\indexing.py:1093, in _LocIndexer._getitem_iterable(self, key, axis)
1090 self._validate_key(key, axis)
1092 # A collection of keys
-> 1093 keyarr, indexer = self._get_listlike_indexer(key, axis)
1094 return self.obj._reindex_with_indexers(
1095 {axis: [keyarr, indexer]}, copy=True, allow_dups=True
1096 )
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\pandas\core\indexing.py:1314, in _LocIndexer._get_listlike_indexer(self, key, axis)
1311 else:
1312 keyarr, indexer, new_indexer = ax._reindex_non_unique(keyarr)
-> 1314 self._validate_read_indexer(keyarr, indexer, axis)
1316 if needs_i8_conversion(ax.dtype) or isinstance(
1317 ax, (IntervalIndex, CategoricalIndex)
1318 ):
1319 # For CategoricalIndex take instead of reindex to preserve dtype.
1320 # For IntervalIndex this is to map integers to the Intervals they match to.
1321 keyarr = ax.take(indexer)
File ~\.conda\envs\gemgis_deltaneset\lib\site-packages\pandas\core\indexing.py:1377, in _LocIndexer._validate_read_indexer(self, key, indexer, axis)
1374 raise KeyError(f"None of [{key}] are in the [{axis_name}]")
1376 not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())
-> 1377 raise KeyError(f"{not_found} not in index")
KeyError: '[21] not in index'
Indeed, there is no lithology with key 21 in the lith_block:
21 in geo_model.solutions.lith_block
> False
To Reproduce Provide detailed steps to reproduce the behavior:
See above.
Expected behavior
- Prior to calling numpy_to_vtk, make sure to check all lith_block IDs are actually in the lith_block
- Only apply numpy_to_vkt (and colouring) with the lith_block IDs that are present.
- Plot the topography with lith_block colours dynamically -> If a higher resolution topography would should more of the lith_blocks, than these should become visible after rerunning the models.
Desktop (please complete the following information):
- OS: Windows 11
- GemPy Version: '2.2.11'