geost.models.model_dataset.ModelDataset.get_top_bottom#

ModelDataset.get_top_bottom(condition: GeotopUnits | DataArray) → Dataset[source]#

Derive the top and bottom depths at each x,y-location of a specified condition from a voxelmodel or layermodel. The condition can be a boolean DataArray or Dataset that indicates which voxels or layers to include in the thickness calculation.

Parameters:

condition (GeotopUnits | xr.DataArray) – A boolean DataArray or Dataset indicating which voxels or layers to include. The condition should have the same dimensions as the model.

Returns:

xarray.Dataset containing “top” and “bottom” DataArrays for each horizontal x,y-location in the model.

Return type:

xr.Dataset

Examples

Determine the top and bottom in a voxelmodel Dataset where “strat” equals 1100 and “lith” equals 1:

>>> top_bottom = voxelmodel.gst.get_top_bottom(
...    (voxelmodel["lithology"] == 1) & (voxelmodel["strat"] == 1100)
... )

Or in a layermodel for a subset of units where the a value is smaller than 10:

>>> top_bottom = layermodel.gst.get_top_bottom(
...    (layermodel["layer"].isin(["B", "D"])) & (layermodel["value"] < 10)
... )

If you are working with GeoTOP, you can also use a GeotopUnits object to specify the condition, for example to get the top and bottom of the “Formatie van Echteld” unit with clay as the lithology:

>>> geotop = geost.read_geotop_from_opendap(bbox=(110_000, 440_000, 120_000, 450_000))
>>> strat_units = geost.bro.geotop_strat_units()
>>> lithok_units = geost.bro.geotop_lithok_units()
>>> top_bottom = geotop.gst.get_top_bottom(
...     strat_units.select_description_contains("Formatie van Echteld")
...     & lithok_units.select_unit("k")
... )