Calculate and store ERTBOX parameters

ERTBOX is a rectangular grid created in RMS that envelopes a field.
From it, we calculate all parameters needed to do distance based localization.
Since more information is written to storage, we do a storage migration and test that it works.

Author: dafeda

src/ert/config/ensemble_config.py

@@ -2,12 +2,11 @@
 
 import logging
 from collections import Counter
+from pathlib import Path
 from typing import Self
 
 from pydantic import BaseModel, Field, model_validator
 
-from ert.field_utils import get_shape
-
 from .ext_param_config import ExtParamConfig
 from .field import Field as FieldConfig
 from .gen_data_config import GenDataConfig
@@ -92,18 +91,13 @@ def from_dict(cls, config_dict: ConfigDict) -> EnsembleConfig:
         gen_kw_list = config_dict.get(ConfigKeys.GEN_KW, [])
         surface_list = config_dict.get(ConfigKeys.SURFACE, [])
         field_list = config_dict.get(ConfigKeys.FIELD, [])
-        global_dims = None
 
-        # When users specify GRID as a separate line in the config,
-        # and not as an option to the FIELD keyword.
         if global_grid_file_path is not None:
-            try:
-                global_dims = get_shape(global_grid_file_path)
-            except Exception as err:
-                raise ConfigValidationError.with_context(
-                    f"Could not read grid file {global_grid_file_path}: {err}",
-                    global_grid_file_path,
-                ) from err
+            global_grid_file_path = Path(global_grid_file_path)
+
+            grid_extension = global_grid_file_path.suffix.lower()
+            if grid_extension not in {".egrid", ".grid"}:
+                raise ConfigValidationError("Only EGRID and GRID formats are supported")
 
         def make_field(field_list: list[str | dict[str, str]]) -> FieldConfig:
             # An example of `field_list` when the keyword `GRID` is set:
@@ -124,26 +118,10 @@ def make_field(field_list: list[str | dict[str, str]]) -> FieldConfig:
 
             # Use field-specific grid if provided,
             # otherwise fall back to global grid.
-            if grid_file_path is not None:
-                try:
-                    dims = get_shape(grid_file_path)
-                except Exception as err:
-                    raise ConfigValidationError.with_context(
-                        f"Could not read grid file {grid_file_path}: {err}",
-                        grid_file_path,
-                    ) from err
-            else:
-                grid_file_path = global_grid_file_path
-                dims = global_dims
-
-            if dims is None:
-                raise ConfigValidationError.with_context(
-                    f"Grid file {grid_file_path} did not contain dimensions",
-                    grid_file_path,
-                )
-            assert grid_file_path is not None
+            if grid_file_path is None:
+                grid_file_path = str(global_grid_file_path)
 
-            return FieldConfig.from_config_list(grid_file_path, dims, field_list)
+            return FieldConfig.from_config_list(grid_file_path, field_list)
 
         gen_kw_cfgs = [
             cfg for g in gen_kw_list for cfg in GenKwConfig.from_config_list(g)

src/ert/config/field.py

@@ -11,9 +11,19 @@
 import networkx as nx
 import numpy as np
 import xarray as xr
+import xtgeo  # type: ignore
 from pydantic import field_serializer
 
-from ert.field_utils import FieldFileFormat, Shape, read_field, read_mask, save_field
+from ert.field_utils import (
+    ErtboxParameters,
+    FieldFileFormat,
+    Shape,
+    calculate_ertbox_parameters,
+    get_shape,
+    read_field,
+    read_mask,
+    save_field,
+)
 from ert.substitutions import substitute_runpath_name
 from ert.utils import log_duration
 
@@ -83,9 +93,7 @@ def adjust_graph_for_masking(
 
 class Field(ParameterConfig):
     type: Literal["field"] = "field"
-    nx: int
-    ny: int
-    nz: int
+    ertbox_params: ErtboxParameters
     file_format: FieldFileFormat
     output_transformation: str | None
     input_transformation: str | None
@@ -115,20 +123,14 @@ def metadata(self) -> list[ParameterMetadata]:
                 key=self.name,
                 transformation=self.output_transformation,
                 dimensionality=3,
-                userdata={
-                    "data_origin": "FIELD",
-                    "nx": self.nx,
-                    "ny": self.ny,
-                    "nz": self.nz,
-                },
+                userdata={"data_origin": "FIELD", "ertbox_params": self.ertbox_params},
             )
         ]
 
     @classmethod
     def from_config_list(
         cls,
         grid_file_path: str,
-        dims: Shape,
         config_list: list[str | dict[str, str]],
     ) -> Self:
         name = cast(str, config_list[0])
@@ -198,11 +200,32 @@ def from_config_list(
         assert file_format is not None
 
         assert init_files is not None
+
+        grid_extension = Path(grid_file_path).suffix.lower()
+
+        try:
+            if grid_extension == ".egrid":
+                grid = xtgeo.grid_from_file(grid_file_path)
+                ertbox_params = calculate_ertbox_parameters(grid)
+            else:
+                dims = get_shape(grid_file_path)
+
+                if dims is None:
+                    raise ConfigValidationError.with_context(
+                        f"Grid file {grid_file_path} did not contain dimensions",
+                        grid_file_path,
+                    )
+
+                ertbox_params = ErtboxParameters(dims.nx, dims.ny, dims.nz)
+        except Exception as err:
+            raise ConfigValidationError.with_context(
+                f"Could not read grid file {grid_file_path}: {err}",
+                grid_file_path,
+            ) from err
+
         return cls(
             name=name,
-            nx=dims.nx,
-            ny=dims.ny,
-            nz=dims.nz,
+            ertbox_params=ertbox_params,
             file_format=file_format,
             output_transformation=output_transform,
             input_transformation=init_transform,
@@ -217,7 +240,7 @@ def from_config_list(
 
     def __len__(self) -> int:
         if self.mask_file is None:
-            return self.nx * self.ny * self.nz
+            return self.ertbox_params.nx * self.ertbox_params.ny * self.ertbox_params.nz
 
         # Uses int() to convert to standard python int for mypy
         return int(np.size(self.mask) - np.count_nonzero(self.mask))
@@ -236,7 +259,11 @@ def read_from_runpath(
                             run_path / file_name,
                             self.name,
                             self.mask,
-                            Shape(self.nx, self.ny, self.nz),
+                            Shape(
+                                self.ertbox_params.nx,
+                                self.ertbox_params.ny,
+                                self.ertbox_params.nz,
+                            ),
                         ),
                         self.input_transformation,
                     ),
@@ -332,10 +359,22 @@ def mask(self) -> Any:
 
     def load_parameter_graph(self) -> nx.Graph:  # type: ignore
         parameter_graph = create_flattened_cube_graph(
-            px=self.nx, py=self.ny, pz=self.nz
+            px=self.ertbox_params.nx, py=self.ertbox_params.ny, pz=self.ertbox_params.nz
         )
         return adjust_graph_for_masking(G=parameter_graph, mask=self.mask.flatten())
 
+    @property
+    def nx(self) -> int:
+        return self.ertbox_params.nx
+
+    @property
+    def ny(self) -> int:
+        return self.ertbox_params.ny
+
+    @property
+    def nz(self) -> int:
+        return self.ertbox_params.nz
+
 
 TRANSFORM_FUNCTIONS = {
     "LN": np.log,

src/ert/field_utils/__init__.py

@@ -1,11 +1,21 @@
 from __future__ import annotations
 
 from .field_file_format import FieldFileFormat
-from .field_utils import Shape, get_shape, read_field, read_mask, save_field
+from .field_utils import (
+    ErtboxParameters,
+    Shape,
+    calculate_ertbox_parameters,
+    get_shape,
+    read_field,
+    read_mask,
+    save_field,
+)
 
 __all__ = [
+    "ErtboxParameters",
     "FieldFileFormat",
     "Shape",
+    "calculate_ertbox_parameters",
     "get_shape",
     "read_field",
     "read_mask",

src/ert/field_utils/field_utils.py

@@ -1,18 +1,21 @@
 from __future__ import annotations
 
+import math
 import os
 from pathlib import Path
 from typing import TYPE_CHECKING, Any, NamedTuple, TypeAlias
 
 import numpy as np
 import resfo
+from pydantic.dataclasses import dataclass
 
 from .field_file_format import ROFF_FORMATS, FieldFileFormat
 from .grdecl_io import export_grdecl, import_bgrdecl, import_grdecl
 from .roff_io import export_roff, import_roff
 
 if TYPE_CHECKING:
     import numpy.typing as npt
+    import xtgeo  # type: ignore
 
 _PathLike: TypeAlias = str | os.PathLike[str]
 
@@ -96,6 +99,116 @@ def get_shape(
     return shape
 
 
+@dataclass(frozen=True)
+class ErtboxParameters:
+    nx: int
+    ny: int
+    nz: int
+    xlength: float | None = None
+    ylength: float | None = None
+    xinc: float | None = None
+    yinc: float | None = None
+    rotation_angle: float | None = None
+    origin: tuple[float, float] | None = None
+
+
+def calculate_ertbox_parameters(
+    grid: xtgeo.Grid, left_handed: bool = False
+) -> ErtboxParameters:
+    """Calculate ERTBOX grid parameters from an XTGeo grid.
+
+    Extracts geometric parameters including dimensions, cell increments,
+    rotation angle, and origin coordinates needed for ERTBOX.
+
+    Args:
+        grid: XTGeo Grid3D object
+        left_handed: If True, use left-handed coordinate system (default: False)
+
+    Returns:
+        ErtboxParameters with grid dimensions, increments, rotation, and origin
+    """
+
+    (nx, ny, nz) = grid.dimensions
+
+    corner_indices = []
+
+    if left_handed:
+        origin_cell = (1, 1, 1)
+        x_direction_cell = (nx, 1, 1)
+        y_direction_cell = (1, ny, 1)
+    else:
+        origin_cell = (1, ny, 1)
+        x_direction_cell = (nx, ny, 1)
+        y_direction_cell = (1, 1, 1)
+
+    corner_indices = [origin_cell, x_direction_cell, y_direction_cell]
+
+    # List with 3 elements, where each element contains the coordinates
+    # for all 8 corners of a single grid cell.
+    coord_cell = []
+
+    for corner_index in corner_indices:
+        # Get real-world (x,y,z) coordinates for all 8 corners of this grid cell
+        # Returns 24 values: [x0,y0,z0, x1,y1,z1, ..., x7,y7,z7]
+        coord = grid.get_xyz_cell_corners(ijk=corner_index, activeonly=False)
+        coord_cell.append(coord)
+
+    if left_handed:
+        # Origin: cell (1,1,1), corner 0
+        x0 = coord_cell[0][0]
+        y0 = coord_cell[0][1]
+
+        # X-direction: cell (nx,1,1), corner 1
+        x1 = coord_cell[1][3]
+        y1 = coord_cell[1][4]
+
+        # Y-direction: cell (1,ny,1), corner 2
+        x2 = coord_cell[2][6]
+        y2 = coord_cell[2][7]
+    else:
+        # Origin: cell (1,ny,1), corner 2
+        x0 = coord_cell[0][6]
+        y0 = coord_cell[0][7]
+
+        # X-direction: cell (nx,ny,1), corner 3
+        x1 = coord_cell[1][9]
+        y1 = coord_cell[1][10]
+
+        # Y-direction: cell (1,1,1), corner 0
+        x2 = coord_cell[2][0]
+        y2 = coord_cell[2][1]
+
+    deltax1 = x1 - x0
+    deltay1 = y1 - y0
+
+    deltax2 = x2 - x0
+    deltay2 = y2 - y0
+
+    xlength = math.sqrt(deltax1**2 + deltay1**2)
+    ylength = math.sqrt(deltax2**2 + deltay2**2)
+    xinc = xlength / nx
+    yinc = ylength / ny
+
+    if math.fabs(deltax1) < 0.00001:
+        angle = 90.0 if deltay1 > 0 else -90.0
+    elif deltax1 > 0:
+        angle = math.atan(deltay1 / deltax1) * 180.0 / math.pi
+    elif deltax1 < 0:
+        angle = (math.atan(deltay1 / deltax1) + math.pi) * 180.0 / math.pi
+
+    return ErtboxParameters(
+        nx=nx,
+        ny=ny,
+        nz=nz,
+        xlength=xlength,
+        ylength=ylength,
+        xinc=xinc,
+        yinc=yinc,
+        rotation_angle=angle,
+        origin=(x0, y0),
+    )
+
+
 def read_field(
     field_path: _PathLike,
     field_name: str,

src/ert/gui/tools/plot/plot_window.py

@@ -247,7 +247,7 @@ def updatePlot(self, layer: int | None = None) -> None:
             if "FIELD" in key_def.metadata["data_origin"]:
                 plot_widget.showLayerWidget.emit(True)
 
-                layers = key_def.metadata["nz"]
+                layers = key_def.metadata["ertbox_params"]["nz"]
                 plot_widget.updateLayerWidget.emit(layers)
 
                 if layer is None:

src/ert/storage/local_storage.py

@@ -30,7 +30,7 @@
 
 logger = logging.getLogger(__name__)
 
-_LOCAL_STORAGE_VERSION = 13
+_LOCAL_STORAGE_VERSION = 14
 
 
 class _Migrations(BaseModel):
@@ -493,6 +493,7 @@ def _migrate(self, version: int) -> None:
             to11,
             to12,
             to13,
+            to14,
         )
 
         try:
@@ -535,6 +536,7 @@ def _migrate(self, version: int) -> None:
                     10: to11,
                     11: to12,
                     12: to13,
+                    13: to14,
                 }
                 for from_version in range(version, _LOCAL_STORAGE_VERSION):
                     migrations[from_version].migrate(self.path)