Add surface set model (#604)

CHANGELOG.md

@@ -6,6 +6,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [UNRELEASED] - YYYY-MM-DD
 ### Changed
+- [#604](https://github.com/equinor/webviz-subsurface/pull/604) - Consolidates surface loading and statistical calculation of surfaces by introducing a shared
+SurfaceSetModel. Refactored SurfaceViewerFMU to use SurfaceSetModel.
 - [#586](https://github.com/equinor/webviz-subsurface/pull/586) - Added phase ratio vs pressure and density vs pressure plots. Added unit and density functions to PVT library. Refactored code and added checklist for plots to be viewed in PVT plot plugin. Improved the layout.
 - [#599](https://github.com/equinor/webviz-subsurface/pull/599) - Fixed an issue in ParameterAnalysis where the plugin did not initialize without FIELD vectors 
 

mypy.ini

@@ -43,3 +43,6 @@ ignore_errors=True
 
 [mypy-tests.unit_tests.model_tests.test_ensemble_set_model]
 ignore_errors=True
+
+[mypy-tests.unit_tests.model_tests.test_surface_set_model]
+ignore_errors=True

setup.py

@@ -74,7 +74,7 @@
         "statsmodels>=0.12.1",  # indirect dependency through https://plotly.com/python/linear-fits/
         "webviz-config>=0.2.7",
         "webviz-subsurface-components>=0.3.0",
-        "xtgeo>=2.8",
+        "xtgeo>=2.14",
     ],
     extras_require={"tests": TESTS_REQUIRE},
     setup_requires=["setuptools_scm~=3.2"],

tests/integration_tests/test_surface_selector.py

@@ -28,7 +28,7 @@
 
 return_value = {
     "name": "lowerreek",
-    "attr": "oilthickness",
+    "attribute": "oilthickness",
     "date": "20030101_20010601",
 }
 

tests/unit_tests/model_tests/test_surface_set_model.py

@@ -0,0 +1,62 @@
+from pathlib import Path
+
+import pytest
+import xtgeo
+
+from webviz_subsurface._models.surface_set_model import SurfaceSetModel
+from webviz_subsurface._datainput.fmu_input import find_surfaces
+
+
+@pytest.mark.usefixtures("app")
+def test_surface_set_model(testdata_folder):
+    ensemble_paths = {
+        "iter-0": str(
+            Path(testdata_folder / "reek_history_match" / "realization-*" / "iter-0")
+        )
+    }
+
+    surface_table = find_surfaces(ensemble_paths)
+    surface_table = surface_table.drop("ENSEMBLE", axis=1)
+
+    smodel = SurfaceSetModel(surface_table)
+    assert set(smodel.attributes) == set(
+        [
+            "average_pressure",
+            "average_swat",
+            "average_permz",
+            "amplitude_max",
+            "average_soil",
+            "average_permx",
+            "average_sgas",
+            "amplitude_rms",
+            "amplitude_min",
+            "facies_fraction_channel",
+            "perm_average",
+            "oilthickness",
+            "poro_average",
+            "zonethickness_average",
+            "timeshift",
+            "ds_extracted_horizons",
+            "facies_fraction_crevasse",
+            "stoiip",
+            "average_poro",
+        ]
+    )
+    assert set(smodel.names_in_attribute("ds_extracted_horizons")) == set(
+        ["topupperreek", "baselowerreek", "toplowerreek", "topmidreek"]
+    )
+    real_surf = smodel.get_realization_surface(
+        attribute="ds_extracted_horizons", name="topupperreek", realization=0
+    )
+    assert isinstance(real_surf, xtgeo.RegularSurface)
+    assert real_surf.values.mean() == pytest.approx(1706.35, 0.00001)
+    stat_surf = smodel.calculate_statistical_surface(
+        attribute="ds_extracted_horizons", name="topupperreek"
+    )
+    assert isinstance(stat_surf, xtgeo.RegularSurface)
+    assert stat_surf.values.mean() == pytest.approx(1706.42, 0.00001)
+
+    stat_surf = smodel.calculate_statistical_surface(
+        attribute="ds_extracted_horizons", name="topupperreek", realizations=[2, 4, 5]
+    )
+    assert stat_surf.values.mean() == pytest.approx(1707.07, 0.00001)

webviz_subsurface/_datainput/fmu_input.py

@@ -1,7 +1,7 @@
 from pathlib import Path
 import warnings
 import glob
-from typing import Union, Optional
+from typing import Union, Optional, List
 
 import pandas as pd
 from fmu.ensemble import ScratchEnsemble, EnsembleSet
@@ -129,7 +129,11 @@ def find_sens_type(senscase: str) -> Optional[str]:
 @CACHE.memoize(timeout=CACHE.TIMEOUT)
 @webvizstore
 def find_surfaces(
-    ensemble_paths: dict, suffix: str = "*.gri", delimiter: str = "--"
+    ensemble_paths: dict,
+    surface_folder: str = "share/results/maps",
+    surface_files: Optional[List] = None,
+    suffix: str = "*.gri",
+    delimiter: str = "--",
 ) -> pd.DataFrame:
     """Reads surface file names stored in standard FMU format, and returns a dictionary
     on the following format:
@@ -143,24 +147,27 @@ def find_surfaces(
     """
     # Create list of all files in all realizations in all ensembles
     files = []
-    for ens, path in ensemble_paths.items():
-        path = Path(path)
+    for _, ensdf in get_realizations(ensemble_paths=ensemble_paths).groupby("ENSEMBLE"):
         ens_files = []
-        for realpath in glob.glob(str(path / "share" / "results" / "maps" / suffix)):
-            stem = Path(realpath).stem.split(delimiter)
-
-            if len(stem) >= 2:
-                ens_files.append(
-                    {
-                        "ENSEMBLE": ens,
-                        "path": realpath,
-                        "name": stem[0],
-                        "attribute": stem[1],
-                        "date": stem[2] if len(stem) >= 3 else None,
-                    }
-                )
+        for _real_no, realdf in ensdf.groupby("REAL"):
+            runpath = realdf.iloc[0]["RUNPATH"]
+            for realpath in glob.glob(str(Path(runpath) / surface_folder / suffix)):
+                filename = Path(realpath)
+                if surface_files and filename.name not in surface_files:
+                    continue
+                stem = filename.stem.split(delimiter)
+                if len(stem) >= 2:
+                    ens_files.append(
+                        {
+                            "path": realpath,
+                            "name": stem[0],
+                            "attribute": stem[1],
+                            "date": stem[2] if len(stem) >= 3 else None,
+                            **realdf.iloc[0],
+                        }
+                    )
         if not ens_files:
-            warnings.warn(f"No surfaces found for ensemble located at {path}.")
+            warnings.warn(f"No surfaces found for ensemble located at {runpath}.")
         else:
             files.extend(ens_files)
 

webviz_subsurface/_models/__init__.py

@@ -1,3 +1,4 @@
 from .ensemble_model import EnsembleModel
 from .ensemble_set_model import EnsembleSetModel
 from .surface_leaflet_model import SurfaceLeafletModel
+from .surface_set_model import SurfaceSetModel

webviz_subsurface/_models/surface_set_model.py

@@ -0,0 +1,219 @@
+from typing import List, Tuple, Callable, Optional, Any
+from pathlib import Path
+import warnings
+import json
+import io
+
+import numpy as np
+import pandas as pd
+import xtgeo
+from webviz_config.webviz_store import webvizstore
+from webviz_config.common_cache import CACHE
+
+
+class SurfaceSetModel:
+    """Class to load and calculate statistical surfaces from an FMU Ensemble"""
+
+    def __init__(self, surface_table: pd.DataFrame):
+        self._surface_table = surface_table
+
+    @property
+    def attributes(self) -> list:
+        """Returns surface attributes"""
+        return list(self._surface_table["attribute"].unique())
+
+    def names_in_attribute(self, attribute: str) -> list:
+        """Returns surface names for a given attribute"""
+        return list(
+            self._surface_table.loc[self._surface_table["attribute"] == attribute][
+                "name"
+            ].unique()
+        )
+
+    def dates_in_attribute(self, attribute: str) -> list:
+        """Returns surface dates for a given attribute"""
+        return list(
+            self._surface_table.loc[self._surface_table["attribute"] == attribute][
+                "date"
+            ].unique()
+        )
+
+    def get_realization_surface(
+        self, name: str, attribute: str, realization: int, date: Optional[str] = None
+    ) -> xtgeo.RegularSurface:
+        """Returns a Xtgeo surface instance of a single realization surface"""
+
+        columns = ["name", "attribute", "REAL"]
+        column_values = [name, attribute, realization]
+        if date is not None:
+            columns.append("date")
+            column_values.append(date)
+
+        df = self._filter_surface_table(
+            name=name, attribute=attribute, date=date, realizations=[realization]
+        )
+        if len(df.index) == 0:
+            warnings.warn(
+                f"No surface found for name: {name}, attribute: {attribute}, date: {date}, "
+                f"realization: {realization}"
+            )
+            return xtgeo.RegularSurface()
+        if len(df.index) > 1:
+            warnings.warn(
+                f"Multiple surfaces found for name: {name}, attribute: {attribute}, date: {date}, "
+                f"realization: {realization}. Returning first surface"
+            )
+        return xtgeo.surface_from_file(get_stored_surface_path(df.iloc[0]["path"]))
+
+    def _filter_surface_table(
+        self,
+        name: str,
+        attribute: str,
+        date: Optional[str] = None,
+        realizations: Optional[List[int]] = None,
+    ) -> pd.DataFrame:
+        """Returns a dataframe of surfaces for the provided filters"""
+        columns: List[str] = ["name", "attribute"]
+        column_values: List[Any] = [name, attribute]
+        if date is not None:
+            columns.append("date")
+            column_values.append(date)
+        if realizations is not None:
+            columns.append("REAL")
+            column_values.append(realizations)
+        df = self._surface_table.copy()
+        for filt, col in zip(column_values, columns):
+            if isinstance(filt, list):
+                df = df.loc[df[col].isin(filt)]
+            else:
+                df = df.loc[df[col] == filt]
+        return df
+
+    @CACHE.memoize(timeout=CACHE.TIMEOUT)
+    def calculate_statistical_surface(
+        self,
+        name: str,
+        attribute: str,
+        calculation: Optional[str] = "Mean",
+        date: Optional[str] = None,
+        realizations: Optional[List[int]] = None,
+    ) -> xtgeo.RegularSurface:
+        """Returns a Xtgeo surface instance for a calculated surface"""
+        df = self._filter_surface_table(
+            name=name, attribute=attribute, date=date, realizations=realizations
+        )
+        return surface_from_json(
+            json.load(save_statistical_surface(sorted(list(df["path"])), calculation))
+        )
+
+    def webviz_store_statistical_calculation(
+        self,
+        calculation: Optional[str] = "Mean",
+        realizations: Optional[List[int]] = None,
+    ) -> Tuple[Callable, list]:
+        """Returns a tuple of functions to calculate statistical surfaces for
+        webviz store"""
+        df = (
+            self._surface_table.loc[self._surface_table["REAL"].isin(realizations)]
+            if realizations is not None
+            else self._surface_table
+        )
+        stored_functions_args = []
+        for _attr, attr_df in df.groupby("attribute"):
+            for _name, name_df in attr_df.groupby("name"):
+
+                if name_df["date"].isnull().values.all():
+                    stored_functions_args.append(
+                        {
+                            "fns": sorted(list(name_df["path"].unique())),
+                            "calculation": calculation,
+                        }
+                    )
+                else:
+                    for _date, date_df in name_df.groupby("date"):
+                        stored_functions_args.append(
+                            {
+                                "fns": sorted(list(date_df["path"].unique())),
+                                "calculation": calculation,
+                            }
+                        )
+
+        return (
+            save_statistical_surface,
+            stored_functions_args,
+        )
+
+    def webviz_store_realization_surfaces(self) -> Tuple[Callable, list]:
+        """Returns a tuple of functions to store all realization surfaces for
+        webviz store"""
+        return (
+            get_stored_surface_path,
+            [{"runpath": path} for path in list(self._surface_table["path"])],
+        )
+
+
+@webvizstore
+def get_stored_surface_path(runpath: Path) -> Path:
+    """Returns path of a stored surface"""
+    return Path(runpath)
+
+
+@webvizstore
+def save_statistical_surface(fns: List[str], calculation: str) -> io.BytesIO:
+    """Wrapper function to store a calculated surface as BytesIO"""
+    surfaces = xtgeo.Surfaces(fns)
+    if len(surfaces.surfaces) == 0:
+        surface = xtgeo.RegularSurface()
+    elif calculation in ["Mean", "StdDev", "Min", "Max", "P10", "P90"]:
+        # Suppress numpy warnings when surfaces have undefined z-values
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", "All-NaN slice encountered")
+            warnings.filterwarnings("ignore", "Mean of empty slice")
+            warnings.filterwarnings("ignore", "Degrees of freedom <= 0 for slice")
+            surface = get_statistical_surface(surfaces, calculation)
+    else:
+        surface = xtgeo.RegularSurface()
+    return io.BytesIO(surface_to_json(surface).encode())
+
+
+# pylint: disable=too-many-return-statements
+def get_statistical_surface(
+    surfaces: xtgeo.Surfaces, calculation: str
+) -> xtgeo.RegularSurface:
+    """Calculates a statistical surface from a list of Xtgeo surface instances"""
+    if calculation == "Mean":
+        return surfaces.apply(np.nanmean, axis=0)
+    if calculation == "StdDev":
+        return surfaces.apply(np.nanstd, axis=0)
+    if calculation == "Min":
+        return surfaces.apply(np.nanmin, axis=0)
+    if calculation == "Max":
+        return surfaces.apply(np.nanmax, axis=0)
+    if calculation == "P10":
+        return surfaces.apply(np.nanpercentile, 10, axis=0)
+    if calculation == "P90":
+        return surfaces.apply(np.nanpercentile, 90, axis=0)
+    return xtgeo.RegularSurface()
+
+
+def surface_to_json(surface: xtgeo.RegularSurface) -> str:
+    """Returns a json represention of a Xtgeo surface instance"""
+    return json.dumps(
+        {
+            "ncol": surface.ncol,
+            "nrow": surface.nrow,
+            "xori": surface.xori,
+            "yori": surface.yori,
+            "rotation": surface.rotation,
+            "xinc": surface.xinc,
+            "yinc": surface.yinc,
+            "values": surface.values.copy().filled(np.nan).tolist(),
+        }
+    )
+
+
+def surface_from_json(surfaceobj: dict) -> xtgeo.RegularSurface:
+    """Returns a Xtgeo surface instance from a json surface representation"""
+    surface = xtgeo.RegularSurface(**surfaceobj)
+    surface.values = surfaceobj["values"]
+    return surface

webviz_subsurface/_private_plugins/surface_selector.py

@@ -330,7 +330,7 @@ def _set_data(
             dates_in_attr = self._dates_in_attr(attr)
             if dates_in_attr and date and not date in dates_in_attr:
                 raise PreventUpdate
-            return json.dumps({"name": name, "attr": attr, "date": date})
+            return json.dumps({"name": name, "attribute": attr, "date": date})
 
 
 def prev_value(current_value: str, options: List[str]) -> str: