From 9facaedab1cf74e015c77abff3e9b9f879ca6c52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Wed, 15 Jan 2025 08:08:23 +0100 Subject: [PATCH] ENH: add method get_property_between_surfaces() for Grid class --- mypy.ini | 2 +- src/lib/include/xtgeo/grid3d.hpp | 28 ++++++ src/lib/include/xtgeo/regsurf.hpp | 14 +++ src/lib/include/xtgeo/xtgeo.h | 34 ------- src/lib/src/grdcp3d_calc_xyz.c | 84 ----------------- src/lib/src/grdcp3d_midpoint.c | 54 ----------- src/lib/src/grid3d/grid.cpp | 54 +++++++++++ src/lib/src/grid3d/grid_surf_oper.cpp | 116 ++++++++++++++++++++++++ src/lib/src/regsurf/utilities.cpp | 72 +++++++++++++++ src/xtgeo/grid3d/_grid_etc1.py | 85 ++++++++++++++++- src/xtgeo/grid3d/grid.py | 25 +++++ tests/test_grid3d/test_grid.py | 46 +++++++++- tests/test_internal/test_cell_grid3d.py | 58 +++++++++++- tests/test_internal/test_regsurf.py | 95 ++++++++++++++++++- 14 files changed, 589 insertions(+), 178 deletions(-) delete mode 100644 src/lib/src/grdcp3d_calc_xyz.c delete mode 100644 src/lib/src/grdcp3d_midpoint.c create mode 100644 src/lib/src/grid3d/grid_surf_oper.cpp diff --git a/mypy.ini b/mypy.ini index 81c8eb2a6..fba1bf9ab 100644 --- a/mypy.ini +++ b/mypy.ini @@ -5,7 +5,7 @@ ignore_missing_imports = True strict_equality = True warn_redundant_casts = True warn_unused_configs = True -warn_unused_ignores = True +warn_unused_ignores = False exclude = tests|docs|examples|bin diff --git a/src/lib/include/xtgeo/grid3d.hpp b/src/lib/include/xtgeo/grid3d.hpp index 01d042cc1..df6383075 100644 --- a/src/lib/include/xtgeo/grid3d.hpp +++ b/src/lib/include/xtgeo/grid3d.hpp @@ -3,6 +3,7 @@ #include // should be included first according to pybind11 docs #include #include +#include #include #include @@ -63,6 +64,30 @@ get_depth_in_cell(const double x, const std::array &corners, int option); +py::array_t +grid_assign_value_between_surfaces(const size_t ncol, + const size_t nrow, + const size_t nlay, + const py::array_t &xmid, + const py::array_t &ymid, + const py::array_t &zmid, + const size_t top_ncol, + const size_t top_nrow, + const double top_xori, + const double top_yori, + const double top_xinc, + const double top_yinc, + const double top_rotation, + const py::array_t &top_values, + const size_t bot_ncol, + const size_t bot_nrow, + const double bot_xori, + const double bot_yori, + const double bot_xinc, + const double bot_yinc, + const double bot_rotation, + const py::array_t &bot_values); + inline void init(py::module &m) { @@ -83,6 +108,9 @@ init(py::module &m) "Determine if a XY point is inside a cell, top or base."); m_grid3d.def("get_depth_in_cell", &get_depth_in_cell, "Determine the interpolated cell face Z from XY, top or base."); + m_grid3d.def("grid_assign_value_between_surfaces", + &grid_assign_value_between_surfaces, + "Make a property that is one if cell center is between two surfaces."); } } // namespace xtgeo::grid3d diff --git a/src/lib/include/xtgeo/regsurf.hpp b/src/lib/include/xtgeo/regsurf.hpp index 3335dfead..1d86ef331 100644 --- a/src/lib/include/xtgeo/regsurf.hpp +++ b/src/lib/include/xtgeo/regsurf.hpp @@ -30,6 +30,18 @@ get_xy_from_ij(const size_t i, const size_t nrow, const double angle_deg); +double +get_z_from_xy(const double x, + const double y, + const double xori, + const double yori, + const double xinc, + const double yinc, + const size_t ncol, + const size_t nrow, + const double angle_deg, + const py::array_t &values); + std::tuple find_cell_range(const double xmin, const double xmax, @@ -80,6 +92,8 @@ init(py::module &m) "Get the outer corners of a regular surface."); m_regsurf.def("get_xy_from_ij", &get_xy_from_ij, "Get the XY coordinates from the grid indices."); + m_regsurf.def("get_z_from_xy", &get_z_from_xy, + "Get the Z value in a regsurf from a x y point."); m_regsurf.def("find_cell_range", &find_cell_range, "Find the range of regular surface 2D nodes within a box."); m_regsurf.def("sample_grid3d_layer", &sample_grid3d_layer, diff --git a/src/lib/include/xtgeo/xtgeo.h b/src/lib/include/xtgeo/xtgeo.h index b4ec72c0e..1b7778ae5 100644 --- a/src/lib/include/xtgeo/xtgeo.h +++ b/src/lib/include/xtgeo/xtgeo.h @@ -1059,26 +1059,6 @@ extern "C" long n_swig_np_dbl_inplace_v1, // ntot, metric m); - void grdcp3d_calc_xyz(long ncol, - long nrow, - long nlay, - - double *swig_np_dbl_in_v1, // *coordsv, - long n_swig_np_dbl_in_v1, // ncoord, - float *swig_np_flt_in_v1, // *zcornsv, - long n_swig_np_flt_in_v1, // nzcorn, - int *swig_np_int_in_v1, // *actnumsv, - long n_swig_np_int_in_v1, // nactnum, - - int option, - - double *swig_np_dbl_aout_v1, // xarr - long n_swig_np_dbl_aout_v1, - double *swig_np_dbl_aout_v2, // yarr - long n_swig_np_dbl_aout_v2, - double *swig_np_dbl_aout_v3, // zarr - long n_swig_np_dbl_aout_v3); - void grd3d_conv_roxapi_grid(int nx, int ny, int nz, @@ -1394,20 +1374,6 @@ extern "C" long n_swig_np_int_inplace_v1 // nactnum ); - void grdcp3d_midpoint(long i, - long j, - long k, - long ncol, - long nrow, - long nlay, - double *coordsv, - long ncoord, - float *zcornsv, - long nzcorn, - double *x, - double *y, - double *z); - void grd3d_midpoint(int i, int j, int k, diff --git a/src/lib/src/grdcp3d_calc_xyz.c b/src/lib/src/grdcp3d_calc_xyz.c deleted file mode 100644 index 7af6cb4b1..000000000 --- a/src/lib/src/grdcp3d_calc_xyz.c +++ /dev/null @@ -1,84 +0,0 @@ -/* - *************************************************************************************** - * - * NAME: - * grdcp3d_calc_xyz.c - * - * DESCRIPTION: - * Get X Y Z vector per cell - * - * ARGUMENTS: - * nx..nz grid dimensions - * coordsv Coordinates - * zcornsv ZCORN array (pointer) of input - * p_x_v .. p_z_v Return arrays for X Y Z - * option 0: use all cells, 1: make undef if ACTNUM=0 - * - * RETURNS: - * Void, update arrays - * - * TODO/ISSUES/BUGS: - * Make proper return codes - * - * LICENCE: - * cf. XTGeo LICENSE - *************************************************************************************** - */ -#include -#include "common.h" - -void -grdcp3d_calc_xyz(long ncol, - long nrow, - long nlay, - double *coordsv, - long ncoord, - float *zcornsv, - long nzcorn, - int *actnumsv, - long nact, - int option, - double *xarr, // xarr - long nxarr, - double *yarr, // yarr - long nyarr, - double *zarr, // zarr - long nzarr) -{ - long ntotv[4] = { nact, nxarr, nyarr, nzarr }; - if (x_verify_vectorlengths(ncol, nrow, nlay, ncoord, nzcorn, ntotv, 4, - XTGFORMAT2) != 0) { - throw_exception("Errors in array lengths checks in grdcp3d_calc_xyz"); - return; - } - - double xv, yv, zv; - for (long i = 0; i < ncol; i++) { - for (long j = 0; j < nrow; j++) { - for (long k = 0; k < nlay; k++) { - /* Offset (i,j,k) by 1 because starting from 0. */ - long ic = x_ijk2ic(i + 1, j + 1, k + 1, ncol, nrow, nlay, 0); - if (ic < 0) { - throw_exception("Loop through grid resulted in index outside " - "boundary in grdcp3d_calc_xyz"); - return; - } - - /* If we want to mask inactive cells */ - if (option == 1 && actnumsv[ic] == 0) { - xarr[ic] = UNDEF; - yarr[ic] = UNDEF; - zarr[ic] = UNDEF; - continue; - } - - grdcp3d_midpoint(i, j, k, ncol, nrow, nlay, coordsv, ncoord, zcornsv, - nzcorn, &xv, &yv, &zv); - - xarr[ic] = xv; - yarr[ic] = yv; - zarr[ic] = zv; - } - } - } -} diff --git a/src/lib/src/grdcp3d_midpoint.c b/src/lib/src/grdcp3d_midpoint.c deleted file mode 100644 index 7be46b240..000000000 --- a/src/lib/src/grdcp3d_midpoint.c +++ /dev/null @@ -1,54 +0,0 @@ -/* - *************************************************************************************** - * - * NAME: - * grdcp3d_midpoint.c - * - * DESCRIPTION: - * Find the midpoint is a specific cell (cell count IJK with base 1) - * - * ARGUMENTS: - * i, j, k i Cell IJK - * ncol, nrow, nlay i Dimensions - * coordsv i Coordinate vector (with numpy dimensions) - * zcornsv i ZCORN vector (with numpy dimensions) - * x, y, z o XYZ output - * - * RETURNS: - * Void - * - * TODO/ISSUES/BUGS: - * None - * - * LICENCE: - * CF XTGeo's LICENSE - *************************************************************************************** - */ -#include - -void -grdcp3d_midpoint(long i, - long j, - long k, - long ncol, - long nrow, - long nlay, - double *coordsv, - long ncoord, - float *zcornsv, - long nzcorn, - double *x, - double *y, - double *z) - -{ - double c[24]; - - /* Get the (x, y, z) coordinates of all 8 corners in `c` */ - grdcp3d_corners(i, j, k, ncol, nrow, nlay, coordsv, ncoord, zcornsv, nzcorn, c); - - /* Find centroid via arithmetic mean (average of each dimensional coord) */ - *x = 0.125 * (c[0] + c[3] + c[6] + c[9] + c[12] + c[15] + c[18] + c[21]); - *y = 0.125 * (c[1] + c[4] + c[7] + c[10] + c[13] + c[16] + c[19] + c[22]); - *z = 0.125 * (c[2] + c[5] + c[8] + c[11] + c[14] + c[17] + c[20] + c[23]); -} diff --git a/src/lib/src/grid3d/grid.cpp b/src/lib/src/grid3d/grid.cpp index 78ac7ea1c..88b07a65a 100644 --- a/src/lib/src/grid3d/grid.cpp +++ b/src/lib/src/grid3d/grid.cpp @@ -55,6 +55,60 @@ grid_cell_volumes(const size_t ncol, return cellvols; } +/* + * Get cell centers for a grid. + * + * @param ncol Grid dimensions ncol/nx + * @param nrow Grid dimensions nrow/ny + * @param nlay Grid dimensions nlay/nz + * @param coordsv Grid Z coordinates vector + * @param zcornsv Grid Z corners vector + * @param actumsv Active cells vector + * @param asmasked Process grid cells as masked (return NaN for inactive cells) + * @return Arrays with the X, Y, Z coordinates of the cell centers + */ +std::tuple, py::array_t, py::array_t> +grid_cell_centers(const size_t ncol, + const size_t nrow, + const size_t nlay, + const py::array_t &coordsv, + const py::array_t &zcornsv, + const py::array_t &actnumsv, + const bool asmasked = false) +{ + pybind11::array_t xmid({ ncol, nrow, nlay }); + pybind11::array_t ymid({ ncol, nrow, nlay }); + pybind11::array_t zmid({ ncol, nrow, nlay }); + auto xmid_ = xmid.mutable_unchecked<3>(); + auto ymid_ = ymid.mutable_unchecked<3>(); + auto zmid_ = zmid.mutable_unchecked<3>(); + auto actnumsv_ = actnumsv.unchecked<3>(); + + for (size_t i = 0; i < ncol; i++) { + for (size_t j = 0; j < nrow; j++) { + for (size_t k = 0; k < nlay; k++) { + size_t idx = i * nrow * nlay + j * nlay + k; + if (asmasked && actnumsv_(i, j, k) == 0) { + xmid_(i, j, k) = std::numeric_limits::quiet_NaN(); + ymid_(i, j, k) = std::numeric_limits::quiet_NaN(); + zmid_(i, j, k) = std::numeric_limits::quiet_NaN(); + continue; + } + auto cr = + grid3d::cell_corners(i, j, k, ncol, nrow, nlay, coordsv, zcornsv); + + xmid_(i, j, k) = 0.125 * (cr[0] + cr[3] + cr[6] + cr[9] + cr[12] + + cr[15] + cr[18] + cr[21]); + ymid_(i, j, k) = 0.125 * (cr[1] + cr[4] + cr[7] + cr[10] + cr[13] + + cr[16] + cr[19] + cr[22]); + zmid_(i, j, k) = 0.125 * (cr[2] + cr[5] + cr[8] + cr[11] + cr[14] + + cr[17] + cr[20] + cr[23]); + } + } + } + return std::make_tuple(xmid, ymid, zmid); +} + /* * Compute cell height above ffl (free fluid level), as input to water saturation. Will * return hbot, htop, hmid (bottom of cell, top of cell, midpoint), but compute method diff --git a/src/lib/src/grid3d/grid_surf_oper.cpp b/src/lib/src/grid3d/grid_surf_oper.cpp new file mode 100644 index 000000000..6c82b4761 --- /dev/null +++ b/src/lib/src/grid3d/grid_surf_oper.cpp @@ -0,0 +1,116 @@ +// File: grid_surf_oper.cpp focus on grid or gridproperty that are associated with +// and/or modified by surface(s) input. This file is part of the xtgeo library. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +namespace xtgeo::grid3d { + +/* + * Create a parameter that is 1 if the cell center is between the top and bottom + * surfaces + * + * @param ncol Grid dimensions ncol/nx + * @param nrow Grid dimensions nrow/ny + * @param nlay Grid dimensions nlay/nz + * @param xmid Grid X midpoints + * @param ymid Grid Y midpoints + * @param zmid Grid Z midpoints + * @param top_ncol Top surface ncol + * @param top_nrow Top surface nrow + * @param top_xori Top surface x origin + * @param top_yori Top surface y origin + * @param top_xinc Top surface x increment + * @param top_yinc Top surface y increment + * @param top_rotation Top surface rotation + * @param top_values Top surface values + * @param bot_ncol Bottom surface ncol + * @param bot_nrow Bottom surface nrow + * @param bot_xori Bottom surface x origin + * @param bot_yori Bottom surface y origin + * @param bot_xinc Bottom surface x increment + * @param bot_yinc Bottom surface y increment + * @param bot_rotation Bottom surface rotation + * @param bot_values Bottom surface values + * @return An int array containing 1 if the cell is between the surfaces, 0 otherwise + */ +py::array_t +grid_assign_value_between_surfaces(const size_t ncol, + const size_t nrow, + const size_t nlay, + const py::array_t &xmid, + const py::array_t &ymid, + const py::array_t &zmid, + const size_t top_ncol, + const size_t top_nrow, + const double top_xori, + const double top_yori, + const double top_xinc, + const double top_yinc, + const double top_rotation, + const py::array_t &top_values, + const size_t bot_ncol, + const size_t bot_nrow, + const double bot_xori, + const double bot_yori, + const double bot_xinc, + const double bot_yinc, + const double bot_rotation, + const py::array_t &bot_values) +{ + pybind11::array_t result({ ncol, nrow, nlay }); + auto result_ = result.mutable_unchecked<3>(); + + // Access the np array without bounds checking for optimizing speed + auto xmid_ = xmid.unchecked<3>(); + auto ymid_ = ymid.unchecked<3>(); + auto zmid_ = zmid.unchecked<3>(); + + for (size_t i = 0; i < ncol; i++) { + for (size_t j = 0; j < nrow; j++) { + for (size_t k = 0; k < nlay; k++) { + // for every cell, project the center to the top and bottom surfaces + double xm = xmid_(i, j, k); + double ym = ymid_(i, j, k); + double zm = zmid_(i, j, k); + + // check if zm is NaN which can occur for inactive cells + if (std ::isnan(zm)) { + result_(i, j, k) = 0; + continue; + } + + double top_z = + regsurf::get_z_from_xy(xm, ym, top_xori, top_yori, top_xinc, top_yinc, + top_ncol, top_nrow, top_rotation, top_values); + double bot_z = + regsurf::get_z_from_xy(xm, ym, bot_xori, bot_yori, bot_xinc, bot_yinc, + bot_ncol, bot_nrow, bot_rotation, bot_values); + + // check if top_z or bot_z is NaN + if (std::isnan(top_z) || std::isnan(bot_z)) { + result_(i, j, k) = 0; + continue; + } + + if (zm >= top_z && zm <= bot_z) { // inside + result_(i, j, k) = 1; + } else { + result_(i, j, k) = 0; + } + } + } + } + return result; +} // grid_assign_value_between_surfaces + +} // namespace xtgeo::grid3d diff --git a/src/lib/src/regsurf/utilities.cpp b/src/lib/src/regsurf/utilities.cpp index d729d0691..a8e0804a0 100644 --- a/src/lib/src/regsurf/utilities.cpp +++ b/src/lib/src/regsurf/utilities.cpp @@ -230,4 +230,76 @@ find_cell_range(const double xmin, return { i_min, i_max, j_min, j_max }; } +/* + * Function to get the Z value at a given X, Y position on a regular 2D grid + * + * @param x The x-coordinate of the point + * @param y The y-coordinate of the point + * @param xori The x-coordinate of the regsurf origin + * @param yori The y-coordinate of the regsurf origin + * @param xinc The x-increment + * @param yinc The y-increment + * @param ncol The number of columns + * @param nrow The number of rows + * @param angle_deg The angle of rotation in degrees, anticlockwise from X + * @param values The 2D array of values on the regsurf (where Nan values are masked) + * @return The interpolated Z value at the given X, Y position + */ + +double +get_z_from_xy(const double x, + const double y, + const double xori, + const double yori, + const double xinc, + const double yinc, + const size_t ncol, + const size_t nrow, + const double angle_deg, + const py::array_t &values) +{ + // Convert the angle to radians + double angle_rad = angle_deg * M_PI / 180.0; + + Point p = { x, y, 0 }; + Point p_rel = inverse_rotate_and_translate(p, xori, yori, std::cos(angle_rad), + std::sin(angle_rad)); + + // Find the indices of the grid cell containing the point + int i_temp = static_cast(p_rel.x / xinc); + int j_temp = static_cast(p_rel.y / yinc); + + // Check if the point is outside the grid, and return NaN if it is + if (i_temp < 0 || i_temp >= static_cast(ncol - 1) || j_temp < 0 || + j_temp >= static_cast(nrow - 1)) { + return std::numeric_limits::quiet_NaN(); + } + + // Convert to size_t after validation + size_t i = static_cast(i_temp); + size_t j = static_cast(j_temp); + // Get the values at the corners of the cell + auto values_unchecked = + values.unchecked<2>(); // Access the array without bounds checking (faster) + double z11 = values_unchecked(i, j); + double z12 = values_unchecked(i, j + 1); + double z21 = values_unchecked(i + 1, j); + double z22 = values_unchecked(i + 1, j + 1); + + // Check if any of the corner values are NaN + if (std::isnan(z11) || std::isnan(z12) || std::isnan(z21) || std::isnan(z22)) { + return std::numeric_limits::quiet_NaN(); + } + + // Perform bilinear interpolation + double x1 = i * xinc; + double x2 = (i + 1) * xinc; + double y1 = j * yinc; + double y2 = (j + 1) * yinc; + + return geometry::interpolate_z_4p_regular(p_rel.x, p_rel.y, { x1, y1, z11 }, + { x2, y1, z21 }, { x1, y2, z12 }, + { x2, y2, z22 }); +} // get_z_from_xy + } // namespace xtgeo::regsurf diff --git a/src/xtgeo/grid3d/_grid_etc1.py b/src/xtgeo/grid3d/_grid_etc1.py index 987448af7..d3680467e 100644 --- a/src/xtgeo/grid3d/_grid_etc1.py +++ b/src/xtgeo/grid3d/_grid_etc1.py @@ -11,7 +11,7 @@ import pandas as pd from packaging.version import parse as versionparse -import xtgeo._internal as _internal +import xtgeo._internal as _internal # type: ignore from xtgeo import _cxtgeo from xtgeo.common import null_logger from xtgeo.common.calc import find_flip @@ -27,6 +27,7 @@ from xtgeo.common.types import Dimensions from xtgeo.grid3d import Grid from xtgeo.grid3d.types import METRIC + from xtgeo.surface.regular_surface import RegularSurface from xtgeo.xyz.points import Points logger = null_logger(__name__) @@ -290,6 +291,88 @@ def get_heights_above_ffl( return htop, hbot, hmid +def get_property_between_surfaces( + grid: Grid, + top: RegularSurface, + base: RegularSurface, + value: int = 1, + name: str = "between_surfaces", +) -> GridProperty: + """For a grid, create a grid property with value between two surfaces. + + The value would be zero elsewhere, or if surfaces has inactive nodes. + """ + if not isinstance(value, int) or value < 1: + raise ValueError(f"Value (integer) must be positive, >= 1, got: {value}") + + grid._xtgformat2() + logger.debug("Creating property between surfaces...") + + xmid, ymid, zmid = _internal.grid3d.grid_cell_centers( + grid.ncol, + grid.nrow, + grid.nlay, + grid._coordsv, + grid._zcornsv, + grid._actnumsv, + True, # asmasked, will give Nan for inactive cells; intentional here + ) + + top_ = top + base_ = base + if top.yflip == -1: + top_ = top.copy() + top_.make_lefthanded() + logger.debug("Top surface is right-handed, flipping a copy prior to operation") + if base.yflip == -1: + base_ = base.copy() + base_.make_lefthanded() + logger.debug("Base surface is right-handed, flipping a copy prior to operation") + + diff = base_ - top_ + if (diff.values).all() <= 0: + raise ValueError( + "Top surface must be equal or above base surface for all nodes" + ) + + # array is always 0, 1 integer + array = _internal.grid3d.grid_assign_value_between_surfaces( + grid.ncol, + grid.nrow, + grid.nlay, + xmid, + ymid, + zmid, + top_.ncol, + top_.nrow, + top_.xori, + top_.yori, + top_.xinc, + top_.yinc, + top_.rotation, + top_.values.filled(np.nan), + base_.ncol, + base_.nrow, + base_.xori, + base_.yori, + base_.xinc, + base_.yinc, + base_.rotation, + base_.values.filled(np.nan), + ) + + logger.debug("Creating property between surfaces... done") + + return GridProperty( + ncol=grid.ncol, + nrow=grid.nrow, + nlay=grid.nlay, + name=name, + values=array * value, + discrete=True, + ) + + def get_ijk( self: Grid, names: tuple[str, str, str] = ("IX", "JY", "KZ"), diff --git a/src/xtgeo/grid3d/grid.py b/src/xtgeo/grid3d/grid.py index 19ce693e5..303c51760 100644 --- a/src/xtgeo/grid3d/grid.py +++ b/src/xtgeo/grid3d/grid.py @@ -1648,6 +1648,31 @@ def get_heights_above_ffl( """ return _grid_etc1.get_heights_above_ffl(self, ffl=ffl, option=option) + def get_property_between_surfaces( + self, + top: xtgeo.RegularSurface, + base: xtgeo.RegularSurface, + value: int = 1, + name: str = "between_surfaces", + ) -> GridProperty: + """Returns a 3D GridProperty object with between two surfaces." + + Args: + top: The bounding top surface (RegularSurface object) + base: The bounding base surface (RegularSurface object) + value: An integer > 0 to assign to cells between surfaces, 1 as default + name: Name of the property, default is "between_surfaces" + + Returns: + xtgeo GridProperty object with if cell center is between surfaces, + otherwise 0. Note that the property wil be discrete if input value is an + integer, otherwise it will be continuous. + + .. versionadded:: 4.5 + + """ + return _grid_etc1.get_property_between_surfaces(self, top, base, value, name) + def get_ijk( self, names: tuple[str, str, str] = ("IX", "JY", "KZ"), diff --git a/tests/test_grid3d/test_grid.py b/tests/test_grid3d/test_grid.py index b4567bbfc..ad4beff11 100644 --- a/tests/test_grid3d/test_grid.py +++ b/tests/test_grid3d/test_grid.py @@ -544,6 +544,48 @@ def test_cell_height_above_ffl(testdata_path): assert hmid.values[4, 5, 0] == pytest.approx(22.4055) +def test_get_property_between_surfaces(testdata_path): + """Generate a marker property between two surfaces.""" + grd = xtgeo.grid_from_file(testdata_path / REEKFIL4) + + surf1 = xtgeo.surface_from_grid3d(grd) + surf1.fill() + surf1.values = 1650 + surf2 = surf1.copy() + surf2.values = 1700 + + prop = grd.get_property_between_surfaces(surf1, surf2) + + assert prop.values.sum() == 137269 # verified with similar method in RMS + + # multiply values with 2 + prop2 = grd.get_property_between_surfaces(surf1, surf2, value=2) + assert prop2.values.sum() == 137269 * 2 + + # swap one if the surfaces so yflip becomes -1 + surf1.make_righthanded() + prop2 = grd.get_property_between_surfaces(surf1, surf2) + assert prop2.values.sum() == 137269 + + +def test_get_property_between_surfaces_w_holes(testdata_path): + """Generate a marker property between two surfaces, where surfaces has holes.""" + grd = xtgeo.grid_from_file(testdata_path / REEKFIL4) + + surf1 = xtgeo.surface_from_grid3d(grd) + surf1.fill() + surf1.values = 1650 + surf2 = surf1.copy() + surf2.values = 1700 + + surf1.values.mask[60:70, 70:75] = True + surf2.values.mask[50:70, 60:71] = True + + prop = grd.get_property_between_surfaces(surf1, surf2) + + assert prop.values.sum() == 131130 # verified manually in RMS + + def test_bad_egrid_ends_before_kw(tmp_path): egrid_file = tmp_path / "test.egrid" with open(egrid_file, "wb") as fh: @@ -758,7 +800,7 @@ def test_get_vtk_geometries_box(show_plot): if show_plot: try: - import pyvista as pv + import pyvista as pv # type: ignore except ImportError: warnings.warn("show_plot is active but no pyvista installed") return @@ -781,7 +823,7 @@ def test_get_vtk_geometries_emerald(show_plot, testdata_path): if show_plot: try: - import pyvista as pv + import pyvista as pv # type: ignore except ImportError: warnings.warn("show_plot is active but no pyvista installed") return diff --git a/tests/test_internal/test_cell_grid3d.py b/tests/test_internal/test_cell_grid3d.py index 17450bab4..843b6fb66 100644 --- a/tests/test_internal/test_cell_grid3d.py +++ b/tests/test_internal/test_cell_grid3d.py @@ -14,11 +14,24 @@ import xtgeo import xtgeo._internal as _internal # type: ignore -from xtgeo.common.log import null_logger +from xtgeo.common.log import functimer, null_logger logger = null_logger(__name__) +@pytest.fixture(scope="module", name="get_drogondata") +def fixture_get_drogondata(testdata_path): + grid = xtgeo.grid_from_file(f"{testdata_path}/3dgrids/drogon/2/geogrid.roff") + poro = xtgeo.gridproperty_from_file( + f"{testdata_path}/3dgrids/drogon/2/geogrid--phit.roff" + ) + facies = xtgeo.gridproperty_from_file( + f"{testdata_path}/3dgrids/drogon/2/geogrid--facies.roff" + ) + + return grid, poro, facies + + @pytest.mark.parametrize( "cell, expected_firstcorner, expected_lastcorner", [ @@ -131,3 +144,46 @@ def test_get_depth_in_cell(testdata_path, x, y, cell, position, expected): assert np.isnan(depth) else: assert depth == pytest.approx(expected) + + +@functimer(output="info") +def test_grid_cell_centers(get_drogondata): + """Test cell centers from a grid; assertions are manually verified in RMS.""" + + grid, _, _ = get_drogondata # total cells 899944 + + xcor, ycor, zcor = _internal.grid3d.grid_cell_centers( + grid.ncol, + grid.nrow, + grid.nlay, + grid._coordsv, + grid._zcornsv, + grid._actnumsv, + True, + ) + + assert isinstance(xcor, np.ndarray) + assert isinstance(ycor, np.ndarray) + assert isinstance(zcor, np.ndarray) + + assert xcor.shape == (grid.ncol, grid.nrow, grid.nlay) + + assert np.nanmean(xcor) == pytest.approx(461461.6, abs=0.1) + assert np.nanmean(ycor) == pytest.approx(5933128.2, abs=0.1) + assert np.nanmean(zcor) == pytest.approx(1731.50, abs=0.1) + + assert np.nanstd(xcor) == pytest.approx(2260.1489, abs=0.1) + assert np.nanstd(ycor) == pytest.approx(2945.992, abs=0.1) + assert np.nanstd(zcor) == pytest.approx(67.27, abs=0.1) + + assert np.nansum(xcor) == pytest.approx(301686967304.8, abs=0.1) # RMS 301686967112 + assert np.nansum(ycor) == pytest.approx(3878865619174.5, abs=0.1) # 3878865622519 + assert np.nansum(zcor) == pytest.approx(1131994843.5, abs=0.1) # 1131994843 + + assert xcor[75, 23, 37] == pytest.approx(461837.19, abs=0.1) + assert ycor[75, 23, 37] == pytest.approx(5937352, abs=0.1) + assert zcor[75, 23, 37] == pytest.approx(1930.57, abs=0.1) + + assert np.isnan(xcor[62, 33, 37]) + assert np.isnan(ycor[62, 33, 37]) + assert np.isnan(zcor[62, 33, 37]) diff --git a/tests/test_internal/test_regsurf.py b/tests/test_internal/test_regsurf.py index 140f54814..eb4de505b 100644 --- a/tests/test_internal/test_regsurf.py +++ b/tests/test_internal/test_regsurf.py @@ -14,7 +14,7 @@ import xtgeo import xtgeo._internal as _internal # type: ignore -from xtgeo.common.log import null_logger +from xtgeo.common.log import functimer, null_logger logger = null_logger(__name__) @@ -301,3 +301,96 @@ def sample_grid(): # check if the top layer is exactly the same for all threads assert np.array_equal(top, keep_top_store["keep_top"]) + + +@pytest.mark.parametrize( + "x, y, rotation, expected", + [ + (2.746, 1.262, 0, 17.73800), + (2.096, 4.897, 0, 17.473), + (2.746, 2.262, 0, np.nan), + (1.9999999999, 1.999999999, 0, 14.0), + (1999, 1999, 0, np.nan), # far outside + (2.746, 1.262, 10, 18.3065), + (2.096, 4.897, 10, 21.9457), + (2.746, 1.262, 20, 18.319), + (2.096, 4.897, 20, 25.751), + (-0.470, -3.354, 210, np.nan), + (-0.013, -5.148, 210, 19.963), # rms gets 19.959 + (-1.433, -5.359, 210, 27.448), # rms gets 27.452 + ], +) +def test_get_z_from_xy_simple(x, y, rotation, expected): + """Test values are manually verified by inspection in RMS""" + values = np.array(range(30)).reshape(5, 6).astype(float) + + values[2, 3] = np.nan + values[4, 5] = np.nan + + surf = xtgeo.RegularSurface( + xori=0, yori=0, xinc=1, yinc=1, ncol=5, nrow=6, rotation=rotation, values=values + ) + + z_value = _internal.regsurf.get_z_from_xy( + x, + y, + surf.xori, + surf.yori, + surf.xinc, + surf.yinc, + surf.ncol, + surf.nrow, + surf.rotation, + surf.values, + ) + + if np.isnan(expected): + assert np.isnan(z_value) + else: + assert z_value == pytest.approx(expected, abs=0.001) + + +@functimer +def test_get_z_from_xy(get_drogondata): + _, _, _, srf = get_drogondata + + surf = srf.copy() + + z_value = _internal.regsurf.get_z_from_xy( + 460103.00, + 5934855.00, + surf.xori, + surf.yori, + surf.xinc, + surf.yinc, + surf.ncol, + surf.nrow, + surf.rotation, + surf.values, + ) + + assert z_value == pytest.approx(1594.303, abs=0.001) + + # make ntimes random points and check the time (cf. @functimer with debug logging) + ntimes = 100000 + xmin = 458000 + xmax = 468000 + ymin = 5927000 + ymax = 5938000 + logger.debug("Start random points loop... %s", ntimes) + for _ in range(ntimes): + x = np.random.uniform(xmin, xmax) + y = np.random.uniform(ymin, ymax) + _internal.regsurf.get_z_from_xy( + x, + y, + surf.xori, + surf.yori, + surf.xinc, + surf.yinc, + surf.ncol, + surf.nrow, + surf.rotation, + surf.values, + ) + logger.debug("End random points loop")