Source code for ogstools.meshlib.boundary_set

# Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
#            Distributed under a Modified BSD License.
#            See accompanying file LICENSE.txt or
#            http://www.opengeosys.org/project/license
#

import tempfile
from abc import ABC, abstractmethod
from collections import namedtuple
from itertools import pairwise
from pathlib import Path

import pandas as pd

from ogstools.meshlib.boundary import Layer, LocationFrame, Raster
from ogstools.meshlib.boundary_subset import Surface




class BoundarySet(ABC):
    """
    Abstract base class representing a collection of boundaries with constraints.

    A BoundarySet is composed of multiple boundaries linked with constraints:
    - Free of gaps and overlaps.
    - Distinguished by markers to identify different boundaries.
    - Adherence to rules of `piecewise linear complex (PLC)`.
    """



    @abstractmethod
    def bounds(self) -> list:
        return []




    @abstractmethod
    def filenames(self) -> list[Path]:
        return []






class LayerSet(BoundarySet):
    """
    Collection of geological layers stacked to represent subsurface arrangements.

    In a geological information system, multiple layers can be stacked vertically to
    represent the subsurface arrangement. This class provides methods to manage and
    work with layered geological data.
    """



    def __init__(self, layers: list[Layer]):
        """
        Initializes a LayerSet. It checks if the list of provided layers are given in a top to bottom order.
        In neighboring layers, layers share the same surface (upper bottom == low top).
        """

        for upper, lower in pairwise(layers):
            if upper.bottom != lower.top:
                msg = "Layerset is not consistent."
                raise ValueError(msg)
        self.layers = layers




    def bounds(self) -> list:
        return list(self.layers[0].top.mesh.bounds)




    def filenames(self) -> list[Path]:
        layer_filenames = [layer.bottom.filename for layer in self.layers]
        layer_filenames.insert(0, self.layers[0].top.filename)  # file interface
        return layer_filenames




    @classmethod
    def from_pandas(cls, df: pd.DataFrame) -> "LayerSet":
        """Create a LayerSet from a Pandas DataFrame."""
        Row = namedtuple("Row", ["material_id", "mesh", "resolution"])
        surfaces = [
            Row(
                material_id=surface._asdict()["material_id"],
                mesh=surface._asdict()["filename"],
                resolution=surface._asdict()["resolution"],
            )
            for surface in df.itertuples(index=False)
        ]
        base_layer = [
            Layer(
                top=Surface(top.mesh, material_id=top.material_id),
                bottom=Surface(bottom.mesh, material_id=bottom.material_id),
                num_subdivisions=bottom.resolution,
            )
            for top, bottom in pairwise(surfaces)
        ]
        return cls(layers=base_layer)




    def create_raster(self, resolution: float) -> tuple[Path, Path]:
        """
        Create raster representations for the LayerSet.

        This method generates raster files at a specified resolution for each layer's
        top and bottom boundaries and returns paths to the raster files.
        """
        bounds = self.layers[0].top.mesh.bounds
        raster_set = self.create_rasters(resolution=resolution)

        locFrame = LocationFrame(
            xmin=bounds[0], xmax=bounds[1], ymin=bounds[2], ymax=bounds[3]
        )
        # Raster needs to be finer then asc files. Otherwise Mesh2Raster fails
        raster = Raster(locFrame, resolution=resolution * 0.95)
        raster_vtu = Path(tempfile.mkstemp(".vtu", "raster")[1])
        raster.as_vtu(raster_vtu)
        rastered_layers_txt = Path(
            tempfile.mkstemp(".txt", "rastered_layers")[1]
        )
        with rastered_layers_txt.open("w") as file:
            file.write("\n".join(str(item) for item in raster_set))
        return raster_vtu, rastered_layers_txt




    def create_rasters(self, resolution: float) -> list[Path]:
        """
        For each surface a (temporary) raster file with given resolution is created.
        """
        rasters = [self.layers[0].top.create_raster_file(resolution=resolution)]
        for layer in self.layers:
            r = layer.create_raster(resolution=resolution)
            rasters.extend(r[1:])

        return list(dict.fromkeys(rasters))