Added distance-based localization functions for 2D fields

added snapshot test
Functions:
- calc_rho_for_2d_grid_layer  to calculate elements of the RHO matrix.

Author: oddvarlia

src/ert/config/surface_config.py

@@ -10,6 +10,11 @@
 from pydantic import field_serializer
 from surfio import IrapHeader, IrapSurface
 
+from ert.field_utils import (
+    calc_rho_for_2d_grid_layer,
+    transform_local_ellipse_angle_to_local_coords,
+    transform_positions_to_local_field_coordinates,
+)
 from ert.substitutions import substitute_runpath_name
 
 from ._str_to_bool import str_to_bool
@@ -248,3 +253,55 @@ def load_parameter_graph(self) -> nx.Graph[int]:
         this flattening process"""
 
         return create_flattened_cube_graph(px=self.ncol, py=self.nrow, pz=1)
+
+    def calc_rho_for_2d_grid_layer(
+        self,
+        obs_xpos: npt.NDArray[np.float64],
+        obs_ypos: npt.NDArray[np.float64],
+        obs_main_range: npt.NDArray[np.float64],
+        obs_perp_range: npt.NDArray[np.float64],
+        obs_anisotropy_angle: npt.NDArray[np.float64],
+    ) -> npt.NDArray[np.float64]:
+        """Function to calculate scaling values to be used in the RHO matrix
+        for distance-based localization.
+
+        Args:
+            obs_xpos: x-coordinates in global coordinates of observations
+            obs_ypos: y-coordinates in global coordinates of observations
+            obs_main_range: Size of influence ellipse main principal direction.
+            obs_perp_range: Size of influence ellipse second principal direction.
+            obs_anisotropy_angle: Rotation angle anticlock wise of main principal
+              direction of influence ellipse relative to global coordinate
+              system's x-axis.
+
+        Returns:
+            Scaling values (elements of the RHO matrix) as a numpy array
+              of shape=(nx,ny,nobservations)
+
+        """
+        # Transform observation positions to local surface coordinates
+        xpos, ypos = transform_positions_to_local_field_coordinates(
+            (self.xori, self.yori), self.rotation, obs_xpos, obs_ypos
+        )
+        # Transform ellipse orientation to local surface coordinates
+        rotation_angle_of_localization_ellipse = (
+            transform_local_ellipse_angle_to_local_coords(
+                self.rotation, obs_anisotropy_angle
+            )
+        )
+
+        # Assume the coordinate system is not flipped.
+        # This means the right_handed_grid_indexing is False
+        assert self.yflip == 1
+        return calc_rho_for_2d_grid_layer(
+            self.ncol,
+            self.nrow,
+            self.xinc,
+            self.yinc,
+            xpos,
+            ypos,
+            obs_main_range,
+            obs_perp_range,
+            rotation_angle_of_localization_ellipse,
+            right_handed_grid_indexing=False,
+        )

tests/ert/unit_tests/config/snapshots/test_surface_config/test_calc_rho_for_2d_grid_layer/xpos0-ypos0-main_range0-perp_range0-anisotropy_angle0/testdata_rho_for_one_grid_layer.txt

@@ -0,0 +1,49 @@
+[[[0.50498824 0.00651928 0.         ... 0.99997919 0.         0.        ]
+  [0.54097227 0.00890534 0.         ... 0.99786346 0.         0.        ]
+  [0.57488894 0.01182687 0.         ... 0.99172024 0.         0.        ]
+  ...
+  [0.         0.         0.06072022 ... 0.         0.         0.99464177]
+  [0.         0.         0.05473404 ... 0.         0.         0.99848823]
+  [0.         0.         0.04896573 ... 0.         0.         0.99964644]]
+
+ [[0.54097227 0.00602875 0.         ... 0.99786346 0.         0.        ]
+  [0.57926199 0.00832618 0.         ... 0.99981325 0.         0.        ]
+  [0.61536138 0.01116308 0.         ... 0.99761887 0.         0.        ]
+  ...
+  [0.         0.         0.06174645 ... 0.         0.         0.99613059]
+  [0.         0.         0.0554226  ... 0.         0.         0.99685459]
+  [0.         0.         0.04936423 ... 0.         0.         0.9948815 ]]
+
+ [[0.57488894 0.00553249 0.         ... 0.99172024 0.         0.        ]
+  [0.61536138 0.00773125 0.         ... 0.99761887 0.         0.        ]
+  [0.65353538 0.01047068 0.         ... 0.99948267 0.         0.        ]
+  ...
+  [0.         0.         0.06243868 ... 0.         0.         0.99136881]
+  [0.         0.         0.05580123 ... 0.         0.         0.98902204]
+  [0.         0.         0.04947859 ... 0.         0.         0.98407544]]
+
+ ...
+
+ [[0.         0.         0.04947859 ... 0.         0.46753824 0.        ]
+  [0.         0.         0.05580123 ... 0.         0.82297363 0.        ]
+  [0.         0.         0.06243868 ... 0.         0.99434669 0.        ]
+  ...
+  [0.         0.01047068 0.         ... 0.18232878 0.         0.        ]
+  [0.         0.00773125 0.         ... 0.16524114 0.         0.        ]
+  [0.         0.00553249 0.         ... 0.14890707 0.         0.        ]]
+
+ [[0.         0.         0.04936423 ... 0.         0.82503208 0.        ]
+  [0.         0.         0.0554226  ... 0.         0.99794506 0.        ]
+  [0.         0.         0.06174645 ... 0.         0.82297363 0.        ]
+  ...
+  [0.         0.01116308 0.         ... 0.19644278 0.         0.        ]
+  [0.         0.00832618 0.         ... 0.17895554 0.         0.        ]
+  [0.         0.00602875 0.         ... 0.16213323 0.         0.        ]]
+
+ [[0.         0.         0.04896573 ... 0.         0.99976955 0.        ]
+  [0.         0.         0.05473404 ... 0.         0.82503208 0.        ]
+  [0.         0.         0.06072022 ... 0.         0.46753824 0.        ]
+  ...
+  [0.         0.01182687 0.         ... 0.21060644 0.         0.        ]
+  [0.         0.00890534 0.         ... 0.19282388 0.         0.        ]
+  [0.         0.00651928 0.         ... 0.17560682 0.         0.        ]]]

tests/ert/unit_tests/config/test_surface_config.py

@@ -4,7 +4,7 @@
 import numpy as np
 import pytest
 import xtgeo
-from surfio import IrapSurface
+from surfio import IrapHeader, IrapSurface
 
 from ert.config import ConfigValidationError, SurfaceConfig
 from ert.config.parameter_config import InvalidParameterFile
@@ -339,3 +339,105 @@ def test_surface_create_storage_datasets_raises_surface_mismatch_error_when_the_
     )
     with pytest.raises(InvalidParameterFile, match=expected_error_msg):
         next(storage_dataset_iterator)
+
+
+@pytest.fixture
+def surface_for_dl():
+    # Return a SurfaceConfig object for test purpose
+
+    # Create a synthetic surface grid with rotation
+    nx = 100
+    ny = 120
+    xori = 1000.0
+    yori = 2000.0
+    xsize = 500.0
+    ysize = 600.0
+    xinc = xsize / nx
+    yinc = ysize / ny
+    rotation = 0.0
+    return SurfaceConfig(
+        type="surface",
+        name="MySurface",
+        forward_init=True,
+        update=True,
+        ncol=nx,
+        nrow=ny,
+        xori=xori,
+        yori=yori,
+        xinc=xinc,
+        yinc=yinc,
+        rotation=rotation,
+        yflip=1,
+        forward_init_file="dummy.txt",
+        output_file=Path("dummy.txt"),
+        base_surface_path="dummy.txt",
+    )
+
+
+@pytest.mark.parametrize(
+    "xpos, ypos, main_range, perp_range, anisotropy_angle",
+    [
+        (
+            [1050.0, 1250.0, 1250.0, 1250.0, 1000.0, 1500.0, 1000.0],  # xpos
+            [2050.0, 2300.0, 2300.0, 2300.0, 2000.0, 2000.0, 2600.0],  # ypos
+            [100.0, 300.0, 300.0, 600.0, 1000.0, 300.0, 300.0],  # main_range
+            [100.0, 100.0, 100.0, 200.0, 100.0, 10.0, 100.0],  # perp_range
+            [0.0, 35.0, 135.0, -135.0, 45.0, -45.0, -60.0],  # angle
+        ),
+    ],
+)
+def test_calc_rho_for_2d_grid_layer(
+    snapshot,
+    surface_for_dl,
+    xpos: list[float],
+    ypos: list[float],
+    main_range: list[float],
+    perp_range: list[float],
+    anisotropy_angle: list[float],
+):
+    write_surface_file = False
+
+    xposition = np.array(xpos)
+    yposition = np.array(ypos)
+    mainrange = np.array(main_range)
+    perprange = np.array(perp_range)
+    angles = np.array(anisotropy_angle)
+
+    #   Dimension of rho_for_one_grid_layer is (nx,ny,nobs)
+    rho_for_one_grid_layer = surface_for_dl.calc_rho_for_2d_grid_layer(
+        xposition,
+        yposition,
+        mainrange,
+        perprange,
+        angles,
+    )
+    # Ensure -0 and +0 will be 0
+    rho_for_one_grid_layer = np.where(
+        rho_for_one_grid_layer == 0, 0.0, rho_for_one_grid_layer
+    )
+
+    snapshot.assert_match(
+        str(rho_for_one_grid_layer) + "\n", "testdata_rho_for_one_grid_layer.txt"
+    )
+    if write_surface_file:
+        # Write surface of rho for visualization
+
+        for obs_indx in range(len(xposition)):
+            surf = IrapSurface(
+                header=IrapHeader(
+                    ncol=surface_for_dl.ncol,
+                    nrow=surface_for_dl.nrow,
+                    xori=surface_for_dl.xori,
+                    yori=surface_for_dl.yori,
+                    xinc=surface_for_dl.xinc,
+                    yinc=surface_for_dl.yinc,
+                    rot=surface_for_dl.rotation,
+                    xrot=surface_for_dl.xori,
+                    yrot=surface_for_dl.yori,
+                ),
+                values=rho_for_one_grid_layer[:, :, obs_indx],
+            )
+
+            file_path = Path("tmp_2d_rho_for_obs_" + str(obs_indx) + ".txt")
+            print(f"Write file: {file_path}")
+            surf.to_ascii_file(file_path)