Add get_simbox_params

Author: oddvarlia

docs/get_rms_simbox.rst

@@ -0,0 +1,30 @@
+rms.get_rms_simbox
+============================
+
+The function extract the information from a geomodel grid necessary to
+create a box grid representing the RMS simulation box. The purpose is to save the necessary
+information to either create a box grid with same resolution, size, location and orientation
+as the a xy-regular geogrid. This information can be used in ERT to calculate cell center points
+for field parameters where (x,y) coordinates is relevant for the field parameter like in
+distance-based localization.
+
+Example of use
+----------------
+
+The function is called from a python job in RMS and use the rmsapi to get the information.
+
+Example from Drogon where a small yaml file is written for each of the three zones.
+Here zone_index is counted from 0 from top of a (multizone-) geogrid.
+Drogon has only single zone grids, therefore zone_index = 0 for all three grids.
+
+.. code-block:: python
+
+    from fmu.tools.rms.get_rms_simbox_params import get_simbox_param, write_simbox
+
+    if __name__ == "__main__":
+        grid_model_names = ["Geogrid_Valysar", "Geogrid_Therys", "Geogrid_Volon"]
+        for grid_model_name in grid_model_names:
+            filename = "simbox_" + grid_model_name + ".yml"
+            zone_index = 0
+            simbox = get_simbox_param(project, grid_model_name, zone_index)
+            write_simbox(filename, simbox)

docs/index.rst

@@ -21,6 +21,7 @@ Contents:
    update_petro_real
    export_and_import_fields
    copy_rms_param_to_ertbox_grid
+   get_rms_simbox
    utilities
    rmsvolumetrics2csv
    ensembles

src/fmu/tools/rms/get_rms_simbox_params.py

@@ -0,0 +1,67 @@
+"""
+Extract simbox information necessary to be able to both create a box grid
+with the simbox size, resolution and orientation and as input to calculate
+cell center position of each grid cell in simbox grid. Note, only that depth
+of the simbox is irrelevant where this is meant to be used in distance-based
+localization and therefore
+not saved.
+"""
+# TODO: NOTE: A bug in RMS14.2.2 cause the origin to be wrong.
+# Have to verify that newer versions of RMS is corrected before the function
+# below can be used. The RMS bug is only for rotated grids not for
+# rotation = 0. The error in RMS is that origin is set equal to MinX, MinY
+# and not the local origin for the rotated coordinate system aligned with
+# the rotated grid, typically the lower left corner point of the grid.
+
+from typing import Any
+
+import yaml
+
+
+def get_simbox_param(project: Any, grid_model_name: str, zone_index: int) -> dict:
+    if grid_model_name in project.grid_models:
+        grid_model = project.grid_models[grid_model_name]
+        grid3D = grid_model.get_grid()
+        cell_increments = grid3D.simbox.cell_increments
+        simbox_indexer = grid3D.simbox_indexer
+        ijk_handedness = simbox_indexer.ijk_handedness
+        zonation = simbox_indexer.zonation
+        layer_range = list(zonation[zone_index][0])
+        nlayers = len(layer_range)
+        dimensions = simbox_indexer.dimensions
+        zone_name = grid3D.zone_names[zone_index]
+        simbox = {
+            "name": str(zone_name),
+            "origin": [float(grid3D.origin[0]), float(grid3D.origin[1])],
+            "rotation": float(grid3D.rotation),
+        }
+        simbox["size"] = [
+            float(cell_increments["x_increment"] * dimensions[0]),
+            float(cell_increments["y_increment"] * dimensions[1]),
+            float(cell_increments["z_increments"][zone_index] * nlayers),
+        ]
+        simbox["dimensions"] = [
+            int(dimensions[0]),
+            int(dimensions[1]),
+            int(nlayers),
+        ]
+        simbox["handedness"] = str(ijk_handedness)
+        return simbox
+    raise ValueError(f"Unknown grid model {grid_model_name}")
+
+
+def write_simbox(filename: str, simbox: dict) -> None:
+    simbox_data_yml = yaml.dump(simbox, default_flow_style=False, sort_keys=False)
+    print(f"Write file:  {filename}")
+    with open(filename, "w") as file:
+        file.write(simbox_data_yml)
+        file.write("\n")
+
+
+# if __name__ == "__main__":
+#    grid_model_names = ["Geogrid_Valysar", "Geogrid_Therys", "Geogrid_Volon"]
+#    for grid_model_name in grid_model_names:
+#        filename = "tmp_simbox_" + grid_model_name + ".yml"
+#        zone_index = 0
+#        simbox = get_simbox_param(project, grid_model_name, zone_index)
+#        write_simbox(filename, simbox)

tests/rms/rms_simbox/ref_simbox_angle0_zone0.txt

@@ -0,0 +1,15 @@
+name: ZoneA
+origin:
+- 10000.0
+- 20000.0
+rotation: 0.0
+size:
+- 2500.0
+- 3500.0
+- 50.0
+dimensions:
+- 50
+- 70
+- 10
+handedness: left
+

tests/rms/rms_simbox/ref_simbox_angle0_zone1.txt

@@ -0,0 +1,15 @@
+name: ZoneB
+origin:
+- 10000.0
+- 20000.0
+rotation: 0.0
+size:
+- 2500.0
+- 3500.0
+- 175.0
+dimensions:
+- 50
+- 70
+- 35
+handedness: left
+

tests/rms/test_get_rms_simbox_params.py

@@ -0,0 +1,159 @@
+import contextlib
+import filecmp
+import shutil
+from os.path import isdir
+from pathlib import Path
+from typing import Any
+
+import pytest
+
+with contextlib.suppress(ImportError):
+    import rmsapi
+
+import xtgeo
+
+from fmu.tools.rms.get_rms_simbox_params import get_simbox_param, write_simbox
+
+# ======================================================================================
+# settings to create RMS project!
+
+REMOVE_RMS_PROJECT_AFTER_TEST = False
+
+TMPD = Path("TMP")
+TMPD.mkdir(parents=True, exist_ok=True)
+
+
+PROJNAME = "tmp_project_get_rms_simbox_params.rmsxxx"
+PRJ = str(TMPD / PROJNAME)
+RESULTDIR = TMPD / "rms_simbox"
+RESULTDIR.mkdir(parents=True, exist_ok=True)
+REFERENCE_DIR = Path("tests/rms/rms_simbox")
+
+NX = 50
+NY = 70
+NZ_ZONEA = 10
+NZ_ZONEB = 35
+XINC = 50.0
+YINC = 50.0
+ZINC = 5.0
+ORIGIN = (10000.0, 20000.0, 2000.0)
+ROTATION = 50.0
+ZONEA = "ZoneA"
+ZONEB = "ZoneB"
+
+
+ZONEA_NUMBER = 1
+ZONEB_NUMBER = 2
+
+
+GRID_MODEL_NAME = "Geogrid"
+OUTPUT_FILE_PREFIX = "simbox"
+OUTPUT_REF_FILE_PREFIX = "ref_simbox"
+
+
+@pytest.mark.skipunlessroxar
+@pytest.mark.parametrize(
+    "rotation, flip, zone_index_list",
+    [
+        (
+            0.0,
+            1,
+            [0, 1],
+        ),
+        # TODO  Extend with more cases with rotation when RMS bug is fixed
+    ],
+)
+def test_get_rms_simbox_params(
+    rotation: float, flip: int, zone_index_list: list[int]
+) -> None:
+    """Create a tmp RMS project for testing, populate with basic data."""
+    project = create_project()
+
+    create_grids(project, rotation, flip)
+    for zone_index in zone_index_list:
+        simbox_output_file_name = Path(RESULTDIR) / Path(
+            OUTPUT_FILE_PREFIX
+            + "_angle"
+            + str(int(rotation))
+            + "_zone"
+            + str(zone_index)
+            + ".txt"
+        )
+        simbox_dict = get_simbox_param(project, GRID_MODEL_NAME, zone_index)
+        write_simbox(simbox_output_file_name, simbox_dict)
+
+        # Verify that original is equal to the new params
+        reference_filename = Path(REFERENCE_DIR) / Path(
+            OUTPUT_REF_FILE_PREFIX
+            + "_angle"
+            + str(int(rotation))
+            + "_zone"
+            + str(zone_index)
+            + ".txt"
+        )
+        compare_results_with_reference(simbox_output_file_name, reference_filename)
+
+    project.close()
+
+    if REMOVE_RMS_PROJECT_AFTER_TEST:
+        print("\n******* Teardown RMS project!\n")
+        if isdir(PRJ):
+            print("Remove existing project!")
+            shutil.rmtree(PRJ)
+        if isdir(RESULTDIR):
+            print("Remove temporary files")
+            shutil.rmtree(RESULTDIR)
+
+
+def create_project():
+    """Create a tmp RMS project for testing, populate with basic data."""
+
+    prj1 = str(PRJ)
+
+    print("\n******** Setup RMS project!\n")
+    if isdir(prj1):
+        print("Remove existing project! (1)")
+        shutil.rmtree(prj1)
+
+    project = rmsapi.Project.create()
+
+    rox = xtgeo.RoxUtils(project)
+    print("rmsapi version is", rox.roxversion)
+    print("RMS version is", rox.rmsversion(rox.roxversion))
+    assert "1." in rox.roxversion
+
+    project.save_as(prj1)
+    return project
+
+
+def create_grids(project: Any, rotation: float, flip: int = 1):
+    nx = NX
+    ny = NY
+    nz_zoneA = NZ_ZONEA
+    nz_zoneB = NZ_ZONEB
+    nz = nz_zoneA + nz_zoneB
+    dimension = (nx, ny, nz)
+    increment = (XINC, YINC, ZINC)
+    origin = ORIGIN
+
+    subgrid_dict = {
+        ZONEA: nz_zoneA,
+        ZONEB: nz_zoneB,
+    }
+    geogrid = xtgeo.create_box_grid(
+        dimension,
+        origin=origin,
+        increment=increment,
+        rotation=rotation,
+        flip=flip,
+    )
+    geogrid.set_subgrids(subgrid_dict)
+    geogrid.to_roxar(project, GRID_MODEL_NAME)
+
+
+def compare_results_with_reference(filename, reference_filename):
+    # Compare with reference
+    check = filecmp.cmp(filename, reference_filename)
+    if check:
+        print("Check OK for multi zone grid")
+    assert check