diff --git a/.github/workflows/subsurface.yml b/.github/workflows/subsurface.yml index d491fd1a1..a59b76003 100644 --- a/.github/workflows/subsurface.yml +++ b/.github/workflows/subsurface.yml @@ -10,7 +10,7 @@ on: - published schedule: # Run CI daily and check that tests are working with latest dependencies - - cron: '0 0 * * *' + - cron: "0 0 * * *" jobs: webviz-subsurface: @@ -18,134 +18,133 @@ jobs: if: github.event_name != 'push' || github.ref == 'refs/heads/master' || contains(github.event.head_commit.message, '[deploy test]') runs-on: ubuntu-latest env: - PYTHONWARNINGS: default # We want to see e.g. DeprecationWarnings + PYTHONWARNINGS: default # We want to see e.g. DeprecationWarnings strategy: fail-fast: false matrix: - python-version: ['3.6', '3.7', '3.8', '3.9', '3.10'] + python-version: ["3.6", "3.7", "3.8", "3.9", "3.10"] steps: - - - name: ๐Ÿงน Remove unused pre-installed software - run: | - # https://github.com/actions/virtual-environments/issues/751 - # https://github.com/actions/virtual-environments/issues/709 - sudo apt-get purge p7zip* yarn ruby-full ghc* php7* - sudo apt-get autoremove - sudo apt-get clean - df -h - - - name: ๐Ÿ“– Checkout commit locally - uses: actions/checkout@v2 - - - name: ๐Ÿ Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v2 - with: - python-version: ${{ matrix.python-version }} - - - name: ๐Ÿ“ฆ Install webviz-subsurface with dependencies - run: | - pip install --upgrade pip - if [[ $(pip freeze) ]]; then - pip freeze | grep -vw "pip" | xargs pip uninstall -y - fi - pip install "bleach<5" # https://github.com/equinor/webviz-config/issues/586 - pip install "werkzeug<2.1" # ...while waiting for https://github.com/plotly/dash/issues/1992 - pip install . - - # Testing against our latest release (including pre-releases) - # Temporarily do not use prerelase of wsc due to https://github.com/equinor/webviz-subsurface-components/issues/1007 - pip install --pre --upgrade webviz-config webviz-core-components # webviz-subsurface-components - - - name: ๐Ÿ“ฆ Install test dependencies - run: | - pip install .[tests] - wget https://chromedriver.storage.googleapis.com/$(wget https://chromedriver.storage.googleapis.com/LATEST_RELEASE -q -O -)/chromedriver_linux64.zip - unzip chromedriver_linux64.zip - export PATH=$PATH:$PWD - - - name: ๐Ÿงพ List all installed packages - run: pip freeze - - - name: ๐Ÿ•ต๏ธ Check code style & linting - run: | - black --check webviz_subsurface tests setup.py - pylint webviz_subsurface tests setup.py - bandit -r -c ./bandit.yml webviz_subsurface tests setup.py - isort --check-only webviz_subsurface tests setup.py - mypy --package webviz_subsurface - - - name: ๐Ÿค– Run tests - env: - # If you want the CI to (temporarily) run against your fork of the testdada, - # change the value her from "equinor" to your username. - TESTDATA_REPO_OWNER: equinor - # If you want the CI to (temporarily) run against another branch than master, - # change the value her from "master" to the relevant branch name. - TESTDATA_REPO_BRANCH: master - run: | - git clone --depth 1 --branch $TESTDATA_REPO_BRANCH https://github.com/$TESTDATA_REPO_OWNER/webviz-subsurface-testdata.git - # Copy any clientside script to the test folder before running tests - mkdir ./tests/assets && cp ./webviz_subsurface/_assets/js/* ./tests/assets - pytest ./tests --headless --forked --testdata-folder ./webviz-subsurface-testdata - rm -rf ./tests/assets - webviz docs --portable ./docs_build --skip-open - - - name: ๐Ÿณ Build Docker example image - if: matrix.python-version != '3.7' # https://github.com/statsmodels/statsmodels/issues/8110 - run: | - pip install --pre webviz-config-equinor - export SOURCE_URL_WEBVIZ_SUBSURFACE=https://github.com/$GITHUB_REPOSITORY - export GIT_POINTER_WEBVIZ_SUBSURFACE=$GITHUB_REF - webviz build ./webviz-subsurface-testdata/webviz_examples/webviz-full-demo.yml --portable ./example_subsurface_app --theme equinor - rm -rf ./webviz-subsurface-testdata - pushd example_subsurface_app - docker build -t webviz/example_subsurface_image:equinor-theme . - popd - - - name: ๐Ÿณ Update Docker Hub example image - if: github.event_name != 'schedule' && github.ref == 'refs/heads/master' && matrix.python-version == '3.6' - run: | - echo ${{ secrets.dockerhub_webviz_token }} | docker login --username webviz --password-stdin - docker push webviz/example_subsurface_image:equinor-theme - - - name: ๐Ÿณ Update review/test Docker example image - if: github.ref != 'refs/heads/master' && contains(github.event.head_commit.message, '[deploy test]') && matrix.python-version == '3.6' - run: | - docker tag webviz/example_subsurface_image:equinor-theme ${{ secrets.review_docker_registry_url }}/${{ secrets.review_container_name }} - - echo ${{ secrets.review_docker_registry_token }} | docker login ${{ secrets.review_docker_registry_url }} --username ${{ secrets.review_docker_registry_username }} --password-stdin - docker push ${{ secrets.review_docker_registry_url }}/${{ secrets.review_container_name }} - - - name: ๐Ÿšข Build and deploy Python package - if: github.event_name == 'release' && matrix.python-version == '3.6' - env: - TWINE_USERNAME: __token__ - TWINE_PASSWORD: ${{ secrets.pypi_webviz_token }} - run: | - python -m pip install --upgrade setuptools wheel twine - python setup.py sdist bdist_wheel - twine upload dist/* - - - name: ๐Ÿ“š Update GitHub pages - if: github.event_name == 'release' && matrix.python-version == '3.6' - run: | - cp -R ./docs_build ../docs_build - - git config --local user.email "webviz-github-action" - git config --local user.name "webviz-github-action" - git fetch origin gh-pages - git checkout --track origin/gh-pages - git clean -f -f -d -x - git rm -r * - - cp -R ../docs_build/* . - - git add . - - if git diff-index --quiet HEAD; then - echo "No changes in documentation. Skip documentation deploy." - else - git commit -m "Update Github Pages" - git push "https://${{ github.actor }}:${{ secrets.GITHUB_TOKEN }}@github.com/${{ github.repository }}.git" gh-pages - fi + - name: ๐Ÿงน Remove unused pre-installed software + run: | + # https://github.com/actions/virtual-environments/issues/751 + # https://github.com/actions/virtual-environments/issues/709 + sudo apt-get purge p7zip* yarn ruby-full ghc* php7* + sudo apt-get autoremove + sudo apt-get clean + df -h + + - name: ๐Ÿ“– Checkout commit locally + uses: actions/checkout@v2 + + - name: ๐Ÿ Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + + - name: ๐Ÿ“ฆ Install webviz-subsurface with dependencies + run: | + pip install --upgrade pip + if [[ $(pip freeze) ]]; then + pip freeze | grep -vw "pip" | xargs pip uninstall -y + fi + pip install "bleach<5" # https://github.com/equinor/webviz-config/issues/586 + pip install "werkzeug<2.1" # ...while waiting for https://github.com/plotly/dash/issues/1992 + pip install . + + # Testing against our latest release (including pre-releases) + # Temporarily do not use prerelase of wsc due to https://github.com/equinor/webviz-subsurface-components/issues/1007 + pip install --pre --upgrade webviz-config webviz-core-components # webviz-subsurface-components + + - name: ๐Ÿ“ฆ Install test dependencies + run: | + pip install .[tests] + wget https://chromedriver.storage.googleapis.com/$(wget https://chromedriver.storage.googleapis.com/LATEST_RELEASE -q -O -)/chromedriver_linux64.zip + unzip chromedriver_linux64.zip + export PATH=$PATH:$PWD + + - name: ๐Ÿงพ List all installed packages + run: pip freeze + + - name: ๐Ÿ•ต๏ธ Check code style & linting + run: | + black --check webviz_subsurface tests setup.py + pylint webviz_subsurface tests setup.py + bandit -r -c ./bandit.yml webviz_subsurface tests setup.py + isort --check-only webviz_subsurface tests setup.py + mypy --package webviz_subsurface + + - name: ๐Ÿค– Run tests + env: + # If you want the CI to (temporarily) run against your fork of the testdada, + # change the value her from "equinor" to your username. + TESTDATA_REPO_OWNER: hanskallekleiv + # If you want the CI to (temporarily) run against another branch than master, + # change the value her from "master" to the relevant branch name. + TESTDATA_REPO_BRANCH: more-grid + run: | + git clone --depth 1 --branch $TESTDATA_REPO_BRANCH https://github.com/$TESTDATA_REPO_OWNER/webviz-subsurface-testdata.git + # Copy any clientside script to the test folder before running tests + mkdir ./tests/assets && cp ./webviz_subsurface/_assets/js/* ./tests/assets + pytest ./tests --headless --forked --testdata-folder ./webviz-subsurface-testdata + rm -rf ./tests/assets + webviz docs --portable ./docs_build --skip-open + + - name: ๐Ÿณ Build Docker example image + if: matrix.python-version != '3.7' # https://github.com/statsmodels/statsmodels/issues/8110 + run: | + pip install --pre webviz-config-equinor + export SOURCE_URL_WEBVIZ_SUBSURFACE=https://github.com/$GITHUB_REPOSITORY + export GIT_POINTER_WEBVIZ_SUBSURFACE=$GITHUB_REF + webviz build ./webviz-subsurface-testdata/webviz_examples/webviz-full-demo.yml --portable ./example_subsurface_app --theme equinor + rm -rf ./webviz-subsurface-testdata + pushd example_subsurface_app + docker build -t webviz/example_subsurface_image:equinor-theme . + popd + + - name: ๐Ÿณ Update Docker Hub example image + if: github.event_name != 'schedule' && github.ref == 'refs/heads/master' && matrix.python-version == '3.6' + run: | + echo ${{ secrets.dockerhub_webviz_token }} | docker login --username webviz --password-stdin + docker push webviz/example_subsurface_image:equinor-theme + + - name: ๐Ÿณ Update review/test Docker example image + if: github.ref != 'refs/heads/master' && contains(github.event.head_commit.message, '[deploy test]') && matrix.python-version == '3.6' + run: | + docker tag webviz/example_subsurface_image:equinor-theme ${{ secrets.review_docker_registry_url }}/${{ secrets.review_container_name }} + + echo ${{ secrets.review_docker_registry_token }} | docker login ${{ secrets.review_docker_registry_url }} --username ${{ secrets.review_docker_registry_username }} --password-stdin + docker push ${{ secrets.review_docker_registry_url }}/${{ secrets.review_container_name }} + + - name: ๐Ÿšข Build and deploy Python package + if: github.event_name == 'release' && matrix.python-version == '3.6' + env: + TWINE_USERNAME: __token__ + TWINE_PASSWORD: ${{ secrets.pypi_webviz_token }} + run: | + python -m pip install --upgrade setuptools wheel twine + python setup.py sdist bdist_wheel + twine upload dist/* + + - name: ๐Ÿ“š Update GitHub pages + if: github.event_name == 'release' && matrix.python-version == '3.6' + run: | + cp -R ./docs_build ../docs_build + + git config --local user.email "webviz-github-action" + git config --local user.name "webviz-github-action" + git fetch origin gh-pages + git checkout --track origin/gh-pages + git clean -f -f -d -x + git rm -r * + + cp -R ../docs_build/* . + + git add . + + if git diff-index --quiet HEAD; then + echo "No changes in documentation. Skip documentation deploy." + else + git commit -m "Update Github Pages" + git push "https://${{ github.actor }}:${{ secrets.GITHUB_TOKEN }}@github.com/${{ github.repository }}.git" gh-pages + fi diff --git a/setup.py b/setup.py index efc576e0b..bef3292cb 100644 --- a/setup.py +++ b/setup.py @@ -40,6 +40,8 @@ "webviz_config_plugins": [ "BhpQc = webviz_subsurface.plugins:BhpQc", "DiskUsage = webviz_subsurface.plugins:DiskUsage", + "EclipseGridViewer = webviz_subsurface.plugins:EclipseGridViewer", + "GridViewer = webviz_subsurface.plugins:GridViewer", "GroupTree = webviz_subsurface.plugins:GroupTree", "HistoryMatch = webviz_subsurface.plugins:HistoryMatch", "HorizonUncertaintyViewer = webviz_subsurface.plugins:HorizonUncertaintyViewer", @@ -87,6 +89,7 @@ "dash>=2.0.0", "dash_bootstrap_components>=0.10.3", "dash-daq>=0.5.0", + # "dash-vtk>=0.0.9", "dataclasses>=0.8; python_version<'3.7'", "defusedxml>=0.6.0", "ecl2df>=0.15.0; sys_platform=='linux'", @@ -99,12 +102,16 @@ "pillow>=6.1", "pyarrow>=5.0.0", "pyscal>=0.7.5", + "pyvista>=0.33.3", "scipy>=1.2", "statsmodels>=0.12.1", # indirect dependency through https://plotly.com/python/linear-fits/ - "webviz-config>=0.3.8", - "webviz-core-components>=0.5.6", + "vtk>=9.2.0rc2", + "webviz_config@git+https://github.com/equinor/webviz-config@50f0d20f5b6fd6c2d7ee98b17e229e625776ff82", + "webviz-core-components", "webviz-subsurface-components>=0.4.12", - "xtgeo>=2.14", + "webviz_vtk@git+https://github.com/equinor/webviz-vtk", + "xtgeo@git+https://github.com/sigurdp/xtgeo/@sigurdp/vtk-esg", + # "xtgeo>=2.18.0a1", ], extras_require={"tests": TESTS_REQUIRE}, setup_requires=["setuptools_scm~=3.2"], diff --git a/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/__init__.py b/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/_utils.py b/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/_utils.py new file mode 100644 index 000000000..0dd0598aa --- /dev/null +++ b/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/_utils.py @@ -0,0 +1,36 @@ +# pylint: skip-file +# type: ignore +import numpy as np +import pyvista as pv + + +def create_explicit_structured_grid( + ni: int, nj: int, nk: int, si: float, sj: float, sk: float +) -> pv.ExplicitStructuredGrid: + + si = float(si) + sj = float(sj) + sk = float(sk) + + # create raw coordinate grid + grid_ijk = np.mgrid[ + : (ni + 1) * si : si, : (nj + 1) * sj : sj, : (nk + 1) * sk : sk + ] + + # repeat array along each Cartesian axis for connectivity + for axis in range(1, 4): + grid_ijk = grid_ijk.repeat(2, axis=axis) + + # slice off unnecessarily doubled edge coordinates + grid_ijk = grid_ijk[:, 1:-1, 1:-1, 1:-1] + + # reorder and reshape to VTK order + corners = grid_ijk.transpose().reshape(-1, 3) + + dims = np.array([ni, nj, nk]) + 1 + + grid = pv.ExplicitStructuredGrid(dims, corners) + grid = grid.compute_connectivity() + grid.ComputeFacesConnectivityFlagsArray() + + return grid diff --git a/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/test_explicit_structured_grid_accessor.py b/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/test_explicit_structured_grid_accessor.py new file mode 100644 index 000000000..28d03d4f9 --- /dev/null +++ b/tests/unit_tests/plugin_tests/test_eclipse_grid_viewer/test_explicit_structured_grid_accessor.py @@ -0,0 +1,96 @@ +# pylint: skip-file +# type: ignore +import pytest +from vtk.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import vtkCellLocator, vtkExplicitStructuredGrid + +from webviz_subsurface.plugins._eclipse_grid_viewer._explicit_structured_grid_accessor import ( + ExplicitStructuredGridAccessor, +) + +from ._utils import create_explicit_structured_grid + +ES_GRID_ACCESSOR = ExplicitStructuredGridAccessor( + create_explicit_structured_grid(5, 4, 3, 20.0, 10.0, 5.0) +) + +CROP_FIRST_CELL = [0, 0], [0, 0], [0, 0] +EXPECTED_FIRST_CELL_ORIGINAL_INDEX = [0] +CROP_LAST_CELL = [5, 5], [4, 4], [3, 3] +EXPECTED_LAST_CELL_ORIGINAL_INDEX = [59] +CROP_BOX = [2, 3], [2, 3], [1, 2] +EXPECTED_CROP_BOX_ORIGINAL_INDEX = [32, 33, 37, 38, 52, 53, 57, 58] + + +@pytest.mark.parametrize( + "crop_range, expected_cells", + [ + (CROP_FIRST_CELL, EXPECTED_FIRST_CELL_ORIGINAL_INDEX), + (CROP_LAST_CELL, EXPECTED_LAST_CELL_ORIGINAL_INDEX), + (CROP_BOX, EXPECTED_CROP_BOX_ORIGINAL_INDEX), + ], +) +def test_crop(crop_range, expected_cells) -> None: + cropped_grid = ES_GRID_ACCESSOR.crop(*crop_range) + assert isinstance(cropped_grid, vtkExplicitStructuredGrid) + assert cropped_grid.GetCellData().HasArray("vtkOriginalCellIds") == 1 + assert set( + vtk_to_numpy(cropped_grid.GetCellData().GetAbstractArray("vtkOriginalCellIds")) + ) == set(expected_cells) + + _polys, _points, indices = ES_GRID_ACCESSOR.extract_skin(cropped_grid) + assert set(indices) == set(expected_cells) + + +RAY_FROM_TOP = [ + [50, 15, 15], + [50, 15, -5], +] +RAY_FROM_BOTTOM = [ + [50, 15, -5], + [50, 15, 20], +] +RAY_FROM_I = [[50, -7, 13], [50, 45, 13]] +RAY_FROM_J = [[-12, 5, 13], [110, 5, 13]] + + +@pytest.mark.parametrize( + "ray, expected_cell_id_and_ijk", + [ + (RAY_FROM_TOP, (47, [2, 1, 2])), + (RAY_FROM_BOTTOM, (7, [2, 1, 0])), + (RAY_FROM_I, (42, [2, 0, 2])), + (RAY_FROM_J, (40, [0, 0, 2])), + ], +) +def test_find_closest_cell_to_ray(ray, expected_cell_id_and_ijk) -> None: + cell_ijk = ES_GRID_ACCESSOR.find_closest_cell_to_ray(ES_GRID_ACCESSOR.es_grid, ray) + assert cell_ijk == expected_cell_id_and_ijk + + +@pytest.mark.parametrize( + "ray, expected_cell_id_and_ijk", + [ + (RAY_FROM_TOP, (27, [2, 1, 1])), + (RAY_FROM_BOTTOM, (7, [2, 1, 0])), + (RAY_FROM_I, (52, [2, 2, 2])), + (RAY_FROM_J, (43, [3, 0, 2])), + ], +) +def test_find_closest_cell_to_ray_with_blanked_cells( + ray, expected_cell_id_and_ijk +) -> None: + grid = ES_GRID_ACCESSOR.es_grid.copy() + + cellLocator = vtkCellLocator() + cellLocator.SetDataSet(grid) + cellLocator.BuildLocator() + cellIds = vtkIdList() + cellLocator.FindCellsAlongLine((6.0, 6.0, 12.0), (67.0, 12.0, 12.0), 0.001, cellIds) + for i in range(cellIds.GetNumberOfIds()): + id = cellIds.GetId(i) + grid.BlankCell(id) + + cell_ijk = ES_GRID_ACCESSOR.find_closest_cell_to_ray(grid, ray) + assert cell_ijk == expected_cell_id_and_ijk diff --git a/tmp_dashvtk/dash_vtk-0.0.9.tar.gz b/tmp_dashvtk/dash_vtk-0.0.9.tar.gz new file mode 100644 index 000000000..57f3e96a6 Binary files /dev/null and b/tmp_dashvtk/dash_vtk-0.0.9.tar.gz differ diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/__init__.py b/webviz_subsurface/_providers/ensemble_grid_provider/__init__.py new file mode 100644 index 000000000..3519237a4 --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/__init__.py @@ -0,0 +1,3 @@ +from .ensemble_grid_provider import EnsembleGridProvider +from .ensemble_grid_provider_factory import EnsembleGridProviderFactory +from .grid_viz_service import GridVizService, CellFilter, PropertySpec, Ray, PickResult diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/_roff_file_discovery.py b/webviz_subsurface/_providers/ensemble_grid_provider/_roff_file_discovery.py new file mode 100644 index 000000000..57eb36f0f --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/_roff_file_discovery.py @@ -0,0 +1,117 @@ +import glob +import os +import re +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, List, Optional, Union, Tuple + +from fmu.ensemble import ScratchEnsemble + + +@dataclass(frozen=True) +class GridParameterFileInfo: + path: str + real: int + name: str + attribute: str + datestr: Optional[str] + + +@dataclass(frozen=True) +class GridParameterIdent: + name: str + attribute: str + datestr: Optional[str] + + +@dataclass(frozen=True) +class GridFileInfo: + path: str + real: int + name: str + + +@dataclass(frozen=True) +class GridIdent: + name: str + + +def _discover_ensemble_realizations_fmu(ens_path: str) -> Dict[int, str]: + """Returns dict indexed by realization number and with runpath as value""" + scratch_ensemble = ScratchEnsemble("dummyEnsembleName", paths=ens_path).filter("OK") + real_dict = {i: r.runpath() for i, r in scratch_ensemble.realizations.items()} + return real_dict + + +def _discover_ensemble_realizations(ens_path: str) -> Dict[int, str]: + # Much faster than FMU impl above, but is it risky? + # Do we need to check for OK-file? + real_dict: Dict[int, str] = {} + + realidxregexp = re.compile(r"realization-(\d+)") + globbed_real_dirs = sorted(glob.glob(str(ens_path))) + for real_dir in globbed_real_dirs: + realnum: Optional[int] = None + for path_comp in reversed(real_dir.split(os.path.sep)): + realmatch = re.match(realidxregexp, path_comp) + if realmatch: + realnum = int(realmatch.group(1)) + break + + if realnum is not None: + real_dict[realnum] = real_dir + + return real_dict + + +def ident_from_filename( + filename: str, +) -> Union[GridIdent, GridParameterIdent]: + """Split the stem part of the roff filename into grid name, attribute and + optionally date part""" + delimiter: str = "--" + parts = Path(filename).stem.split(delimiter) + if len(parts) == 1: + return GridIdent(name=parts[0]) + + return GridParameterIdent( + name=parts[0], attribute=parts[1], datestr=parts[2] if len(parts) >= 3 else None + ) + + +def discover_per_realization_grid_files( + ens_path: str, grid_name: str, attribute_filter: List[str] = None +) -> Tuple[List[GridFileInfo], List[GridParameterFileInfo]]: + rel_folder: str = "share/results/grids" + suffix: str = "*.roff" + + grid_parameters_info: List[GridParameterFileInfo] = [] + grid_info: List[GridFileInfo] = [] + real_dict = _discover_ensemble_realizations_fmu(ens_path) + for realnum, runpath in sorted(real_dict.items()): + globbed_filenames = glob.glob(str(Path(runpath) / rel_folder / suffix)) + for filename in sorted(globbed_filenames): + ident = ident_from_filename(filename) + if ident.name != grid_name: + continue + if isinstance(ident, GridParameterIdent): + if ( + attribute_filter is not None + and ident.attribute not in attribute_filter + ): + continue + grid_parameters_info.append( + GridParameterFileInfo( + path=filename, + real=realnum, + name=ident.name, + attribute=ident.attribute, + datestr=ident.datestr, + ) + ) + else: + grid_info.append( + GridFileInfo(path=filename, real=realnum, name=ident.name) + ) + + return grid_info, grid_parameters_info diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/_xtgeo_to_vtk_explicit_structured_grid.py b/webviz_subsurface/_providers/ensemble_grid_provider/_xtgeo_to_vtk_explicit_structured_grid.py new file mode 100644 index 000000000..4da1324eb --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/_xtgeo_to_vtk_explicit_structured_grid.py @@ -0,0 +1,101 @@ +import xtgeo +import numpy as np + +# pylint: disable=no-name-in-module, +from vtkmodules.vtkCommonDataModel import vtkExplicitStructuredGrid +from vtkmodules.vtkCommonDataModel import vtkCellArray +from vtkmodules.vtkCommonDataModel import vtkDataSetAttributes +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.util.numpy_support import numpy_to_vtkIdTypeArray +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.util import vtkConstants + +from webviz_subsurface._utils.perf_timer import PerfTimer + +# from ._xtgeo_to_explicit_structured_grid_hack import ( +# _clean_vtk_ug, +# _vtk_esg_to_ug, +# _vtk_ug_to_esg, +# ) + + +# ----------------------------------------------------------------------------- +def xtgeo_grid_to_vtk_explicit_structured_grid( + xtg_grid: xtgeo.Grid, +) -> vtkExplicitStructuredGrid: + + print("entering xtgeo_grid_to_vtk_explicit_structured_grid()") + t = PerfTimer() + + pt_dims, vertex_arr, conn_arr, inactive_arr = xtg_grid.get_vtk_esg_geometry_data() + vertex_arr[:, 2] *= -1 + print(f"get_vtk_esg_geometry_data() took {t.lap_s():.2f}s") + + print(f"pt_dims={pt_dims}") + print(f"vertex_arr.shape={vertex_arr.shape}") + print(f"vertex_arr.dtype={vertex_arr.dtype}") + print(f"conn_arr.shape={conn_arr.shape}") + print(f"conn_arr.dtype={conn_arr.dtype}") + + vtk_esgrid = _create_vtk_esgrid_from_verts_and_conn(pt_dims, vertex_arr, conn_arr) + print(f"create vtk_esgrid : {t.lap_s():.2f}s") + + # Make sure we hide the inactive cells. + # First we let VTK allocate cell ghost array, then we obtain a numpy view + # on the array and write to that (we're actually modifying the native VTK array) + ghost_arr_vtk = vtk_esgrid.AllocateCellGhostArray() + ghost_arr_np = vtk_to_numpy(ghost_arr_vtk) + ghost_arr_np[inactive_arr] = vtkDataSetAttributes.HIDDENCELL + print(f"hide {len(inactive_arr)} inactive cells : {t.lap_s():.2f}s") + + print(f"memory used by vtk_esgrid: {vtk_esgrid.GetActualMemorySize()/1024.0:.2f}MB") + + print(f"xtgeo_grid_to_vtk_explicit_structured_grid() - DONE: {t.elapsed_s():.2f}s") + + # print("==================================================================") + # print(pv.ExplicitStructuredGrid(vtk_esgrid)) + # print("==================================================================") + + return vtk_esgrid + + +# ----------------------------------------------------------------------------- +def _create_vtk_esgrid_from_verts_and_conn( + point_dims: np.ndarray, vertex_arr_np: np.ndarray, conn_arr_np: np.ndarray +) -> vtkExplicitStructuredGrid: + + print("_create_vtk_esgrid_from_verts_and_conn() - entering") + + t = PerfTimer() + + vertex_arr_np = vertex_arr_np.reshape(-1, 3) + points_vtkarr = numpy_to_vtk(vertex_arr_np, deep=1) + vtk_points = vtkPoints() + vtk_points.SetData(points_vtkarr) + print(f"_create_vtk_esgrid_from_verts_and_conn() - vtk_points: {t.lap_s():.2f}s") + + # conn_idarr = numpy_to_vtk(conn_arr_np, deep=1, array_type=vtkConstants.VTK_ID_TYPE) + conn_idarr = numpy_to_vtkIdTypeArray(conn_arr_np, deep=1) + vtk_cellArray = vtkCellArray() + vtk_cellArray.SetData(8, conn_idarr) + print(f"_create_vtk_esgrid_from_verts_and_conn() - vtk_cellArray: {t.lap_s():.2f}s") + + vtk_esgrid = vtkExplicitStructuredGrid() + vtk_esgrid.SetDimensions(point_dims) + vtk_esgrid.SetPoints(vtk_points) + vtk_esgrid.SetCells(vtk_cellArray) + print(f"_create_vtk_esgrid_from_verts_and_conn() - vtk_esgrid: {t.lap_s():.2f}s") + + vtk_esgrid.ComputeFacesConnectivityFlagsArray() + print(f"_create_vtk_esgrid_from_verts_and_conn() - conn flags: {t.lap_s():.2f}s") + + # print(pv.ExplicitStructuredGrid(vtk_esgrid)) + # ugrid = _vtk_esg_to_ug(vtk_esgrid) + # ugrid = _clean_vtk_ug(ugrid) + # vtk_esgrid = _vtk_ug_to_esg(ugrid) + # print(pv.ExplicitStructuredGrid(vtk_esgrid)) + + print(f"_create_vtk_esgrid_from_verts_and_conn() - DONE: {t.elapsed_s():.2f}s") + + return vtk_esgrid diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider.py b/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider.py new file mode 100644 index 000000000..ee9e18509 --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider.py @@ -0,0 +1,46 @@ +import abc +from typing import List, Optional + +import xtgeo +import numpy as np + + +class EnsembleGridProvider(abc.ABC): + @abc.abstractmethod + def provider_id(self) -> str: + """Returns string ID of the provider.""" + + @abc.abstractmethod + def static_property_names(self) -> List[str]: + """Returns list of all available static properties.""" + + @abc.abstractmethod + def dynamic_property_names(self) -> List[str]: + """Returns list of all available dynamic properties.""" + + @abc.abstractmethod + def dates_for_dynamic_property(self, property_name: str) -> Optional[List[str]]: + """Returns list of all available dates for a given dynamic property.""" + + @abc.abstractmethod + def realizations(self) -> List[int]: + """Returns list of all available realizations.""" + + @abc.abstractmethod + def get_3dgrid( + self, + realization: int, + ) -> Optional[xtgeo.Grid]: + """Returns grid for specified realization""" + + @abc.abstractmethod + def get_static_property_values( + self, property_name: str, realization: int + ) -> Optional[np.ndarray]: + """Returns 1d cell values for a given static property""" + + @abc.abstractmethod + def get_dynamic_property_values( + self, property_name: str, property_date: str, realization: int + ) -> Optional[np.ndarray]: + """Returns 1d cell values for a given dynamic property""" diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider_factory.py b/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider_factory.py new file mode 100644 index 000000000..acbd7abee --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/ensemble_grid_provider_factory.py @@ -0,0 +1,116 @@ +import hashlib +import logging +import os +from pathlib import Path +from typing import List + +from webviz_config.webviz_factory import WebvizFactory +from webviz_config.webviz_factory_registry import WEBVIZ_FACTORY_REGISTRY +from webviz_config.webviz_instance_info import WebvizRunMode + +from webviz_subsurface._utils.perf_timer import PerfTimer + +from .provider_impl_roff import ProviderImplRoff +from ._roff_file_discovery import discover_per_realization_grid_files +from .ensemble_grid_provider import EnsembleGridProvider + +LOGGER = logging.getLogger(__name__) + + +class EnsembleGridProviderFactory(WebvizFactory): + def __init__( + self, + root_storage_folder: Path, + allow_storage_writes: bool, + avoid_copying_grid_data: bool, + ) -> None: + self._storage_dir = Path(root_storage_folder) / __name__ + self._allow_storage_writes = allow_storage_writes + self._avoid_copying_grid_data = avoid_copying_grid_data + + LOGGER.info( + f"EnsembleGridProviderFactory init: storage_dir={self._storage_dir}" + ) + + if self._allow_storage_writes: + os.makedirs(self._storage_dir, exist_ok=True) + + @staticmethod + def instance() -> "EnsembleGridProviderFactory": + """Static method to access the singleton instance of the factory.""" + + factory = WEBVIZ_FACTORY_REGISTRY.get_factory(EnsembleGridProviderFactory) + if not factory: + app_instance_info = WEBVIZ_FACTORY_REGISTRY.app_instance_info + storage_folder = app_instance_info.storage_folder + allow_writes = app_instance_info.run_mode != WebvizRunMode.PORTABLE + dont_copy_grid_data = ( + app_instance_info.run_mode == WebvizRunMode.NON_PORTABLE + ) + + factory = EnsembleGridProviderFactory( + root_storage_folder=storage_folder, + allow_storage_writes=allow_writes, + avoid_copying_grid_data=dont_copy_grid_data, + ) + + # Store the factory object in the global factory registry + WEBVIZ_FACTORY_REGISTRY.set_factory(EnsembleGridProviderFactory, factory) + + return factory + + def create_from_roff_files( + self, ens_path: str, grid_name: str, attribute_filter: List[str] = None + ) -> EnsembleGridProvider: + timer = PerfTimer() + string_to_hash = ( + f"{ens_path}_{grid_name}" + if attribute_filter is None + else f"{ens_path}_{grid_name}_{'_'.join([str(attr) for attr in attribute_filter])}" + ) + storage_key = f"ens__{_make_hash_string(string_to_hash)}" + provider = ProviderImplRoff.from_backing_store(self._storage_dir, storage_key) + if provider: + LOGGER.info( + f"Loaded grid provider from backing store in {timer.elapsed_s():.2f}s (" + f"ens_path={ens_path})" + ) + return provider + + # We can only import data from data source if storage writes are allowed + if not self._allow_storage_writes: + raise ValueError(f"Failed to load grid provider for {ens_path}") + + LOGGER.info(f"Importing/copying grid data for: {ens_path}") + + timer.lap_s() + grid_info, grid_parameters_info = discover_per_realization_grid_files( + ens_path, grid_name, attribute_filter + ) + + # As an optimization, avoid copying the grid data into the backing store, + # typically when we're running in non-portable mode + ProviderImplRoff.write_backing_store( + self._storage_dir, + storage_key, + grid_geometries_info=grid_info, + grid_parameters_info=grid_parameters_info, + avoid_copying_grid_data=self._avoid_copying_grid_data, + ) + et_write_s = timer.lap_s() + + provider = ProviderImplRoff.from_backing_store(self._storage_dir, storage_key) + if not provider: + raise ValueError(f"Failed to load/create grid provider for {ens_path}") + + LOGGER.info( + f"Saved grid provider to backing store in {timer.elapsed_s():.2f}s (" + f" write={et_write_s:.2f}s, ens_path={ens_path})" + ) + + return provider + + +def _make_hash_string(string_to_hash: str) -> str: + # There is no security risk here and chances of collision should be very slim + return hashlib.md5(string_to_hash.encode()).hexdigest() # nosec diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/grid_viz_service.py b/webviz_subsurface/_providers/ensemble_grid_provider/grid_viz_service.py new file mode 100644 index 000000000..a0725c728 --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/grid_viz_service.py @@ -0,0 +1,775 @@ +import dataclasses +import logging +from dataclasses import dataclass +from pathlib import Path +from typing import Any, Callable, Dict, List, Optional, Tuple + +import numpy as np +import xtgeo +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkCommonCore import reference, vtkIdList, vtkPoints +from vtkmodules.vtkCommonDataModel import ( + vtkCellArray, + vtkCellLocator, + vtkExplicitStructuredGrid, + vtkGenericCell, + vtkLine, + vtkPlane, + vtkPolyData, + vtkStaticCellLocator, + vtkUnstructuredGrid, +) +from vtkmodules.vtkFiltersCore import ( + vtkAppendPolyData, + vtkExplicitStructuredGridCrop, + vtkExplicitStructuredGridToUnstructuredGrid, + vtkPlaneCutter, + vtkClipPolyData, + vtkUnstructuredGridToExplicitStructuredGrid, +) +from vtkmodules.vtkFiltersGeneral import vtkBoxClipDataSet +from vtkmodules.vtkFiltersGeometry import vtkExplicitStructuredGridSurfaceFilter + +from webviz_subsurface._utils.perf_timer import PerfTimer + +# Requires updated xtgeo +from ._xtgeo_to_vtk_explicit_structured_grid import ( + xtgeo_grid_to_vtk_explicit_structured_grid, +) +from .ensemble_grid_provider import EnsembleGridProvider + +LOGGER = logging.getLogger(__name__) + +_GRID_VIZ_SERVICE_INSTANCE: Optional["GridVizService"] = None + + +@dataclass +class PropertySpec: + prop_name: str + prop_date: Optional[str] + + +@dataclass +class CellFilter: + i_min: int + i_max: int + j_min: int + j_max: int + k_min: int + k_max: int + + +@dataclass +class SurfacePolys: + point_arr: np.ndarray + poly_arr: np.ndarray + + +@dataclass +class PropertyScalars: + value_arr: np.ndarray + # min_value: float + # max_value: float + + +@dataclass +class Ray: + origin: List[float] + end: List[float] + # direction: List[float] + + +@dataclass +class PickResult: + cell_index: int + cell_i: int + cell_j: int + cell_k: int + intersection_point: List[float] + cell_property_value: Optional[float] + + +# ============================================================================= +class GridWorker: + # ----------------------------------------------------------------------------- + def __init__(self, full_esgrid: vtkExplicitStructuredGrid) -> None: + self._full_esgrid = full_esgrid + + self._cached_cell_filter: Optional[CellFilter] = None + self._cached_original_cell_indices: Optional[np.ndarray] = None + + # ----------------------------------------------------------------------------- + def get_full_esgrid(self) -> vtkExplicitStructuredGrid: + return self._full_esgrid + + # ----------------------------------------------------------------------------- + def get_cached_original_cell_indices( + self, cell_filter: Optional[CellFilter] + ) -> Optional[np.ndarray]: + if self._cached_original_cell_indices is None: + return None + + if cell_filter == self._cached_cell_filter: + return self._cached_original_cell_indices + + return None + + # ----------------------------------------------------------------------------- + def set_cached_original_cell_indices( + self, cell_filter: Optional[CellFilter], original_cell_indices: np.ndarray + ) -> None: + # Make copy of the cell filter + self._cached_cell_filter = dataclasses.replace(cell_filter) + self._cached_original_cell_indices = original_cell_indices + + +# ============================================================================= +class GridVizService: + # ----------------------------------------------------------------------------- + def __init__(self) -> None: + self._id_to_provider_dict: Dict[str, EnsembleGridProvider] = {} + self._key_to_worker_dict: Dict[str, GridWorker] = {} + + # ----------------------------------------------------------------------------- + @staticmethod + def instance() -> "GridVizService": + global _GRID_VIZ_SERVICE_INSTANCE + if not _GRID_VIZ_SERVICE_INSTANCE: + LOGGER.debug("Initializing GridVizService instance") + _GRID_VIZ_SERVICE_INSTANCE = GridVizService() + + return _GRID_VIZ_SERVICE_INSTANCE + + # ----------------------------------------------------------------------------- + def register_provider(self, provider: EnsembleGridProvider) -> None: + provider_id = provider.provider_id() + LOGGER.debug(f"Adding grid provider with id={provider_id}") + + existing_provider = self._id_to_provider_dict.get(provider_id) + if existing_provider: + # Issue a warning if there already is a provider registered with the same + # id AND if the actual provider instance is different. + # This wil happen until the provider factory gets caching. + if existing_provider is not provider: + LOGGER.warning( + f"Provider with id={provider_id} ignored, the id is already present" + ) + return + + self._id_to_provider_dict[provider_id] = provider + + # ----------------------------------------------------------------------------- + def get_surface( + self, + provider_id: str, + realization: int, + property_spec: Optional[PropertySpec], + cell_filter: Optional[CellFilter], + ) -> Tuple[SurfacePolys, Optional[PropertyScalars]]: + + LOGGER.debug( + f"Getting grid surface... " + f"(provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + timer = PerfTimer() + + provider = self._id_to_provider_dict.get(provider_id) + if not provider: + raise ValueError("Could not find provider") + + worker = self._get_or_create_grid_worker(provider_id, realization) + if not worker: + raise ValueError("Could not get grid worker") + + grid = worker.get_full_esgrid() + + if cell_filter: + grid = _calc_cropped_grid(grid, cell_filter) + + polydata = _calc_grid_surface(grid) + + # !!!!!! + # Need to watch out here, think these may go out of scope! + points_np = vtk_to_numpy(polydata.GetPoints().GetData()).ravel() + polys_np = vtk_to_numpy(polydata.GetPolys().GetData()) + original_cell_indices_np = vtk_to_numpy( + polydata.GetCellData().GetAbstractArray("vtkOriginalCellIds") + ) + + surface_polys = SurfacePolys(point_arr=points_np, poly_arr=polys_np) + + property_scalars: Optional[PropertyScalars] = None + if property_spec: + raw_cell_vals = _load_property_values(provider, realization, property_spec) + if raw_cell_vals is not None: + mapped_cell_vals = raw_cell_vals[original_cell_indices_np] + property_scalars = PropertyScalars(value_arr=mapped_cell_vals) + + worker.set_cached_original_cell_indices(cell_filter, original_cell_indices_np) + + LOGGER.debug( + f"Got grid surface in {timer.elapsed_s():.2f}s " + f"(provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + + return surface_polys, property_scalars + + # ----------------------------------------------------------------------------- + def get_mapped_property_values( + self, + provider_id: str, + realization: int, + property_spec: PropertySpec, + cell_filter: Optional[CellFilter], + ) -> Optional[PropertyScalars]: + + LOGGER.debug( + f"Getting property values... " + f"(provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + timer = PerfTimer() + + provider = self._id_to_provider_dict.get(provider_id) + if not provider: + raise ValueError("Could not find provider") + + worker = self._get_or_create_grid_worker(provider_id, realization) + if not worker: + raise ValueError("Could not get grid worker") + + raw_cell_vals = _load_property_values(provider, realization, property_spec) + if raw_cell_vals is None: + LOGGER.warning( + f"No cell values found for " + f"prop=({property_spec.prop_name}, {property_spec.prop_name})" + ) + return None + + original_cell_indices_np = worker.get_cached_original_cell_indices(cell_filter) + if original_cell_indices_np is None: + # Must first generate the grid to get the original cell indices + grid = worker.get_full_esgrid() + if cell_filter: + grid = _calc_cropped_grid(grid, cell_filter) + + polydata = _calc_grid_surface(grid) + original_cell_indices_np = vtk_to_numpy( + polydata.GetCellData().GetAbstractArray("vtkOriginalCellIds") + ) + worker.set_cached_original_cell_indices( + cell_filter, original_cell_indices_np + ) + + mapped_cell_vals = raw_cell_vals[original_cell_indices_np] + + LOGGER.debug( + f"Got property values in {timer.elapsed_s():.2f}s " + f"(provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + + return PropertyScalars(value_arr=mapped_cell_vals) + + # ----------------------------------------------------------------------------- + def cut_along_polyline( + self, + provider_id: str, + realization: int, + polyline_xy: List[float], + property_spec: Optional[PropertySpec], + ) -> Tuple[SurfacePolys, Optional[PropertyScalars]]: + + LOGGER.debug( + f"Cutting along polyline... " + f"(provider_id={provider_id}, real={realization})" + ) + timer = PerfTimer() + + provider = self._id_to_provider_dict.get(provider_id) + if not provider: + raise ValueError("Could not find provider") + + worker = self._get_or_create_grid_worker(provider_id, realization) + if not worker: + raise ValueError("Could not get grid worker") + + esgrid = worker.get_full_esgrid() + + bounds = esgrid.GetBounds() + min_z = bounds[4] + max_z = bounds[5] + + num_points_in_polyline = int(len(polyline_xy) / 2) + + ugrid = _vtk_esg_to_ug(esgrid) + + # !!!!!!!!!!!!!! + # Requires VTK 9.2-ish + # ugrid = _extract_intersected_ugrid(ugrid, polyline_xy, 10.0) + + cutter_alg = vtkPlaneCutter() + cutter_alg.SetInputDataObject(ugrid) + + # cell_locator = vtkStaticCellLocator() + # cell_locator.SetDataSet(esgrid) + # cell_locator.BuildLocator() + + # box_clip_alg = vtkBoxClipDataSet() + # box_clip_alg.SetInputDataObject(ugrid) + + append_alg = vtkAppendPolyData() + et_setup_s = timer.lap_s() + + et_cut_s = 0.0 + et_clip_s = 0.0 + + for i in range(0, num_points_in_polyline - 1): + x0 = polyline_xy[2 * i] + y0 = polyline_xy[2 * i + 1] + x1 = polyline_xy[2 * (i + 1)] + y1 = polyline_xy[2 * (i + 1) + 1] + fwd_vec = np.array([x1 - x0, y1 - y0, 0.0]) + fwd_vec /= np.linalg.norm(fwd_vec) + right_vec = np.array([fwd_vec[1], -fwd_vec[0], 0]) + + # box_clip_alg.SetBoxClip(x0, x1, y0, y1, min_z, max_z) + # box_clip_alg.Update() + # clipped_ugrid = box_clip_alg.GetOutputDataObject(0) + + # polyline_bounds = _calc_polyline_bounds([x0, y0, x1, y1]) + # polyline_bounds.extend([min_z, max_z]) + # cell_ids = vtkIdList() + # cell_locator.FindCellsWithinBounds(polyline_bounds, cell_ids) + # print(f"{cell_ids.GetNumberOfIds()} {polyline_bounds=}") + + plane = vtkPlane() + plane.SetOrigin([x0, y0, 0]) + plane.SetNormal(right_vec) + + plane_0 = vtkPlane() + plane_0.SetOrigin([x0, y0, 0]) + plane_0.SetNormal(fwd_vec) + + plane_1 = vtkPlane() + plane_1.SetOrigin([x1, y1, 0]) + plane_1.SetNormal(-fwd_vec) + + cutter_alg.SetPlane(plane) + cutter_alg.Update() + + cut_surface_polydata = cutter_alg.GetOutput() + # print(f"{type(cut_surface_polydata)=}") + et_cut_s += timer.lap_s() + + # Used vtkPolyDataPlaneClipper earlier, but it seems that it doesn't + # maintain the original cell IDs that we need for the result mapping. + # May want to check up on any performance degradation! + clipper_0 = vtkClipPolyData() + clipper_0.SetInputDataObject(cut_surface_polydata) + clipper_0.SetClipFunction(plane_0) + clipper_0.Update() + clipped_polydata = clipper_0.GetOutput() + + clipper_1 = vtkClipPolyData() + clipper_1.SetInputDataObject(clipped_polydata) + clipper_1.SetClipFunction(plane_1) + clipper_1.Update() + clipped_polydata = clipper_1.GetOutput() + + append_alg.AddInputData(clipped_polydata) + + et_clip_s += timer.lap_s() + + append_alg.Update() + comb_polydata = append_alg.GetOutput() + et_combine_s = timer.lap_s() + + points_np = vtk_to_numpy(comb_polydata.GetPoints().GetData()).ravel() + polys_np = vtk_to_numpy(comb_polydata.GetPolys().GetData()) + + surface_polys = SurfacePolys(point_arr=points_np, poly_arr=polys_np) + + property_scalars: Optional[PropertyScalars] = None + if property_spec: + raw_cell_vals = _load_property_values(provider, realization, property_spec) + if raw_cell_vals is not None: + original_cell_indices_np = vtk_to_numpy( + comb_polydata.GetCellData().GetAbstractArray("vtkOriginalCellIds") + ) + mapped_cell_vals = raw_cell_vals[original_cell_indices_np] + property_scalars = PropertyScalars(value_arr=mapped_cell_vals) + + LOGGER.debug( + f"Cutting along polyline done in {timer.elapsed_s():.2f}s " + f"setup={et_setup_s:.2f}s, cut={et_cut_s:.2f}s, clip={et_clip_s:.2f}s, combine={et_combine_s:.2f}s, " + f"(provider_id={provider_id}, real={realization})" + ) + + return surface_polys, property_scalars + + """ + dbg_point_arr = [] + dbg_conn_arr = [] + for i in range(0, num_points_in_polyline): + x = polyline_xy[2 * i] + y = polyline_xy[2 * i + 1] + dbg_point_arr.extend([x, y, min_z]) + dbg_point_arr.extend([x, y, max_z]) + if i > 0: + base = 2 * (i - 1) + dbg_conn_arr.extend([4, base, base + 2, base + 3, base + 1]) + + for i in range(0, int(len(dbg_point_arr) / 3)): + print( + f"{i}: {dbg_point_arr[3*i]}, {dbg_point_arr[3*i + 1]}, {dbg_point_arr[3*i + 2]}" + ) + + point_arr_np = np.array(dbg_point_arr).reshape(-1, 3) + conn_arr_np = np.array(dbg_conn_arr) + surface_polys = SurfacePolys(point_arr=point_arr_np, poly_arr=conn_arr_np) + + return surface_polys, None + """ + + # ----------------------------------------------------------------------------- + def ray_pick( + self, + provider_id: str, + realization: int, + ray: Ray, + property_spec: Optional[PropertySpec], + cell_filter: Optional[CellFilter], + ) -> Optional[PickResult]: + + LOGGER.debug( + f"Doing ray pick: " + f"ray.origin={ray.origin}, ray.end={ray.end}, " + f"(provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + timer = PerfTimer() + + provider = self._id_to_provider_dict.get(provider_id) + if not provider: + raise ValueError("Could not find provider") + + worker = self._get_or_create_grid_worker(provider_id, realization) + if not worker: + raise ValueError("Could not get grid worker") + + grid = worker.get_full_esgrid() + if cell_filter: + grid = _calc_cropped_grid(grid, cell_filter) + et_crop_s = timer.lap_s() + + cell_id, isect_pt = _raypick_in_grid(grid, ray) + et_pick_s = timer.lap_s() + if cell_id is None: + return None + + original_cell_id = cell_id + if cell_filter: + # If a cell filter is present, assume picking was done against cropped grid + original_cell_id = ( + grid.GetCellData() + .GetAbstractArray("vtkOriginalCellIds") + .GetValue(cell_id) + ) + + i_ref = reference(0) + j_ref = reference(0) + k_ref = reference(0) + grid.ComputeCellStructuredCoords(cell_id, i_ref, j_ref, k_ref, True) + + cell_property_val: Optional[float] = None + if property_spec: + raw_cell_vals = _load_property_values(provider, realization, property_spec) + if raw_cell_vals is not None: + cell_property_val = raw_cell_vals[original_cell_id] + et_props_s = timer.lap_s() + + LOGGER.debug( + f"Did ray pick in {timer.elapsed_s():.2f}s (" + f"crop={et_crop_s:.2f}s, pick={et_pick_s:.2f}s, props={et_props_s:.2f}s, " + f"provider_id={provider_id}, real={realization}, " + f"{_property_spec_dbg_str(property_spec)}, " + f"{_cell_filter_dbg_str(cell_filter)})" + ) + + return PickResult( + cell_index=original_cell_id, + cell_i=i_ref.get(), + cell_j=j_ref.get(), + cell_k=k_ref.get(), + intersection_point=isect_pt, + cell_property_value=cell_property_val, + ) + + # ----------------------------------------------------------------------------- + def _get_or_create_grid_worker( + self, provider_id: str, realization: int + ) -> Optional[GridWorker]: + + worker_key = f"P{provider_id}__R{realization}" + worker = self._key_to_worker_dict.get(worker_key) + if worker: + return worker + + provider = self._id_to_provider_dict.get(provider_id) + if not provider: + raise ValueError("Could not find provider") + + xtg_grid = provider.get_3dgrid(realization=realization) + vtk_esg = xtgeo_grid_to_vtk_explicit_structured_grid(xtg_grid) + + worker = GridWorker(vtk_esg) + self._key_to_worker_dict[worker_key] = worker + + return worker + + +# ----------------------------------------------------------------------------- +def _calc_cropped_grid( + esgrid: vtkExplicitStructuredGrid, cell_filter: CellFilter +) -> vtkExplicitStructuredGrid: + crop_filter = vtkExplicitStructuredGridCrop() + crop_filter.SetInputData(esgrid) + + # In VTK dimensions correspond to points + crop_filter.SetOutputWholeExtent( + cell_filter.i_min, + cell_filter.i_max + 1, + cell_filter.j_min, + cell_filter.j_max + 1, + cell_filter.k_min, + cell_filter.k_max + 1, + ) + crop_filter.Update() + cropped_grid = crop_filter.GetOutput() + return cropped_grid + + +# ----------------------------------------------------------------------------- +def _raypick_in_grid( + esgrid: vtkExplicitStructuredGrid, ray: Ray +) -> Optional[Tuple[int, List[float]]]: + """Do a ray pick against the specified grid. + Returns None if nothing was hit, otherwise returns the cellId (cell index) of the cell + that was hit and the intersection point + """ + + locator = vtkCellLocator() + locator.SetDataSet(esgrid) + locator.BuildLocator() + + tolerance = 0.0 + t_ref = reference(0.0) + isect_pt = [0.0, 0.0, 0.0] + pcoords = [0.0, 0.0, 0.0] + sub_id_ref = reference(0) + cell_id_ref = reference(0) + cell = vtkGenericCell() + + # From doc for vtkCell it seems that isect_pt will be the actual intersection point + # while pcoords is in parametric coordinates + anyHits = locator.IntersectWithLine( + ray.origin, + ray.end, + tolerance, + t_ref, + isect_pt, + pcoords, + sub_id_ref, + cell_id_ref, + cell, + ) + + if not anyHits: + return None + + return cell_id_ref.get(), isect_pt + + +# ----------------------------------------------------------------------------- +def _calc_grid_surface(esgrid: vtkExplicitStructuredGrid) -> vtkPolyData: + surf_filter = vtkExplicitStructuredGridSurfaceFilter() + surf_filter.SetInputData(esgrid) + surf_filter.PassThroughCellIdsOn() + surf_filter.Update() + + polydata: vtkPolyData = surf_filter.GetOutput() + return polydata + + +# ----------------------------------------------------------------------------- +def _load_property_values( + provider: EnsembleGridProvider, realization: int, property_spec: PropertySpec +) -> Optional[np.ndarray]: + if property_spec.prop_date: + prop_values = provider.get_dynamic_property_values( + property_spec.prop_name, property_spec.prop_date, realization + ) + else: + prop_values = provider.get_static_property_values( + property_spec.prop_name, realization + ) + + return prop_values + + +# ----------------------------------------------------------------------------- +def _vtk_esg_to_ug(vtk_esgrid: vtkExplicitStructuredGrid) -> vtkUnstructuredGrid: + convertFilter = vtkExplicitStructuredGridToUnstructuredGrid() + convertFilter.SetInputData(vtk_esgrid) + convertFilter.Update() + vtk_ugrid = convertFilter.GetOutput() + + return vtk_ugrid + + +# ----------------------------------------------------------------------------- +def _vtk_ug_to_esg(vtk_ugrid: vtkUnstructuredGrid) -> vtkExplicitStructuredGrid: + convertFilter = vtkUnstructuredGridToExplicitStructuredGrid() + convertFilter.SetInputData(vtk_ugrid) + convertFilter.SetInputArrayToProcess(0, 0, 0, 1, "BLOCK_I") + convertFilter.SetInputArrayToProcess(1, 0, 0, 1, "BLOCK_J") + convertFilter.SetInputArrayToProcess(2, 0, 0, 1, "BLOCK_K") + convertFilter.Update() + vtk_esgrid = convertFilter.GetOutput() + + return vtk_esgrid + + +# ----------------------------------------------------------------------------- +def _calc_polyline_bounds(polyline_xy: List[float]) -> List[float]: + num_points = int(len(polyline_xy) / 2) + if num_points < 1: + return None + + min_x = min(polyline_xy[0::2]) + max_x = max(polyline_xy[0::2]) + min_y = min(polyline_xy[1::2]) + max_y = max(polyline_xy[1::2]) + + return [min_x, max_x, min_y, max_y] + + +# ----------------------------------------------------------------------------- +def _extract_intersected_ugrid( + ugrid: vtkUnstructuredGrid, polyline_xy_in: List[float], max_point_dist: float +) -> vtkUnstructuredGrid: + + # Requires VTK 9.2 + from vtkmodules.vtkFiltersCore import vtkExtractCellsAlongPolyLine + + timer = PerfTimer() + + polyline_xy = _resample_polyline(polyline_xy_in, max_point_dist) + et_resample_s = timer.lap_s() + + num_points_in_polyline = int(len(polyline_xy) / 2) + if num_points_in_polyline < 1: + return ugrid + + bounds = ugrid.GetBounds() + min_z = bounds[4] + max_z = bounds[5] + + points = vtkPoints() + lines = vtkCellArray() + + for i in range(0, num_points_in_polyline): + x = polyline_xy[2 * i] + y = polyline_xy[2 * i + 1] + + points.InsertNextPoint([x, y, min_z]) + points.InsertNextPoint([x, y, max_z]) + + line = vtkLine() + line.GetPointIds().SetId(0, 2 * i) + line.GetPointIds().SetId(1, 2 * i + 1) + + lines.InsertNextCell(line) + + polyData = vtkPolyData() + polyData.SetPoints(points) + polyData.SetLines(lines) + + et_build_s = timer.lap_s() + + extractor = vtkExtractCellsAlongPolyLine() + extractor.SetInputData(0, ugrid) + extractor.SetInputData(1, polyData) + extractor.Update() + + ret_grid = extractor.GetOutput(0) + + et_extract_s = timer.lap_s() + + LOGGER.debug( + f"extraction with {num_points_in_polyline} points took {timer.elapsed_s():.2f}s " + f"(resample={et_resample_s:.2f}s, build={et_build_s:.2f}s, extract={et_extract_s:.2f}s)" + ) + + return ret_grid + + +# ----------------------------------------------------------------------------- +def _resample_polyline(polyline_xy: List[float], max_point_dist: float) -> List[float]: + num_points = int(len(polyline_xy) / 2) + if num_points < 2: + return polyline_xy + + ret_polyline = [] + + prev_x = polyline_xy[0] + prev_y = polyline_xy[1] + ret_polyline.extend([prev_x, prev_y]) + + for i in range(1, num_points): + x = polyline_xy[2 * i] + y = polyline_xy[2 * i + 1] + + fwd = [x - prev_x, y - prev_y] + length = np.linalg.norm(fwd) + if length > max_point_dist: + n = int(length / max_point_dist) + delta_t = 1.0 / (n + 1) + for j in range(0, n): + pt_x = prev_x + fwd[0] * (j + 1) * delta_t + pt_y = prev_y + fwd[1] * (j + 1) * delta_t + ret_polyline.extend([pt_x, pt_y]) + + ret_polyline.extend([x, y]) + prev_x = x + prev_y = y + + return ret_polyline + + +# ----------------------------------------------------------------------------- +def _property_spec_dbg_str(property_spec: Optional[PropertySpec]) -> str: + if not property_spec: + return "prop=None" + + return f"prop=({property_spec.prop_name}, {property_spec.prop_date})" + + +# ----------------------------------------------------------------------------- +def _cell_filter_dbg_str(cell_filter: Optional[CellFilter]) -> str: + if not cell_filter: + return "IJK=None" + + return ( + f"I=[{cell_filter.i_min},{cell_filter.i_max}] " + f"J=[{cell_filter.j_min},{cell_filter.j_max}] " + f"K=[{cell_filter.k_min},{cell_filter.k_max}]" + ) diff --git a/webviz_subsurface/_providers/ensemble_grid_provider/provider_impl_roff.py b/webviz_subsurface/_providers/ensemble_grid_provider/provider_impl_roff.py new file mode 100644 index 000000000..a426a07fa --- /dev/null +++ b/webviz_subsurface/_providers/ensemble_grid_provider/provider_impl_roff.py @@ -0,0 +1,356 @@ +import logging +import shutil +import warnings +from enum import Enum +from pathlib import Path +from typing import List, Optional, Set + +import numpy as np +import pandas as pd +import xtgeo + +from webviz_subsurface._utils.perf_timer import PerfTimer + + +from ._roff_file_discovery import GridFileInfo, GridParameterFileInfo +from .ensemble_grid_provider import EnsembleGridProvider + + +LOGGER = logging.getLogger(__name__) + + +# pylint: disable=too-few-public-methods +class Col: + TYPE = "type" + REAL = "real" + ATTRIBUTE = "attribute" + NAME = "name" + DATESTR = "datestr" + ORIGINAL_PATH = "original_path" + REL_PATH = "rel_path" + + +class GridType(str, Enum): + GEOMETRY = "geometry" + STATIC_PROPERTY = "static_property" + DYNAMIC_PROPERTY = "dynamic_property" + + +class ProviderImplRoff(EnsembleGridProvider): + def __init__( + self, provider_id: str, provider_dir: Path, grid_inventory_df: pd.DataFrame + ) -> None: + self._provider_id = provider_id + self._provider_dir = provider_dir + self._inventory_df = grid_inventory_df + + @staticmethod + # pylint: disable=too-many-locals + def write_backing_store( + storage_dir: Path, + storage_key: str, + grid_geometries_info: List[GridFileInfo], + grid_parameters_info: List[GridParameterFileInfo], + avoid_copying_grid_data: bool, + ) -> None: + """If avoid_copying_grid_data if True, the specified grid data will NOT be copied + into the backing store, but will be referenced from their source locations. + Note that this is only useful when running in non-portable mode and will fail + in portable mode. + """ + + timer = PerfTimer() + + do_copy_grid_data_into_store = not avoid_copying_grid_data + + # All data for this provider will be stored inside a sub-directory + # given by the storage key + provider_dir = storage_dir / storage_key + LOGGER.debug(f"Writing grid data backing store to: {provider_dir}") + provider_dir.mkdir(parents=True, exist_ok=True) + + type_arr: List[GridType] = [] + real_arr: List[int] = [] + attribute_arr: List[str] = [] + name_arr: List[str] = [] + datestr_arr: List[str] = [] + rel_path_arr: List[str] = [] + original_path_arr: List[str] = [] + for grid_info in grid_geometries_info: + type_arr.append(GridType.GEOMETRY) + name_arr.append(grid_info.name) + real_arr.append(grid_info.real) + attribute_arr.append("") + datestr_arr.append("") + original_path_arr.append(grid_info.path) + rel_path_in_store = "" + + if do_copy_grid_data_into_store: + rel_path_in_store = _compose_rel_grid_pathstr( + real=grid_info.real, + attribute=None, + name=grid_info.name, + datestr=None, + extension=Path(grid_info.path).suffix, + ) + + rel_path_arr.append(rel_path_in_store) + + for grid_parameter_info in grid_parameters_info: + name_arr.append(grid_parameter_info.name) + real_arr.append(grid_parameter_info.real) + attribute_arr.append(grid_parameter_info.attribute) + if grid_parameter_info.datestr: + datestr_arr.append(grid_parameter_info.datestr) + type_arr.append(GridType.DYNAMIC_PROPERTY) + else: + datestr_arr.append("") + type_arr.append(GridType.STATIC_PROPERTY) + + original_path_arr.append(grid_parameter_info.path) + + rel_path_in_store = "" + if do_copy_grid_data_into_store: + rel_path_in_store = _compose_rel_grid_pathstr( + real=grid_parameter_info.real, + attribute=grid_parameter_info.attribute, + name=grid_parameter_info.name, + datestr=grid_parameter_info.datestr, + extension=Path(grid_parameter_info.path).suffix, + ) + + rel_path_arr.append(rel_path_in_store) + + timer.lap_s() + if do_copy_grid_data_into_store: + LOGGER.debug( + f"Copying {len(original_path_arr)} grid data into backing store..." + ) + _copy_grid_parameters_into_provider_dir( + original_path_arr, rel_path_arr, provider_dir + ) + et_copy_s = timer.lap_s() + + grid_inventory_df = pd.DataFrame( + { + Col.TYPE: type_arr, + Col.REAL: real_arr, + Col.ATTRIBUTE: attribute_arr, + Col.NAME: name_arr, + Col.DATESTR: datestr_arr, + Col.REL_PATH: rel_path_arr, + Col.ORIGINAL_PATH: original_path_arr, + } + ) + + parquet_file_name = provider_dir / "grid_inventory.parquet" + grid_inventory_df.to_parquet(path=parquet_file_name) + + if do_copy_grid_data_into_store: + LOGGER.debug( + f"Wrote grid backing store in: {timer.elapsed_s():.2f}s (" + f"copy={et_copy_s:.2f}s)" + ) + else: + LOGGER.debug( + f"Wrote grid backing store without copying grid data in: " + f"{timer.elapsed_s():.2f}s" + ) + + @staticmethod + def from_backing_store( + storage_dir: Path, + storage_key: str, + ) -> Optional["ProviderImplRoff"]: + + provider_dir = storage_dir / storage_key + parquet_file_name = provider_dir / "grid_inventory.parquet" + + try: + grid_inventory_df = pd.read_parquet(path=parquet_file_name) + return ProviderImplRoff(storage_key, provider_dir, grid_inventory_df) + except FileNotFoundError: + return None + + def provider_id(self) -> str: + return self._provider_id + + def static_property_names(self) -> List[str]: + return sorted( + list( + self._inventory_df.loc[ + self._inventory_df[Col.TYPE] == GridType.STATIC_PROPERTY + ][Col.ATTRIBUTE].unique() + ) + ) + + def dynamic_property_names(self) -> List[str]: + return sorted( + list( + self._inventory_df.loc[ + self._inventory_df[Col.TYPE] == GridType.DYNAMIC_PROPERTY + ][Col.ATTRIBUTE].unique() + ) + ) + + def dates_for_dynamic_property(self, property_name: str) -> Optional[List[str]]: + print(property_name) + dates = sorted( + list( + self._inventory_df.loc[ + (self._inventory_df[Col.TYPE] == GridType.DYNAMIC_PROPERTY) + & (self._inventory_df[Col.ATTRIBUTE] == property_name) + ][Col.DATESTR].unique() + ) + ) + if len(dates) == 1 and not bool(dates[0]): + return None + + return dates + + def realizations(self) -> List[int]: + unique_reals = self._inventory_df[Col.REAL].unique() + + # Sort and strip out any entries with real == -1 + return sorted([r for r in unique_reals if r >= 0]) + + def get_3dgrid(self, realization: int) -> xtgeo.Grid: + df = self._inventory_df.loc[self._inventory_df[Col.TYPE] == GridType.GEOMETRY] + print(df) + df = df.loc[df[Col.REAL] == realization] + + df = df[[Col.REL_PATH, Col.ORIGINAL_PATH]] + fn_list: List[str] = [] + for _index, row in df.iterrows(): + if row[Col.REL_PATH]: + fn_list.append(self._provider_dir / row[Col.REL_PATH]) + else: + fn_list.append(row[Col.ORIGINAL_PATH]) + if len(fn_list) == 0: + LOGGER.warning(f"No grid geometry found for realization {realization}") + return None + if len(fn_list) > 1: + raise ValueError( + f"Multiple grid geometries found for: {realization}" + "Something has gone terribly wrong." + ) + + grid = xtgeo.grid_from_file(fn_list[0]) + return grid + + def get_static_property_values( + self, property_name: str, realization: int + ) -> Optional[np.ndarray]: + fn_list: List[str] = self._locate_static_property( + property_name=property_name, realizations=[realization] + ) + if len(fn_list) == 0: + LOGGER.warning(f"No grid parameter found for realization {realization}") + return None + if len(fn_list) > 1: + raise ValueError( + f"Multiple grid parameters found for: {realization}" + "Something has gone terribly wrong." + ) + grid_property = xtgeo.gridproperty_from_file(fn_list[0]) + fill_value = np.nan if not grid_property.isdiscrete else -1 + return grid_property.get_npvalues1d(order="F", fill_value=fill_value).ravel() + + def get_dynamic_property_values( + self, property_name: str, property_date: str, realization: int + ) -> Optional[np.ndarray]: + fn_list: List[str] = self._locate_dynamic_property( + property_name=property_name, + property_datestr=property_date, + realizations=[realization], + ) + if len(fn_list) == 0: + LOGGER.warning(f"No grid parameter found for realization {realization}") + return None + if len(fn_list) > 1: + raise ValueError( + f"Multiple grid parameters found for: {realization}" + "Something has gone terribly wrong." + ) + grid_property = xtgeo.gridproperty_from_file(fn_list[0]) + return grid_property.get_npvalues1d(order="F").ravel() + + def _locate_static_property( + self, property_name: str, realizations: List[int] + ) -> List[str]: + """Returns list of file names matching the specified filter criteria""" + df = self._inventory_df.loc[ + self._inventory_df[Col.TYPE] == GridType.STATIC_PROPERTY + ] + + df = df.loc[ + (df[Col.ATTRIBUTE] == property_name) & (df[Col.REAL].isin(realizations)) + ] + + df = df[[Col.REL_PATH, Col.ORIGINAL_PATH]] + + # Return file name within backing store if the data was copied there, + # otherwise return the original source file name + fn_list: List[str] = [] + for _index, row in df.iterrows(): + if row[Col.REL_PATH]: + fn_list.append(self._provider_dir / row[Col.REL_PATH]) + else: + fn_list.append(row[Col.ORIGINAL_PATH]) + + return fn_list + + def _locate_dynamic_property( + self, property_name: str, property_datestr: str, realizations: List[int] + ) -> List[str]: + """Returns list of file names matching the specified filter criteria""" + df = self._inventory_df.loc[ + self._inventory_df[Col.TYPE] == GridType.DYNAMIC_PROPERTY + ] + + df = df.loc[ + (df[Col.ATTRIBUTE] == property_name) + & (df[Col.DATESTR] == property_datestr) + & (df[Col.REAL].isin(realizations)) + ] + + df = df[[Col.REL_PATH, Col.ORIGINAL_PATH]] + + # Return file name within backing store if the data was copied there, + # otherwise return the original source file name + fn_list: List[str] = [] + for _index, row in df.iterrows(): + if row[Col.REL_PATH]: + fn_list.append(self._provider_dir / row[Col.REL_PATH]) + else: + fn_list.append(row[Col.ORIGINAL_PATH]) + + return fn_list + + +def _copy_grid_parameters_into_provider_dir( + original_path_arr: List[str], + rel_path_arr: List[str], + provider_dir: Path, +) -> None: + for src_path, dst_rel_path in zip(original_path_arr, rel_path_arr): + shutil.copyfile(src_path, provider_dir / dst_rel_path) + + # full_dst_path_arr = [storage_dir / dst_rel_path for dst_rel_path in store_path_arr] + # with ProcessPoolExecutor() as executor: + # executor.map(shutil.copyfile, original_path_arr, full_dst_path_arr) + + +def _compose_rel_grid_pathstr( + real: int, + name: str, + attribute: Optional[str], + datestr: Optional[str], + extension: str, +) -> str: + """Compose path to grid file, relative to provider's directory""" + if not attribute and not datestr: + return str(Path(f"{real}--{name}{extension}")) + if not datestr: + return str(Path(f"{real}--{name}--{attribute}{extension}")) + return str(Path(f"{real}--{name}--{attribute}--{datestr}{extension}")) diff --git a/webviz_subsurface/_providers/well_provider/_provider_impl_file.py b/webviz_subsurface/_providers/well_provider/_provider_impl_file.py index f1f006ebe..62b61bf8a 100644 --- a/webviz_subsurface/_providers/well_provider/_provider_impl_file.py +++ b/webviz_subsurface/_providers/well_provider/_provider_impl_file.py @@ -3,11 +3,14 @@ from pathlib import Path from typing import Dict, List, Optional +import numpy as np +import pandas as pd import xtgeo from webviz_subsurface._utils.perf_timer import PerfTimer -from .well_provider import WellPath, WellProvider +from .well_provider import WellPath, WellProvider, WellIntersectionPolyLine +from ._simplify_polyline import rdp LOGGER = logging.getLogger(__name__) @@ -131,3 +134,40 @@ def get_well_xtgeo_obj(self, well_name: str) -> xtgeo.Well: ) return well + + def get_polyline_along_well_path_SIMPLIFIED( + self, well_name: str, tvdmin=None, use_rdp=True + ) -> np.array: + """Returns a polyline for the well path along with MD for the well.""" + well = self.get_well_xtgeo_obj(well_name).copy() + if tvdmin is not None: + well.dataframe = well.dataframe[well.dataframe["Z_TVDSS"] >= tvdmin] + + xy_arr = np.stack( + ( + np.array(well.dataframe["X_UTME"].values), + np.array(well.dataframe["Y_UTMN"].values), + ), + axis=-1, + ) + + timer = PerfTimer() + if np.all(xy_arr == xy_arr[0]): + xy_start = xy_arr[0] - 500 + xy_end = xy_arr[0] + 500 + + print( + f"Well is vertical. Returning two points extended in xy from trajectory" + ) + return [xy_start, xy_end] + if not use_rdp: + return xy_arr + simplified_xy_arr = rdp(xy_arr) + + print( + f"Original well polyline has {len(xy_arr)} points. " + f"Simplified well polyline has {len(simplified_xy_arr)} points. ", + f"Time to calculate: {timer.elapsed_ms()}ms", + ) + + return simplified_xy_arr diff --git a/webviz_subsurface/_providers/well_provider/_simplify_polyline.py b/webviz_subsurface/_providers/well_provider/_simplify_polyline.py new file mode 100644 index 000000000..de4d9b6bf --- /dev/null +++ b/webviz_subsurface/_providers/well_provider/_simplify_polyline.py @@ -0,0 +1,27 @@ +import numpy as np + + +def rdp(points, epsilon=5): + # https://towardsdatascience.com/simplify-polylines-with-the-douglas-peucker-algorithm-ac8ed487a4a1 + # get the start and end points + + start = np.tile(np.expand_dims(points[0], axis=0), (points.shape[0], 1)) + end = np.tile(np.expand_dims(points[-1], axis=0), (points.shape[0], 1)) + # find distance from other_points to line formed by start and end + dist_point_to_line = np.abs( + np.cross(end - start, points - start, axis=-1) + ) / np.linalg.norm(end - start, axis=-1) + # get the index of the points with the largest distance + max_idx = np.argmax(dist_point_to_line) + max_value = dist_point_to_line[max_idx] + + result = [] + if max_value > epsilon: + partial_results_left = rdp(points[: max_idx + 1], epsilon) + result += [list(i) for i in partial_results_left if list(i) not in result] + partial_results_right = rdp(points[max_idx:], epsilon) + result += [list(i) for i in partial_results_right if list(i) not in result] + else: + result += [points[0], points[-1]] + + return result diff --git a/webviz_subsurface/_providers/well_provider/well_provider.py b/webviz_subsurface/_providers/well_provider/well_provider.py index 9ddb60bda..ebd8c6132 100644 --- a/webviz_subsurface/_providers/well_provider/well_provider.py +++ b/webviz_subsurface/_providers/well_provider/well_provider.py @@ -14,6 +14,12 @@ class WellPath: md_arr: np.ndarray +@dataclass(frozen=True) +class WellIntersectionPolyLine: + x_arr: np.ndarray + y_arr: np.ndarray + + # Class provides data for wells class WellProvider(abc.ABC): @abc.abstractmethod @@ -31,3 +37,10 @@ def get_well_path(self, well_name: str) -> WellPath: @abc.abstractmethod def get_well_xtgeo_obj(self, well_name: str) -> xtgeo.Well: ... + + @abc.abstractmethod + def get_polyline_along_well_path_SIMPLIFIED( + self, well_name: str, tvdmin=None, use_rdp: bool = True + ) -> WellIntersectionPolyLine: + """Returns a polyline for the well path.""" + ... diff --git a/webviz_subsurface/_providers/well_provider/well_server.py b/webviz_subsurface/_providers/well_provider/well_server.py index 7063855cb..24eb06b27 100644 --- a/webviz_subsurface/_providers/well_provider/well_server.py +++ b/webviz_subsurface/_providers/well_provider/well_server.py @@ -1,10 +1,19 @@ import logging from typing import Dict, List, Optional from urllib.parse import quote +from dataclasses import dataclass import flask import geojson +import numpy as np from dash import Dash +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonDataModel import ( + vtkCellArray, + vtkPolyData, + vtkPolyLine, +) +from vtkmodules.util.numpy_support import vtk_to_numpy from webviz_subsurface._providers.well_provider.well_provider import WellProvider from webviz_subsurface._utils.perf_timer import PerfTimer @@ -16,6 +25,12 @@ _WELL_SERVER_INSTANCE: Optional["WellServer"] = None +@dataclass +class PolyLine: + point_arr: np.ndarray + line_arr: np.ndarray + + class WellServer: def __init__(self, app: Dash) -> None: self._setup_url_rule(app) @@ -31,7 +46,7 @@ def instance(app: Dash) -> "WellServer": return _WELL_SERVER_INSTANCE - def add_provider(self, provider: WellProvider) -> None: + def register_provider(self, provider: WellProvider) -> None: provider_id = provider.provider_id() LOGGER.debug(f"Adding provider with id={provider_id}") @@ -113,3 +128,36 @@ def _handle_wells_request( LOGGER.debug(f"Request handled in: {timer.elapsed_s():.2f}s") return response + + def get_polyline( + self, provider_id: str, well_name: str, tvdmin: float = None + ) -> PolyLine: + provider = self._id_to_provider_dict[provider_id] + well_path = provider.get_well_path(well_name) + xyz_arr = [ + [x, y, z] + for x, y, z in zip(well_path.x_arr, well_path.y_arr, well_path.z_arr) + ] + points = vtkPoints() + for p in xyz_arr: + points.InsertNextPoint(p[0], p[1], -p[2]) + + polyLine = vtkPolyLine() + polyLine.GetPointIds().SetNumberOfIds(len(xyz_arr)) + for i in range(0, len(xyz_arr)): + polyLine.GetPointIds().SetId(i, i) + + # Create a cell array to store the lines in and add the lines to it + cells = vtkCellArray() + cells.InsertNextCell(polyLine) + + # Create a polydata to store everything in + polyData = vtkPolyData() + # Add the points to the dataset + polyData.SetPoints(points) + + # Add the lines to the dataset + polyData.SetLines(cells) + points = vtk_to_numpy(polyData.GetPoints().GetData()) + lines = vtk_to_numpy(polyData.GetLines().GetData()) + return PolyLine(point_arr=points, line_arr=lines) diff --git a/webviz_subsurface/plugins/__init__.py b/webviz_subsurface/plugins/__init__.py index 27a8548ac..c90cf1dcb 100644 --- a/webviz_subsurface/plugins/__init__.py +++ b/webviz_subsurface/plugins/__init__.py @@ -23,6 +23,8 @@ from ._assisted_history_matching_analysis import AssistedHistoryMatchingAnalysis from ._bhp_qc import BhpQc from ._disk_usage import DiskUsage +from ._eclipse_grid_viewer import EclipseGridViewer +from ._grid_viewer import GridViewer from ._group_tree import GroupTree from ._history_match import HistoryMatch from ._horizon_uncertainty_viewer import HorizonUncertaintyViewer diff --git a/webviz_subsurface/plugins/_eclipse_grid_viewer/__init__.py b/webviz_subsurface/plugins/_eclipse_grid_viewer/__init__.py new file mode 100644 index 000000000..6a0678213 --- /dev/null +++ b/webviz_subsurface/plugins/_eclipse_grid_viewer/__init__.py @@ -0,0 +1 @@ +from ._plugin import EclipseGridViewer diff --git a/webviz_subsurface/plugins/_eclipse_grid_viewer/_callbacks.py b/webviz_subsurface/plugins/_eclipse_grid_viewer/_callbacks.py new file mode 100644 index 000000000..780926a36 --- /dev/null +++ b/webviz_subsurface/plugins/_eclipse_grid_viewer/_callbacks.py @@ -0,0 +1,513 @@ +import hashlib +import json +from time import time +from typing import Any, Callable, Dict, List, Optional, Tuple + +import numpy as np +from dash import Input, Output, State, callback, no_update, callback_context, MATCH, ALL +from webviz_vtk.utils.vtk import b64_encode_numpy + +from webviz_subsurface._utils.perf_timer import PerfTimer +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProvider, + GridVizService, + PropertySpec, + CellFilter, + Ray, +) +from webviz_subsurface._providers.well_provider import WellProvider, WellServer + +from ._layout import PROPERTYTYPE, LayoutElements, GRID_DIRECTION + + +def plugin_callbacks( + get_uuid: Callable, + grid_provider: EnsembleGridProvider, + grid_viz_service: GridVizService, + well_provider: WellProvider, + well_server: WellServer, +) -> None: + @callback( + Output(get_uuid(LayoutElements.PROPERTIES), "options"), + Output(get_uuid(LayoutElements.PROPERTIES), "value"), + Input(get_uuid(LayoutElements.INIT_RESTART), "value"), + ) + def _populate_properties( + init_restart: str, + ) -> Tuple[ + List[Dict[str, str]], List[str], List[Dict[str, str]], Optional[List[str]] + ]: + if PROPERTYTYPE(init_restart) == PROPERTYTYPE.INIT: + prop_names = grid_provider.static_property_names() + + else: + prop_names = grid_provider.dynamic_property_names() + + return ( + [{"label": prop, "value": prop} for prop in prop_names], + [prop_names[0]], + ) + + @callback( + Output(get_uuid(LayoutElements.DATES), "options"), + Output(get_uuid(LayoutElements.DATES), "value"), + Input(get_uuid(LayoutElements.PROPERTIES), "value"), + State(get_uuid(LayoutElements.INIT_RESTART), "value"), + State(get_uuid(LayoutElements.DATES), "options"), + ) + def _populate_dates( + property_name: List[str], + init_restart: str, + current_date_options: List, + ) -> Tuple[List[Dict[str, str]], Optional[List[str]]]: + if PROPERTYTYPE(init_restart) == PROPERTYTYPE.INIT: + return [], None + else: + property_name = property_name[0] + dates = grid_provider.dates_for_dynamic_property( + property_name=property_name + ) + dates = dates if dates else [] + current_date_options = current_date_options if current_date_options else [] + if set(dates) == set( + [dateopt["value"] for dateopt in current_date_options] + ): + return no_update, no_update + return ( + ([{"label": prop, "value": prop} for prop in dates]), + [dates[0]] if dates else None, + ) + + @callback( + Output(get_uuid(LayoutElements.VTK_GRID_POLYDATA), "polys"), + Output(get_uuid(LayoutElements.VTK_GRID_POLYDATA), "points"), + Output(get_uuid(LayoutElements.VTK_GRID_CELLDATA), "values"), + Output(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "colorDataRange"), + Input(get_uuid(LayoutElements.PROPERTIES), "value"), + Input(get_uuid(LayoutElements.DATES), "value"), + Input(get_uuid(LayoutElements.REALIZATIONS), "value"), + Input(get_uuid(LayoutElements.GRID_RANGE_STORE), "data"), + State(get_uuid(LayoutElements.INIT_RESTART), "value"), + State(get_uuid(LayoutElements.VTK_GRID_POLYDATA), "polys"), + ) + def _set_geometry_and_scalar( + prop: List[str], + date: List[int], + realizations: List[int], + grid_range: List[List[int]], + proptype: str, + current_polys: str, + ) -> Tuple[Any, Any, Any, List, Any]: + + if PROPERTYTYPE(proptype) == PROPERTYTYPE.INIT: + property_spec = PropertySpec(prop_name=prop[0], prop_date=0) + else: + property_spec = PropertySpec(prop_name=prop[0], prop_date=date[0]) + + triggered = callback_context.triggered[0]["prop_id"] + timer = PerfTimer() + if ( + triggered == "." + or current_polys is None + or get_uuid(LayoutElements.GRID_RANGE_STORE) in triggered + or get_uuid(LayoutElements.REALIZATIONS) in triggered + ): + surface_polys, scalars = grid_viz_service.get_surface( + provider_id=grid_provider.provider_id(), + realization=realizations[0], + property_spec=property_spec, + cell_filter=CellFilter( + i_min=grid_range[0][0], + i_max=grid_range[0][1], + j_min=grid_range[1][0], + j_max=grid_range[1][1], + k_min=grid_range[2][0], + k_max=grid_range[2][1], + ), + ) + + return ( + b64_encode_numpy(surface_polys.poly_arr.astype(np.float32)), + b64_encode_numpy(surface_polys.point_arr.astype(np.float32)), + b64_encode_numpy(scalars.value_arr.astype(np.float32)), + [np.nanmin(scalars.value_arr), np.nanmax(scalars.value_arr)], + ) + else: + scalars = grid_viz_service.get_mapped_property_values( + provider_id=grid_provider.provider_id(), + realization=realizations[0], + property_spec=property_spec, + cell_filter=CellFilter( + i_min=grid_range[0][0], + i_max=grid_range[0][1], + j_min=grid_range[1][0], + j_max=grid_range[1][1], + k_min=grid_range[2][0], + k_max=grid_range[2][1], + ), + ) + return ( + no_update, + no_update, + b64_encode_numpy(scalars.value_arr.astype(np.float32)), + [np.nanmin(scalars.value_arr), np.nanmax(scalars.value_arr)], + ) + + @callback( + Output(get_uuid(LayoutElements.VTK_WELL_INTERSECT_POLYDATA), "points"), + Output(get_uuid(LayoutElements.VTK_WELL_INTERSECT_POLYDATA), "polys"), + Output(get_uuid(LayoutElements.VTK_WELL_INTERSECT_CELL_DATA), "values"), + Output(get_uuid(LayoutElements.VTK_WELL_2D_INTERSECT_POLYDATA), "points"), + Output(get_uuid(LayoutElements.VTK_WELL_2D_INTERSECT_POLYDATA), "polys"), + Output(get_uuid(LayoutElements.VTK_WELL_2D_INTERSECT_CELL_DATA), "values"), + Output(get_uuid(LayoutElements.VTK_INTERSECT_VIEW), "cameraPosition"), + Output(get_uuid(LayoutElements.VTK_INTERSECT_VIEW), "cameraFocalPoint"), + Output(get_uuid(LayoutElements.VTK_INTERSECT_VIEW), "cameraViewUp"), + Output(get_uuid(LayoutElements.VTK_INTERSECT_VIEW), "cameraParallelHorScale"), + Output(get_uuid(LayoutElements.LINEGRAPH), "figure"), + Input(get_uuid(LayoutElements.WELL_SELECT), "value"), + Input(get_uuid(LayoutElements.REALIZATIONS), "value"), + Input(get_uuid(LayoutElements.PROPERTIES), "value"), + Input(get_uuid(LayoutElements.DATES), "value"), + State(get_uuid(LayoutElements.INIT_RESTART), "value"), + ) + def set_well_geometries( + well_names: List[str], + realizations: List[int], + prop: List[str], + date: List[int], + proptype: str, + ) -> Tuple[ + List[Dict[str, str]], List[str], List[Dict[str, str]], Optional[List[str]] + ]: + + if not well_names: + return ( + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + no_update, + ) + polyline_xy = np.array( + well_provider.get_polyline_along_well_path_SIMPLIFIED(well_names[0]) + ) + polyline_xy_full = np.array( + well_provider.get_polyline_along_well_path_SIMPLIFIED( + well_names[0], use_rdp=False + ) + ) + + print(polyline_xy[:, 0], polyline_xy[:, 1]) + print(polyline_xy_full) + + def plotly_xy_plot(xy, xy2): + return { + "data": [ + { + "x": xy[:, 0], + "y": xy[:, 1], + "marker": dict( + size=20, + line=dict(color="MediumPurple", width=8), + ), + }, + {"x": xy2[:, 0], "y": xy2[:, 1]}, + ] + } + + if PROPERTYTYPE(proptype) == PROPERTYTYPE.INIT: + property_spec = PropertySpec(prop_name=prop[0], prop_date=0) + else: + property_spec = PropertySpec(prop_name=prop[0], prop_date=date[0]) + + surface_polys, scalars = grid_viz_service.cut_along_polyline( + provider_id=grid_provider.provider_id(), + realization=realizations[0], + polyline_xy=np.array(polyline_xy).flatten(), + property_spec=property_spec, + ) + + approx_plane_normal = _calc_approx_plane_normal_from_polyline_xy(polyline_xy) + + surf_points_3d = np.asarray(surface_polys.point_arr).reshape(-1, 3) + bb_min = np.min(surf_points_3d, axis=0) + bb_max = np.max(surf_points_3d, axis=0) + bb_radius = np.linalg.norm(bb_max - bb_min) / 2 + + center_pt = (bb_max + bb_min) / 2.0 + eye_pt = center_pt + bb_radius * approx_plane_normal + view_up_vec = [0.0, 0.0, 1.0] + + # Make scale slightly larger so we get some space on each side of the viewport + cameraParallelHorScale = bb_radius * 1.05 + + return ( + b64_encode_numpy(surface_polys.point_arr), + b64_encode_numpy(surface_polys.poly_arr), + b64_encode_numpy(scalars.value_arr) if scalars is not None else no_update, + b64_encode_numpy(surface_polys.point_arr), + b64_encode_numpy(surface_polys.poly_arr), + b64_encode_numpy(scalars.value_arr) if scalars is not None else no_update, + eye_pt, + center_pt, + view_up_vec, + cameraParallelHorScale, + plotly_xy_plot(polyline_xy, polyline_xy_full), + ) + + @callback( + Output(get_uuid(LayoutElements.VTK_WELL_PATH_POLYDATA), "points"), + Output(get_uuid(LayoutElements.VTK_WELL_PATH_POLYDATA), "lines"), + Output(get_uuid(LayoutElements.VTK_WELL_PATH_2D_POLYDATA), "points"), + Output(get_uuid(LayoutElements.VTK_WELL_PATH_2D_POLYDATA), "lines"), + Input(get_uuid(LayoutElements.WELL_SELECT), "value"), + ) + def set_well_geometries( + well_names: List[str], + ) -> Tuple[ + List[Dict[str, str]], List[str], List[Dict[str, str]], Optional[List[str]] + ]: + + if not well_names: + return no_update, no_update, no_update, no_update, no_update + polyline = well_server.get_polyline( + provider_id=well_provider.provider_id(), well_name=well_names[0] + ) + + return ( + b64_encode_numpy(polyline.point_arr.astype(np.float32)), + b64_encode_numpy(polyline.line_arr.astype(np.float32)), + b64_encode_numpy(polyline.point_arr.astype(np.float32)), + b64_encode_numpy(polyline.line_arr.astype(np.float32)), + ) + + @callback( + Output(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "actor"), + Output(get_uuid(LayoutElements.VTK_WELL_INTERSECT_REPRESENTATION), "actor"), + Output(get_uuid(LayoutElements.VTK_WELL_PATH_REPRESENTATION), "actor"), + Output(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "showCubeAxes"), + Input(get_uuid(LayoutElements.Z_SCALE), "value"), + Input(get_uuid(LayoutElements.SHOW_AXES), "value"), + ) + def _set_representation_actor( + z_scale: int, axes_is_on: List[str] + ) -> Tuple[dict, bool]: + show_axes = bool(z_scale == 1 and axes_is_on) + actor = {"scale": (1, 1, z_scale)} + return actor, actor, actor, show_axes + + @callback( + Output(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "property"), + Output(get_uuid(LayoutElements.VTK_WELL_INTERSECT_REPRESENTATION), "property"), + Input(get_uuid(LayoutElements.SHOW_GRID_LINES), "value"), + ) + def _set_representation_property( + show_grid_lines: List[str], + ) -> dict: + properties = {"edgeVisibility": bool(show_grid_lines)} + + return properties, properties + + @callback( + Output(get_uuid(LayoutElements.VTK_VIEW), "triggerResetCamera"), + Input(get_uuid(LayoutElements.REALIZATIONS), "value"), + Input(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "actor"), + ) + def _reset_camera(realizations: List[int], _actor: dict) -> float: + + return time() + + @callback( + Output(get_uuid(LayoutElements.SELECTED_CELL), "children"), + Output(get_uuid(LayoutElements.VTK_PICK_SPHERE), "state"), + Output(get_uuid(LayoutElements.VTK_PICK_REPRESENTATION), "actor"), + Input(get_uuid(LayoutElements.VTK_VIEW), "clickInfo"), + Input(get_uuid(LayoutElements.ENABLE_PICKING), "value"), + Input(get_uuid(LayoutElements.PROPERTIES), "value"), + Input(get_uuid(LayoutElements.DATES), "value"), + Input(get_uuid(LayoutElements.REALIZATIONS), "value"), + Input(get_uuid(LayoutElements.GRID_RANGE_STORE), "data"), + Input(get_uuid(LayoutElements.INIT_RESTART), "value"), + State(get_uuid(LayoutElements.Z_SCALE), "value"), + State(get_uuid(LayoutElements.VTK_PICK_REPRESENTATION), "actor"), + ) + # pylint: disable = too-many-locals, too-many-arguments + def _update_click_info( + click_data: Optional[Dict], + enable_picking: Optional[str], + prop: List[str], + date: List[int], + realizations: List[int], + grid_range: List[List[int]], + proptype: str, + zscale: float, + pick_representation_actor: Optional[Dict], + ) -> Tuple[str, Dict[str, Any], Dict[str, bool]]: + pick_representation_actor = ( + pick_representation_actor if pick_representation_actor else {} + ) + if not click_data: + return no_update, no_update, no_update + if not enable_picking: + pick_representation_actor.update({"visibility": False}) + return "", {}, pick_representation_actor + pick_representation_actor.update({"visibility": True}) + + client_world_pos = click_data["worldPosition"] + client_ray = click_data["ray"] + + # Remove z-scaling from client ray + client_world_pos[2] = client_world_pos[2] / zscale + client_ray[0][2] = client_ray[0][2] / zscale + client_ray[1][2] = client_ray[1][2] / zscale + + ray = Ray(origin=client_ray[0], end=client_ray[1]) + cell_filter = CellFilter( + i_min=grid_range[0][0], + i_max=grid_range[0][1], + j_min=grid_range[1][0], + j_max=grid_range[1][1], + k_min=grid_range[2][0], + k_max=grid_range[2][1], + ) + + if PROPERTYTYPE(proptype) == PROPERTYTYPE.INIT: + property_spec = PropertySpec(prop_name=prop[0], prop_date=0) + else: + property_spec = PropertySpec(prop_name=prop[0], prop_date=date[0]) + + pick_result = grid_viz_service.ray_pick( + provider_id=grid_provider.provider_id(), + realization=realizations[0], + ray=ray, + property_spec=property_spec, + cell_filter=cell_filter, + ) + + pick_sphere_pos = pick_result.intersection_point.copy() + pick_sphere_pos[2] *= zscale + + propname = f"{prop[0]}-{date[0]}" if date else f"{prop[0]}" + return ( + json.dumps( + { + "x": pick_result.intersection_point[0], + "y": pick_result.intersection_point[1], + "z": pick_result.intersection_point[2], + "i": pick_result.cell_i, + "j": pick_result.cell_j, + "k": pick_result.cell_k, + propname: float(pick_result.cell_property_value), + }, + indent=2, + ), + {"center": pick_sphere_pos, "radius": 100}, + pick_representation_actor, + ) + + @callback( + Output(get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), "colorMapPreset"), + Input(get_uuid(LayoutElements.COLORMAP), "value"), + ) + def _set_colormap(colormap: str) -> str: + return colormap + + @callback( + Output( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": MATCH, + "component": "input", + "component2": MATCH, + }, + "value", + ), + Output( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": MATCH, + "component": "slider", + "component2": MATCH, + }, + "value", + ), + Input( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": MATCH, + "component": "input", + "component2": MATCH, + }, + "value", + ), + Input( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": MATCH, + "component": "slider", + "component2": MATCH, + }, + "value", + ), + ) + def _synchronize_crop_slider_and_input( + input_val: int, slider_val: int + ) -> Tuple[Any, Any]: + trigger_id = callback_context.triggered[0]["prop_id"].split(".")[0] + if "slider" in trigger_id: + return slider_val, no_update + return no_update, input_val + + @callback( + Output(get_uuid(LayoutElements.GRID_RANGE_STORE), "data"), + Input( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": ALL, + "component": "input", + "component2": "start", + }, + "value", + ), + Input( + { + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": ALL, + "component": "input", + "component2": "width", + }, + "value", + ), + ) + def _store_grid_range_from_crop_widget( + input_vals: List[int], width_vals: List[int] + ) -> List[List[int]]: + if not input_vals or not width_vals: + return no_update + return [[val, val + width - 1] for val, width in zip(input_vals, width_vals)] + + +def _calc_approx_plane_normal_from_polyline_xy(polyline_xy: List[float]) -> List[float]: + polyline_np = np.asarray(polyline_xy).reshape(-1, 2) + num_points_in_polyline = len(polyline_np) + + aggr_right_vec = np.array([0.0, 0.0]) + for i in range(0, num_points_in_polyline - 1): + p0 = polyline_np[i] + p1 = polyline_np[i + 1] + fwd_vec = p1 - p0 + fwd_vec /= np.linalg.norm(fwd_vec) + right_vec = np.array([fwd_vec[1], -fwd_vec[0]]) + aggr_right_vec += right_vec + + avg_right_vec = aggr_right_vec / np.linalg.norm(aggr_right_vec) + approx_plane_normal = np.array([aggr_right_vec[0], aggr_right_vec[1], 0]) + + return approx_plane_normal diff --git a/webviz_subsurface/plugins/_eclipse_grid_viewer/_layout.py b/webviz_subsurface/plugins/_eclipse_grid_viewer/_layout.py new file mode 100644 index 000000000..2b9a77c2e --- /dev/null +++ b/webviz_subsurface/plugins/_eclipse_grid_viewer/_layout.py @@ -0,0 +1,542 @@ +from enum import Enum +from typing import Callable, Optional, List + +import webviz_vtk +import webviz_core_components as wcc +from dash import dcc, html + +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProvider, + CellFilter, +) +from webviz_subsurface._providers.well_provider import WellProvider + +# pylint: disable = too-few-public-methods +class LayoutElements(str, Enum): + REALIZATIONS = "realization" + INIT_RESTART = "init-restart-select" + PROPERTIES = "properties-select" + DATES = "dates-select" + WELL_SELECT = "well-select" + Z_SCALE = "z-scale" + VTK_VIEW = "vtk-view" + VTK_INTERSECT_VIEW = "vtk-intersect-view" + VTK_GRID_REPRESENTATION = "vtk-grid-representation" + VTK_GRID_POLYDATA = "vtk-grid-polydata" + VTK_GRID_CELLDATA = "vtk-grid-celldata" + VTK_WELL_INTERSECT_REPRESENTATION = "vtk-well-intersect-representation" + VTK_WELL_INTERSECT_POLYDATA = "vtk-well-intersect-polydata" + VTK_WELL_INTERSECT_CELL_DATA = "vtk-well-intersect-celldata" + VTK_WELL_PATH_REPRESENTATION = "vtk-well-path-representation" + VTK_WELL_PATH_POLYDATA = "vtk-well-path-polydata" + VTK_WELL_2D_INTERSECT_REPRESENTATION = "vtk-well-2d-intersect-representation" + VTK_WELL_2D_INTERSECT_POLYDATA = "vtk-well-2d-intersect-polydata" + VTK_WELL_2D_INTERSECT_CELL_DATA = "vtk-well-2d-intersect-celldata" + VTK_WELL_PATH_2D_REPRESENTATION = "vtk-well-2d-path-representation" + VTK_WELL_PATH_2D_POLYDATA = "vtk-well-2d-path-polydata" + STORED_CELL_INDICES_HASH = "stored-cell-indices-hash" + SELECTED_CELL = "selected-cell" + SHOW_GRID_LINES = "show-grid-lines" + COLORMAP = "color-map" + ENABLE_PICKING = "enable-picking" + VTK_PICK_REPRESENTATION = "vtk-pick-representation" + VTK_PICK_SPHERE = "vtk-pick-sphere" + SHOW_AXES = "show-axes" + CROP_WIDGET = "crop-widget" + GRID_RANGE_STORE = "crop-widget-store" + LINEGRAPH = "line-graph" + + +class LayoutTitles(str, Enum): + REALIZATIONS = "Realization" + INIT_RESTART = "Init / Restart" + PROPERTIES = "Property" + DATES = "Date" + WELL_SELECT = "Well" + Z_SCALE = "Z-scale" + SHOW_GRID_LINES = "Show grid lines" + COLORMAP = "Color map" + GRID_FILTERS = "Grid filters" + COLORS = "Colors" + PICKING = "Picking" + ENABLE_PICKING = "Enable readout from picked cell" + SHOW_AXES = "Show axes" + + +class GRID_DIRECTION(str, Enum): + I = "I" + J = "J" + K = "K" + + +COLORMAPS = ["erdc_rainbow_dark", "Viridis (matplotlib)", "BuRd"] + + +class PROPERTYTYPE(str, Enum): + INIT = "Init" + RESTART = "Restart" + + +class LayoutStyle: + MAIN_HEIGHT = "87vh" + SIDEBAR = {"flex": 1, "height": "87vh"} + VTK_VIEW = {"height": "40vh", "marginBottom": "10px"} + + +def plugin_main_layout( + get_uuid: Callable, grid_provider: EnsembleGridProvider, well_names: List[str] +) -> wcc.FlexBox: + initial_grid = grid_provider.get_3dgrid(grid_provider.realizations()[0]) + grid_dimensions = CellFilter( + i_min=0, + j_min=0, + k_min=0, + i_max=initial_grid.dimensions[0] - 1, + j_max=initial_grid.dimensions[1] - 1, + k_max=initial_grid.dimensions[2] - 1, + ) + return wcc.FlexBox( + children=[ + sidebar( + get_uuid=get_uuid, + grid_dimensions=grid_dimensions, + grid_provider=grid_provider, + well_names=well_names, + ), + html.Div( + style={"flex": "5"}, + children=[ + vtk_3d_view(get_uuid=get_uuid), + wcc.FlexBox( + children=[ + html.Div( + style={"flex": 1}, + children=vtk_intersect_view(get_uuid=get_uuid), + ), + html.Div( + style={"flex": 1}, + children=dcc.Graph( + id=get_uuid(LayoutElements.LINEGRAPH) + ), + ), + ] + ), + ], + ), + dcc.Store(id=get_uuid(LayoutElements.STORED_CELL_INDICES_HASH)), + dcc.Store( + id=get_uuid(LayoutElements.GRID_RANGE_STORE), + data=[ + [grid_dimensions.i_min, grid_dimensions.i_max], + [grid_dimensions.j_min, grid_dimensions.j_max], + [grid_dimensions.k_min, grid_dimensions.k_min], + ], + ), + ] + ) + + +def sidebar( + get_uuid: Callable, + grid_dimensions: CellFilter, + grid_provider: EnsembleGridProvider, + well_names: List[str], +) -> wcc.Frame: + + realizations = grid_provider.realizations() + property_options = [] + property_value = None + + if grid_provider.static_property_names(): + property_options.append( + {"label": PROPERTYTYPE.INIT, "value": PROPERTYTYPE.INIT} + ) + property_value = PROPERTYTYPE.INIT + + if grid_provider.dynamic_property_names(): + property_options.append( + {"label": PROPERTYTYPE.RESTART, "value": PROPERTYTYPE.RESTART} + ) + if property_value is None: + property_value = PROPERTYTYPE.RESTART + + return wcc.Frame( + style=LayoutStyle.SIDEBAR, + children=[ + wcc.SelectWithLabel( + id=get_uuid(LayoutElements.REALIZATIONS), + label=LayoutTitles.REALIZATIONS, + options=[{"label": real, "value": real} for real in realizations], + value=[realizations[0]], + multi=False, + ), + wcc.RadioItems( + label=LayoutTitles.INIT_RESTART, + id=get_uuid(LayoutElements.INIT_RESTART), + options=property_options, + value=property_value, + ), + wcc.SelectWithLabel( + id=get_uuid(LayoutElements.PROPERTIES), + label=LayoutTitles.PROPERTIES, + ), + wcc.SelectWithLabel( + id=get_uuid(LayoutElements.DATES), label=LayoutTitles.DATES + ), + wcc.SelectWithLabel( + id=get_uuid(LayoutElements.WELL_SELECT), + label=LayoutTitles.WELL_SELECT, + multi=False, + options=[{"value": well, "label": well} for well in well_names], + value=[well_names[0]], + ), + wcc.Slider( + label=LayoutTitles.Z_SCALE, + id=get_uuid(LayoutElements.Z_SCALE), + min=1, + max=10, + value=1, + step=1, + ), + wcc.Selectors( + label=LayoutTitles.COLORS, + children=[ + wcc.Dropdown( + id=get_uuid(LayoutElements.COLORMAP), + options=[ + {"value": colormap, "label": colormap} + for colormap in COLORMAPS + ], + value=COLORMAPS[0], + clearable=False, + ) + ], + ), + wcc.Selectors( + label="Range filters", + children=[ + crop_widget( + get_uuid=get_uuid, + min_val=grid_dimensions.i_min, + max_val=grid_dimensions.i_max, + direction=GRID_DIRECTION.I, + ), + crop_widget( + get_uuid=get_uuid, + min_val=grid_dimensions.j_min, + max_val=grid_dimensions.j_max, + direction=GRID_DIRECTION.J, + ), + crop_widget( + get_uuid=get_uuid, + min_val=grid_dimensions.k_min, + max_val=grid_dimensions.k_max, + max_width=grid_dimensions.k_min, + direction=GRID_DIRECTION.K, + ), + ], + ), + wcc.Selectors( + label="Options", + children=[ + wcc.Checklist( + id=get_uuid(LayoutElements.SHOW_AXES), + options=[LayoutTitles.SHOW_AXES], + value=[LayoutTitles.SHOW_AXES], + ), + wcc.Checklist( + id=get_uuid(LayoutElements.SHOW_GRID_LINES), + options=[LayoutTitles.SHOW_GRID_LINES], + value=[LayoutTitles.SHOW_GRID_LINES], + ), + ], + ), + wcc.Selectors( + label="Readout", + children=[ + wcc.Checklist( + id=get_uuid(LayoutElements.ENABLE_PICKING), + options=[LayoutTitles.ENABLE_PICKING], + value=[LayoutTitles.ENABLE_PICKING], + ) + ], + ), + html.Pre(id=get_uuid(LayoutElements.SELECTED_CELL), children=[]), + ], + ) + + +def crop_widget( + get_uuid: Callable, + min_val: int, + max_val: int, + direction: str, + max_width: Optional[int] = None, +) -> html.Div: + max_width = max_width if max_width else max_val + return html.Div( + children=[ + html.Div( + style={ + "display": "grid", + "marginBotton": "0px", + "gridTemplateColumns": f"2fr 1fr 8fr", + }, + children=[ + wcc.Label( + children=f"{direction} Start", + style={ + "fontSize": "0.7em", + "fontWeight": "bold", + "marginRight": "5px", + }, + ), + dcc.Input( + style={"width": "50px", "height": "10px"}, + id={ + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": direction, + "component": "input", + "component2": "start", + }, + type="number", + placeholder="Min", + persistence=True, + persistence_type="session", + value=min_val, + min=min_val, + max=max_val, + ), + wcc.Slider( + id={ + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": direction, + "component": "slider", + "component2": "start", + }, + min=min_val, + max=max_val, + value=min_val, + step=1, + marks=None, + ), + ], + ), + html.Div( + style={ + "display": "grid", + "marginTop": "0px", + "padding": "0px", + "gridTemplateColumns": f"2fr 1fr 8fr", + }, + children=[ + wcc.Label( + children=f"Width", + style={ + "fontSize": "0.7em", + "textAlign": "right", + "marginRight": "5px", + }, + ), + dcc.Input( + style={"width": "50px", "height": "10px"}, + id={ + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": direction, + "component": "input", + "component2": "width", + }, + type="number", + placeholder="Min", + persistence=True, + persistence_type="session", + value=max_width, + min=1, + max=max_val, + ), + wcc.Slider( + id={ + "id": get_uuid(LayoutElements.CROP_WIDGET), + "direction": direction, + "component": "slider", + "component2": "width", + }, + min=1, + max=max_val, + value=max_width, + step=1, + marks=None, + ), + ], + ), + ], + ) + + +def vtk_3d_view(get_uuid: Callable) -> webviz_vtk.View: + return webviz_vtk.View( + id=get_uuid(LayoutElements.VTK_VIEW), + style=LayoutStyle.VTK_VIEW, + pickingModes=["click"], + interactorSettings=[ + { + "button": 1, + "action": "Zoom", + "scrollEnabled": True, + }, + { + "button": 3, + "action": "Pan", + }, + { + "button": 2, + "action": "Rotate", + }, + { + "button": 1, + "action": "Pan", + "shift": True, + }, + { + "button": 1, + "action": "Zoom", + "alt": True, + }, + { + "button": 1, + "action": "Roll", + "alt": True, + "shift": True, + }, + ], + children=[ + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_GRID_REPRESENTATION), + showCubeAxes=True, + showScalarBar=True, + children=[ + webviz_vtk.PolyData( + id=get_uuid(LayoutElements.VTK_GRID_POLYDATA), + children=[ + webviz_vtk.CellData( + [ + webviz_vtk.DataArray( + id=get_uuid(LayoutElements.VTK_GRID_CELLDATA), + registration="setScalars", + name="scalar", + ) + ] + ) + ], + ) + ], + property={"edgeVisibility": True}, + ), + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_PICK_REPRESENTATION), + actor={"visibility": False}, + children=[ + webviz_vtk.Algorithm( + id=get_uuid(LayoutElements.VTK_PICK_SPHERE), + vtkClass="vtkSphereSource", + ) + ], + ), + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_WELL_INTERSECT_REPRESENTATION), + actor={"visibility": True}, + children=[ + webviz_vtk.PolyData( + id=get_uuid(LayoutElements.VTK_WELL_INTERSECT_POLYDATA), + children=[ + webviz_vtk.CellData( + [ + webviz_vtk.DataArray( + id=get_uuid( + LayoutElements.VTK_WELL_INTERSECT_CELL_DATA + ), + registration="setScalars", + name="scalar", + ) + ] + ) + ], + ) + ], + ), + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_WELL_PATH_REPRESENTATION), + actor={"visibility": True}, + children=[ + webviz_vtk.PolyData( + id=get_uuid(LayoutElements.VTK_WELL_PATH_POLYDATA), + ) + ], + ), + ], + ) + + +def vtk_intersect_view(get_uuid: Callable) -> webviz_vtk.View: + + # fmt: off + interactorSettings = [ + {"button": 1, "action": "Zoom"}, + {"button": 1, "action": "Zoom", "alt": True}, + {"button": 1, "action": "Pan", "shift": True}, + {"button": 3, "action": "Pan"}, + #{"button": 2, "action": "Rotate", "useFocalPointAsCenterOfRotation": True}, + {"dragEnabled": False, "action": "Zoom", "scrollEnabled": True}, + {"dragEnabled": False, "action": "ZoomToMouse", "scrollEnabled": True, "control": True,}, + ] + # fmt: on + + return webviz_vtk.View( + id=get_uuid(LayoutElements.VTK_INTERSECT_VIEW), + style=LayoutStyle.VTK_VIEW, + pickingModes=["click"], + cameraParallelProjection=True, + autoResetCamera=False, + interactorSettings=interactorSettings, + children=[ + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_WELL_2D_INTERSECT_REPRESENTATION), + actor={"visibility": True}, + property={ + "edgeVisibility": False, + "lighting": False, + }, + children=[ + webviz_vtk.PolyData( + id=get_uuid(LayoutElements.VTK_WELL_2D_INTERSECT_POLYDATA), + children=[ + webviz_vtk.CellData( + [ + webviz_vtk.DataArray( + id=get_uuid( + LayoutElements.VTK_WELL_2D_INTERSECT_CELL_DATA + ), + registration="setScalars", + name="scalar", + ) + ] + ) + ], + ) + ], + ), + webviz_vtk.GeometryRepresentation( + id=get_uuid(LayoutElements.VTK_WELL_PATH_2D_REPRESENTATION), + actor={"visibility": True}, + children=[ + webviz_vtk.PolyData( + id=get_uuid(LayoutElements.VTK_WELL_PATH_2D_POLYDATA), + ) + ], + ), + ], + ) + + +def plotly_xy_plot(xy, xy2): + return {"data": [{}]} diff --git a/webviz_subsurface/plugins/_eclipse_grid_viewer/_plugin.py b/webviz_subsurface/plugins/_eclipse_grid_viewer/_plugin.py new file mode 100644 index 000000000..0bc1c24d0 --- /dev/null +++ b/webviz_subsurface/plugins/_eclipse_grid_viewer/_plugin.py @@ -0,0 +1,71 @@ +from pathlib import Path +from typing import Callable, Dict, List, Tuple +from dash import html +import webviz_core_components as wcc +from webviz_config import WebvizPluginABC, WebvizSettings +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProviderFactory, + GridVizService, +) +from webviz_subsurface._providers.well_provider import ( + WellProvider, + WellProviderFactory, + WellServer, +) + +from ._callbacks import plugin_callbacks +from ._layout import plugin_main_layout + + +class EclipseGridViewer(WebvizPluginABC): + """Eclipse grid viewer""" + + def __init__( + self, + webviz_settings: WebvizSettings, + app, + ensembles: List[str], + grid_name: str, + well_folder: Path = None, + well_suffix: str = ".rmswell", + ) -> None: + super().__init__() + grid_provider_factory = EnsembleGridProviderFactory.instance() + self.grid_provider = grid_provider_factory.create_from_roff_files( + ens_path=webviz_settings.shared_settings["scratch_ensembles"][ensembles[0]], + grid_name=grid_name, + ) + self.grid_viz_service = GridVizService.instance() + self.grid_viz_service.register_provider(self.grid_provider) + factory = WellProviderFactory.instance() + + if well_folder is not None: + self.well_provider = factory.create_from_well_files( + well_folder=str(well_folder), + well_suffix=".rmswell", + md_logname="MDepth", + ) + + self.well_server = WellServer.instance(app) + self.well_server.register_provider(self.well_provider) + else: + self.well_provider = None + self.well_server = None + + plugin_callbacks( + get_uuid=self.uuid, + grid_provider=self.grid_provider, + grid_viz_service=self.grid_viz_service, + well_provider=self.well_provider, + well_server=self.well_server, + ) + + @property + def layout(self) -> wcc.FlexBox: + return plugin_main_layout( + get_uuid=self.uuid, + grid_provider=self.grid_provider, + well_names=self.well_provider.well_names() + if self.well_provider is not None + else [], + ) diff --git a/webviz_subsurface/plugins/_grid_viewer/__init__.py b/webviz_subsurface/plugins/_grid_viewer/__init__.py new file mode 100644 index 000000000..181db1d9d --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/__init__.py @@ -0,0 +1 @@ +from ._plugin import GridViewer diff --git a/webviz_subsurface/plugins/_grid_viewer/_layout_elements.py b/webviz_subsurface/plugins/_grid_viewer/_layout_elements.py new file mode 100644 index 000000000..ae8e20f8d --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/_layout_elements.py @@ -0,0 +1,30 @@ +class ElementIds: + ID = "grid-view" + + class VTKVIEW3D: + ID = "vtk" + VIEW = "vtk-view" + GRID_REPRESENTATION = "grid-representation" + GRID_POLYDATA = "grid-polydata" + GRID_CELLDATA = "grid-celldata" + PICK_REPRESENTATION = "pick-representation" + PICK_SPHERE = "pick-sphere" + + class DataSelectors: + ID = "data-selectors" + REALIZATIONS = "realizations" + STATIC_DYNAMIC = "static-dynamic" + PROPERTIES = "properties" + DATES = "dates" + WELLS = "wells" + + class GridFilter: + ID = "grid-filter" + IJK_CROP_STORE = "ijk-filter-store" + IJK_CROP_WIDGET = "ijk-crop-widget" + + class Settings: + ID = "settings" + ZSCALE = "z-scale" + COLORMAP = "colormap" + SHOW_CUBEAXES = "show-cube-axes" diff --git a/webviz_subsurface/plugins/_grid_viewer/_plugin.py b/webviz_subsurface/plugins/_grid_viewer/_plugin.py new file mode 100644 index 000000000..e241e6a6b --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/_plugin.py @@ -0,0 +1,56 @@ +from typing import List, Optional +from pathlib import Path + +from webviz_config import WebvizPluginABC, WebvizSettings + +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProviderFactory, + GridVizService, + EnsembleGridProvider, +) +from webviz_subsurface._providers.well_provider import ( + WellProvider, + WellProviderFactory, + WellServer, +) + +from ._layout_elements import ElementIds +from .views.view_3d._view_3d import View3D + + +class GridViewer(WebvizPluginABC): + def __init__( + self, + webviz_settings: WebvizSettings, + ensembles: List[str], + grid_name: str, + ): + super().__init__(stretch=True) + + self.ensembles = { + ens_name: webviz_settings.shared_settings["scratch_ensembles"][ens_name] + for ens_name in ensembles + } + + self.add_grid_provider(grid_name=grid_name) + + self.add_store( + ElementIds.GridFilter.IJK_CROP_STORE, + storage_type=WebvizPluginABC.StorageType.SESSION, + ) + print(self._stores) + self.add_view( + View3D( + grid_provider=self.grid_provider, grid_viz_service=self.grid_viz_service + ), + ElementIds.ID, + ) + + def add_grid_provider(self, grid_name: str) -> None: + factory = EnsembleGridProviderFactory.instance() + self.grid_provider: EnsembleGridProvider = factory.create_from_roff_files( + ens_path=list(self.ensembles.values())[0], + grid_name=grid_name, + ) + self.grid_viz_service = GridVizService.instance() + self.grid_viz_service.register_provider(self.grid_provider) diff --git a/webviz_subsurface/plugins/_grid_viewer/_types.py b/webviz_subsurface/plugins/_grid_viewer/_types.py new file mode 100644 index 000000000..2c1cca78e --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/_types.py @@ -0,0 +1,12 @@ +from enum import Enum + + +class PROPERTYTYPE(str, Enum): + STATIC = "Static" + DYNAMIC = "Dynamic" + + +class GRID_DIRECTION(str, Enum): + I = "I" + J = "J" + K = "K" diff --git a/webviz_subsurface/plugins/_grid_viewer/views/__init__.py b/webviz_subsurface/plugins/_grid_viewer/views/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/__init__.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/_view_3d.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/_view_3d.py new file mode 100644 index 000000000..a7b27f82a --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/_view_3d.py @@ -0,0 +1,177 @@ +from typing import List, Tuple, Dict, Optional + +import numpy as np +from webviz_config.webviz_plugin_subclasses import ( + ViewABC, +) + +from dash.development.base_component import Component + +from dash import Input, Output, State, callback, callback_context, no_update, ALL +import webviz_core_components as wcc +from webviz_vtk.utils.vtk import b64_encode_numpy + +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProvider, + GridVizService, + PropertySpec, + CellFilter, + Ray, +) +from webviz_subsurface._utils.perf_timer import PerfTimer +from ..._layout_elements import ElementIds +from ..._types import PROPERTYTYPE +from .settings import DataSettings, GridFilter, Settings +from .view_elements._vtk_view_3d_element import VTKView3D + + +class View3D(ViewABC): + def __init__( + self, grid_provider: EnsembleGridProvider, grid_viz_service: GridVizService + ) -> None: + super().__init__("Grid View") + self.grid_provider = grid_provider + self.grid_viz_service = grid_viz_service + self.vtk_view_3d = VTKView3D() + self.add_view_element(self.vtk_view_3d, ElementIds.VTKVIEW3D.ID) + self.add_settings_group( + DataSettings(grid_provider=grid_provider), ElementIds.DataSelectors.ID + ) + self.add_settings_group( + GridFilter(grid_provider=grid_provider), ElementIds.GridFilter.ID + ) + + def set_callbacks(self) -> None: + @callback( + Output( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.GRID_POLYDATA) + .to_string(), + "polys", + ), + Output( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.GRID_POLYDATA) + .to_string(), + "points", + ), + Output( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.GRID_CELLDATA) + .to_string(), + "values", + ), + Output( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.GRID_REPRESENTATION) + .to_string(), + "colorDataRange", + ), + Output( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.VIEW) + .to_string(), + "triggerResetCamera", + ), + Input( + self.settings_group(ElementIds.DataSelectors.ID) + .component_unique_id(ElementIds.DataSelectors.PROPERTIES) + .to_string(), + "value", + ), + Input( + self.settings_group(ElementIds.DataSelectors.ID) + .component_unique_id(ElementIds.DataSelectors.DATES) + .to_string(), + "value", + ), + Input( + self.settings_group(ElementIds.DataSelectors.ID) + .component_unique_id(ElementIds.DataSelectors.REALIZATIONS) + .to_string(), + "value", + ), + Input( + self.get_store_unique_id(ElementIds.GridFilter.IJK_CROP_STORE), "data" + ), + State( + self.settings_group(ElementIds.DataSelectors.ID) + .component_unique_id(ElementIds.DataSelectors.STATIC_DYNAMIC) + .to_string(), + "value", + ), + State( + self.view_element(ElementIds.VTKVIEW3D.ID) + .component_unique_id(ElementIds.VTKVIEW3D.GRID_POLYDATA) + .to_string(), + "polys", + ), + ) + def _set_geometry_and_scalar( + prop: List[str], + date: List[int], + realizations: List[int], + grid_range: List[List[int]], + proptype: str, + current_polys: str, + ) -> Tuple: + + if PROPERTYTYPE(proptype) == PROPERTYTYPE.STATIC: + property_spec = PropertySpec(prop_name=prop[0], prop_date=0) + else: + property_spec = PropertySpec(prop_name=prop[0], prop_date=date[0]) + + triggered = callback_context.triggered[0]["prop_id"] + timer = PerfTimer() + if ( + triggered == "." + or current_polys is None + or self.get_store_unique_id(ElementIds.GridFilter.IJK_CROP_STORE) + in triggered + or self.settings_group(ElementIds.DataSelectors.ID) + .component_unique_id(ElementIds.DataSelectors.REALIZATIONS) + .to_string() + in triggered + ): + surface_polys, scalars = self.grid_viz_service.get_surface( + provider_id=self.grid_provider.provider_id(), + realization=realizations[0], + property_spec=property_spec, + cell_filter=CellFilter( + i_min=grid_range[0][0], + i_max=grid_range[0][1], + j_min=grid_range[1][0], + j_max=grid_range[1][1], + k_min=grid_range[2][0], + k_max=grid_range[2][1], + ), + ) + + return ( + b64_encode_numpy(surface_polys.poly_arr.astype(np.float32)), + b64_encode_numpy(surface_polys.point_arr.astype(np.float32)), + b64_encode_numpy(scalars.value_arr.astype(np.float32)), + [np.nanmin(scalars.value_arr), np.nanmax(scalars.value_arr)], + timer.elapsed_ms(), + ) + else: + scalars = self.grid_viz_service.get_mapped_property_values( + provider_id=self.grid_provider.provider_id(), + realization=realizations[0], + property_spec=property_spec, + cell_filter=CellFilter( + i_min=grid_range[0][0], + i_max=grid_range[0][1], + j_min=grid_range[1][0], + j_max=grid_range[1][1], + k_min=grid_range[2][0], + k_max=grid_range[2][1], + ), + ) + return ( + no_update, + no_update, + b64_encode_numpy(scalars.value_arr.astype(np.float32)), + [np.nanmin(scalars.value_arr), np.nanmax(scalars.value_arr)], + no_update, + ) diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/__init__.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/__init__.py new file mode 100644 index 000000000..00f236001 --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/__init__.py @@ -0,0 +1,3 @@ +from ._data_selection import DataSettings +from ._grid_filter import GridFilter +from ._settings import Settings diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_data_selection.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_data_selection.py new file mode 100644 index 000000000..03c7bc0ed --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_data_selection.py @@ -0,0 +1,157 @@ +from typing import List, Tuple, Dict, Optional + +from dash.development.base_component import Component +from dash import Input, Output, State, callback, no_update +import webviz_core_components as wcc +from webviz_config.webviz_plugin_subclasses import SettingsGroupABC +from webviz_subsurface._providers.ensemble_grid_provider import EnsembleGridProvider + +from webviz_subsurface.plugins._grid_viewer._types import PROPERTYTYPE +from webviz_subsurface.plugins._grid_viewer._layout_elements import ElementIds + + +def list_to_options(values: List) -> List: + return [{"value": val, "label": val} for val in values] + + +class DataSettings(SettingsGroupABC): + def __init__(self, grid_provider: EnsembleGridProvider) -> None: + super().__init__("Data Selection") + self.grid_provider = grid_provider + self.static_dynamic_options = [] + self.static_dynamic_value = None + + if grid_provider.static_property_names(): + self.static_dynamic_options.append( + {"label": PROPERTYTYPE.STATIC, "value": PROPERTYTYPE.STATIC} + ) + self.static_dynamic_value = PROPERTYTYPE.STATIC + + if grid_provider.dynamic_property_names(): + self.static_dynamic_options.append( + {"label": PROPERTYTYPE.DYNAMIC, "value": PROPERTYTYPE.DYNAMIC} + ) + if self.static_dynamic_value is None: + self.static_dynamic_value = PROPERTYTYPE.DYNAMIC + + def layout(self) -> List[Component]: + + return [ + wcc.SelectWithLabel( + label="Realizations", + id=self.register_component_unique_id( + ElementIds.DataSelectors.REALIZATIONS + ), + multi=False, + options=list_to_options(self.grid_provider.realizations()), + value=[self.grid_provider.realizations()[0]], + ), + wcc.RadioItems( + label="Static / Dynamic", + id=self.register_component_unique_id( + ElementIds.DataSelectors.STATIC_DYNAMIC + ), + options=self.static_dynamic_options, + value=self.static_dynamic_value, + ), + wcc.SelectWithLabel( + label="Property", + id=self.register_component_unique_id( + ElementIds.DataSelectors.PROPERTIES + ), + multi=False, + ), + wcc.SelectWithLabel( + label="Date", + id=self.register_component_unique_id(ElementIds.DataSelectors.DATES), + multi=False, + ), + ] + + def set_callbacks(self) -> None: + @callback( + Output( + self.component_unique_id( + ElementIds.DataSelectors.PROPERTIES + ).to_string(), + "options", + ), + Output( + self.component_unique_id( + ElementIds.DataSelectors.PROPERTIES + ).to_string(), + "value", + ), + Input( + self.component_unique_id( + ElementIds.DataSelectors.STATIC_DYNAMIC + ).to_string(), + "value", + ), + ) + def _populate_properties( + static_dynamic: str, + ) -> Tuple[ + List[Dict[str, str]], List[str], List[Dict[str, str]], Optional[List[str]] + ]: + if PROPERTYTYPE(static_dynamic) == PROPERTYTYPE.STATIC: + prop_names = self.grid_provider.static_property_names() + + else: + prop_names = self.grid_provider.dynamic_property_names() + + return ( + [{"label": prop, "value": prop} for prop in prop_names], + [prop_names[0]], + ) + + @callback( + Output( + self.component_unique_id(ElementIds.DataSelectors.DATES).to_string(), + "options", + ), + Output( + self.component_unique_id(ElementIds.DataSelectors.DATES).to_string(), + "value", + ), + Input( + self.component_unique_id( + ElementIds.DataSelectors.PROPERTIES + ).to_string(), + "value", + ), + State( + self.component_unique_id( + ElementIds.DataSelectors.STATIC_DYNAMIC + ).to_string(), + "value", + ), + State( + self.component_unique_id(ElementIds.DataSelectors.DATES).to_string(), + "options", + ), + ) + def _populate_dates( + property_name: List[str], + static_dynamic: str, + current_date_options: List, + ) -> Tuple[List[Dict[str, str]], Optional[List[str]]]: + if PROPERTYTYPE(static_dynamic) == PROPERTYTYPE.STATIC: + return [], None + else: + property_name = property_name[0] + dates = self.grid_provider.dates_for_dynamic_property( + property_name=property_name + ) + dates = dates if dates else [] + current_date_options = ( + current_date_options if current_date_options else [] + ) + if set(dates) == set( + [dateopt["value"] for dateopt in current_date_options] + ): + return no_update, no_update + return ( + ([{"label": prop, "value": prop} for prop in dates]), + [dates[0]] if dates else None, + ) diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_grid_filter.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_grid_filter.py new file mode 100644 index 000000000..5d1ea0b05 --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_grid_filter.py @@ -0,0 +1,257 @@ +from typing import List, Tuple, Dict, Optional, Any + +from dash.development.base_component import Component +from dash import ( + html, + dcc, + Input, + Output, + callback, + no_update, + ALL, + MATCH, + callback_context, +) +import webviz_core_components as wcc +from webviz_config.webviz_plugin_subclasses import SettingsGroupABC +from webviz_subsurface._providers.ensemble_grid_provider import ( + EnsembleGridProvider, + CellFilter, +) + +from webviz_subsurface.plugins._grid_viewer._types import GRID_DIRECTION +from webviz_subsurface.plugins._grid_viewer._layout_elements import ElementIds + + +def list_to_options(values: List) -> List: + return [{"value": val, "label": val} for val in values] + + +class GridFilter(SettingsGroupABC): + def __init__(self, grid_provider: EnsembleGridProvider) -> None: + super().__init__("Grid IJK Filter") + self.grid_provider = grid_provider + initial_grid = grid_provider.get_3dgrid(grid_provider.realizations()[0]) + self.grid_dimensions = CellFilter( + i_min=0, + j_min=0, + k_min=0, + i_max=initial_grid.dimensions[0] - 1, + j_max=initial_grid.dimensions[1] - 1, + k_max=initial_grid.dimensions[2] - 1, + ) + + def layout(self) -> List[Component]: + + return [ + wcc.Selectors( + label="Range filters", + children=[ + crop_widget( + widget_id=self.get_unique_id().to_string(), + min_val=self.grid_dimensions.i_min, + max_val=self.grid_dimensions.i_max, + direction=GRID_DIRECTION.I, + ), + crop_widget( + widget_id=self.get_unique_id().to_string(), + min_val=self.grid_dimensions.j_min, + max_val=self.grid_dimensions.j_max, + direction=GRID_DIRECTION.J, + ), + crop_widget( + widget_id=self.get_unique_id().to_string(), + min_val=self.grid_dimensions.k_min, + max_val=self.grid_dimensions.k_max, + max_width=self.grid_dimensions.k_min, + direction=GRID_DIRECTION.K, + ), + ], + ), + ] + + def set_callbacks(self) -> None: + @callback( + Output( + { + "id": self.get_unique_id().to_string(), + "direction": MATCH, + "component": "input", + "component2": MATCH, + }, + "value", + ), + Output( + { + "id": self.get_unique_id().to_string(), + "direction": MATCH, + "component": "slider", + "component2": MATCH, + }, + "value", + ), + Input( + { + "id": self.get_unique_id().to_string(), + "direction": MATCH, + "component": "input", + "component2": MATCH, + }, + "value", + ), + Input( + { + "id": self.get_unique_id().to_string(), + "direction": MATCH, + "component": "slider", + "component2": MATCH, + }, + "value", + ), + ) + def _synchronize_crop_slider_and_input( + input_val: int, slider_val: int + ) -> Tuple[Any, Any]: + trigger_id = callback_context.triggered[0]["prop_id"].split(".")[0] + if "slider" in trigger_id: + return slider_val, no_update + return no_update, input_val + + @callback( + Output( + self.get_store_unique_id(ElementIds.GridFilter.IJK_CROP_STORE), "data" + ), + Input( + { + "id": self.get_unique_id().to_string(), + "direction": ALL, + "component": "input", + "component2": "start", + }, + "value", + ), + Input( + { + "id": self.get_unique_id().to_string(), + "direction": ALL, + "component": "input", + "component2": "width", + }, + "value", + ), + ) + def _store_grid_range_from_crop_widget( + input_vals: List[int], width_vals: List[int] + ) -> List[List[int]]: + if not input_vals or not width_vals: + return no_update + return [ + [val, val + width - 1] for val, width in zip(input_vals, width_vals) + ] + + +def crop_widget( + widget_id: str, + min_val: int, + max_val: int, + direction: str, + max_width: Optional[int] = None, +) -> html.Div: + max_width = max_width if max_width else max_val + return html.Div( + children=[ + html.Div( + style={ + "display": "grid", + "marginBotton": "0px", + "gridTemplateColumns": f"2fr 1fr 8fr", + }, + children=[ + wcc.Label( + children=f"{direction} Start", + style={ + "fontSize": "0.7em", + "fontWeight": "bold", + "marginRight": "5px", + }, + ), + dcc.Input( + style={"width": "50px", "height": "10px"}, + id={ + "id": widget_id, + "direction": direction, + "component": "input", + "component2": "start", + }, + type="number", + placeholder="Min", + persistence=True, + persistence_type="session", + value=min_val, + min=min_val, + max=max_val, + ), + wcc.Slider( + id={ + "id": widget_id, + "direction": direction, + "component": "slider", + "component2": "start", + }, + min=min_val, + max=max_val, + value=min_val, + step=1, + marks=None, + ), + ], + ), + html.Div( + style={ + "display": "grid", + "marginTop": "0px", + "padding": "0px", + "gridTemplateColumns": f"2fr 1fr 8fr", + }, + children=[ + wcc.Label( + children=f"Width", + style={ + "fontSize": "0.7em", + "textAlign": "right", + "marginRight": "5px", + }, + ), + dcc.Input( + style={"width": "50px", "height": "10px"}, + id={ + "id": widget_id, + "direction": direction, + "component": "input", + "component2": "width", + }, + type="number", + placeholder="Min", + persistence=True, + persistence_type="session", + value=max_width, + min=1, + max=max_val, + ), + wcc.Slider( + id={ + "id": widget_id, + "direction": direction, + "component": "slider", + "component2": "width", + }, + min=1, + max=max_val, + value=max_width, + step=1, + marks=None, + ), + ], + ), + ], + ) diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_settings.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_settings.py new file mode 100644 index 000000000..52fb7c89c --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/settings/_settings.py @@ -0,0 +1,51 @@ +from typing import List, Tuple, Dict, Optional + +from dash.development.base_component import Component +from dash import Input, Output, State, callback, no_update +import webviz_core_components as wcc +from webviz_config.webviz_plugin_subclasses import SettingsGroupABC +from webviz_subsurface._providers.ensemble_grid_provider import EnsembleGridProvider + +from webviz_subsurface.plugins._grid_viewer._types import PROPERTYTYPE +from webviz_subsurface.plugins._grid_viewer._layout_elements import ElementIds + + +def list_to_options(values: List) -> List: + return [{"value": val, "label": val} for val in values] + + +class Settings(SettingsGroupABC): + def __init__(self) -> None: + super().__init__("Settings") + self.colormaps = ["erdc_rainbow_dark", "Viridis (matplotlib)", "BuRd"] + + def layout(self) -> List[Component]: + + return [ + wcc.Slider( + label="Z Scale", + id=self.register_component_unique_id(ElementIds.Settings.ZSCALE), + min=1, + max=10, + value=1, + step=1, + ), + wcc.Selectors( + label="Color map", + children=[ + wcc.Dropdown( + id=self.register_component_unique_id( + ElementIds.Settings.COLORMAP + ), + options=list_to_options(self.colormaps), + value=self.colormaps[0], + clearable=False, + ) + ], + ), + wcc.Checklist( + id=self.register_component_unique_id(ElementIds.Settings.SHOW_CUBEAXES), + options=["Show bounding box"], + value=["Show bounding box"], + ), + ] diff --git a/webviz_subsurface/plugins/_grid_viewer/views/view_3d/view_elements/_vtk_view_3d_element.py b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/view_elements/_vtk_view_3d_element.py new file mode 100644 index 000000000..49b0f8834 --- /dev/null +++ b/webviz_subsurface/plugins/_grid_viewer/views/view_3d/view_elements/_vtk_view_3d_element.py @@ -0,0 +1,133 @@ +from typing import List, Tuple + +from dash.development.base_component import Component +from webviz_config.webviz_plugin_subclasses import ViewElementABC + +import webviz_vtk +from dash import callback, Input, Output +from webviz_subsurface.plugins._grid_viewer._layout_elements import ElementIds +from ..settings import Settings + + +class VTKView3D(ViewElementABC): + def __init__(self) -> None: + super().__init__() + self.add_settings_group(Settings(), ElementIds.Settings.ID) + + def inner_layout(self) -> Component: + + return webviz_vtk.View( + id=self.register_component_unique_id(ElementIds.VTKVIEW3D.VIEW), + style={"height": "90vh"}, + pickingModes=["click"], + autoResetCamera=True, + interactorSettings=[ + { + "button": 1, + "action": "Zoom", + "scrollEnabled": True, + }, + { + "button": 3, + "action": "Pan", + }, + { + "button": 2, + "action": "Rotate", + }, + { + "button": 1, + "action": "Pan", + "shift": True, + }, + { + "button": 1, + "action": "Zoom", + "alt": True, + }, + { + "button": 1, + "action": "Roll", + "alt": True, + "shift": True, + }, + ], + children=[ + webviz_vtk.GeometryRepresentation( + id=self.register_component_unique_id( + ElementIds.VTKVIEW3D.GRID_REPRESENTATION + ), + showCubeAxes=True, + showScalarBar=True, + children=[ + webviz_vtk.PolyData( + id=self.register_component_unique_id( + ElementIds.VTKVIEW3D.GRID_POLYDATA + ), + children=[ + webviz_vtk.CellData( + [ + webviz_vtk.DataArray( + id=self.register_component_unique_id( + ElementIds.VTKVIEW3D.GRID_CELLDATA + ), + registration="setScalars", + name="scalar", + ) + ] + ) + ], + ) + ], + property={"edgeVisibility": True}, + ), + webviz_vtk.GeometryRepresentation( + id=self.register_component_unique_id( + ElementIds.VTKVIEW3D.PICK_REPRESENTATION + ), + actor={"visibility": False}, + children=[ + webviz_vtk.Algorithm( + id=self.register_component_unique_id( + ElementIds.VTKVIEW3D.PICK_SPHERE + ), + vtkClass="vtkSphereSource", + ) + ], + ), + ], + ) + + def set_callbacks(self) -> None: + @callback( + Output( + self.component_unique_id( + ElementIds.VTKVIEW3D.GRID_REPRESENTATION + ).to_string(), + "actor", + ), + Output( + self.component_unique_id( + ElementIds.VTKVIEW3D.GRID_REPRESENTATION + ).to_string(), + "showCubeAxes", + ), + Input( + self.settings_groups()[0] + .component_unique_id(ElementIds.Settings.ZSCALE) + .to_string(), + "value", + ), + Input( + self.settings_groups()[0] + .component_unique_id(ElementIds.Settings.SHOW_CUBEAXES) + .to_string(), + "value", + ), + ) + def _set_representation_actor( + z_scale: int, axes_is_on: List[str] + ) -> Tuple[dict, bool]: + show_axes = bool(z_scale == 1 and axes_is_on) + actor = {"scale": (1, 1, z_scale)} + return actor, show_axes