Merge pull request #278 from jcrivenaes/grid-axes-swap

Grid axes swap

src/xtgeo/clib/xtg/grd3d_reverse_jrows.c

@@ -0,0 +1,122 @@
+/*
+ ***************************************************************************************
+ *
+ * Do reversing / flipping  of J coordinates in the grid
+ *
+ ***************************************************************************************
+ */
+
+#include "logger.h"
+#include "libxtg.h"
+#include "libxtg_.h"
+
+/*
+ ***************************************************************************************
+ *
+ * NAME:
+ *    grd3d_reverse_jrows.c
+ *
+ * DESCRIPTION:
+ *    Reverting J axes, i.e. grid goes from right-handed to left-handed or vice-versa
+ *
+ * ARGUMENTS:
+ *    nx, ny, nz       i     NCOL, NROW, NLAY dimens
+ *    p_coord_v       i/o    Cornerlines, pillars
+ *    p_zcorn_v       i/o    ZCORN
+ *    p_actnum_v      i/o    ACTNUM values
+ *
+ * RETURNS:
+ *    Function: 0: upon success. Update pointers in-place
+ *
+ * TODO/ISSUES/BUGS:
+ *
+ * LICENCE:
+ *    cf. XTGeo LICENSE
+ ***************************************************************************************
+ */
+
+int grd3d_reverse_jrows (
+    int nx,
+    int ny,
+    int nz,
+    double *p_coord_v,
+    double *p_zcorn_v,
+    int *p_actnum_v
+    )
+
+{
+
+    long ic;
+    int i, j, k;
+    long ntot, ib1, ib2;
+    double *tmpcv;
+    int *itmp;
+
+    logger_init(__FILE__, __FUNCTION__);
+    logger_info(__LINE__, "Flip/swap J axis");
+
+    ntot = (nx + 1) * (ny + 1) * 6;  /*  COORD values */
+    tmpcv = calloc(ntot, sizeof(double));
+
+    logger_info(__LINE__, "J swapping COORD");
+    logger_info(__LINE__, "NX NY NZ %d %d %d", nx, ny, nz);
+    ib2 = 0;
+    for (j = 1; j <= ny + 1; j++) {
+        for (i = 1; i <= nx + 1; i++) {
+
+            int jx = ny + 1 - j + 1;
+
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+0] = p_coord_v[6*((j-1)*(nx+1)+i-1)+0];
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+1] = p_coord_v[6*((j-1)*(nx+1)+i-1)+1];
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+2] = p_coord_v[6*((j-1)*(nx+1)+i-1)+2];
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+3] = p_coord_v[6*((j-1)*(nx+1)+i-1)+3];
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+4] = p_coord_v[6*((j-1)*(nx+1)+i-1)+4];
+            tmpcv[6*((jx-1)*(nx+1)+i-1)+5] = p_coord_v[6*((j-1)*(nx+1)+i-1)+5];
+        }
+    }
+    for (ic = 0; ic < ntot; ic++) p_coord_v[ic] = tmpcv[ic];
+    free(tmpcv);
+
+    /*  ***************************************************************************** */
+    logger_info(__LINE__, "J swapping ZCORN");
+
+    ntot = nx * ny * (nz + 1) * 4;  /*  ZCORN values */
+    tmpcv = calloc(ntot, sizeof(double));
+
+    ib2 = 0;
+    for (k = 1; k <= nz + 1; k++) {
+        for (j = ny; j >= 1; j--) {
+            for (i = 1; i <= nx; i++) {
+                ib1 = x_ijk2ib(i, j, k, nx, ny, nz + 1, 0);
+                tmpcv[4*ib2 + 0] = p_zcorn_v[4*ib1 + 2];
+                tmpcv[4*ib2 + 1] = p_zcorn_v[4*ib1 + 3];
+                tmpcv[4*ib2 + 2] = p_zcorn_v[4*ib1 + 0];
+                tmpcv[4*ib2 + 3] = p_zcorn_v[4*ib1 + 1];
+                ib2++;
+            }
+        }
+    }
+    for (ic = 0; ic < ntot; ic++) p_zcorn_v[ic] = tmpcv[ic];
+    free(tmpcv);
+
+    /*  ***************************************************************************** */
+    logger_info(__LINE__, "J swapping ACTNUM");
+
+    ntot = nx * ny * nz;  /*  ACTNUM values */
+    itmp = calloc(ntot, sizeof(int));
+
+
+    for (k = 1; k <= nz; k++) {
+        for (j = 1; j <= ny; j++) {
+            for (i = 1; i <= nx; i++) {
+                ib1 = x_ijk2ib(i, j, k, nx, ny, nz, 0);
+                ib2 = x_ijk2ib(i, ny - j + 1, k, nx, ny, nz, 0);
+                itmp[ib1] = p_actnum_v[ib2];
+            }
+        }
+    }
+    for (ic = 0; ic < ntot; ic++) p_actnum_v[ic] = itmp[ic];
+    free(itmp);
+
+    return EXIT_SUCCESS;
+}

src/xtgeo/clib/xtg/grd3d_translate.c

@@ -1,5 +1,5 @@
 /*
- *******************************************************************************
+****************************************************************************************
  *
  * NAME:
  *    grd3d_translate.c
@@ -29,10 +29,10 @@
  *    cf. XTGeo LICENSE
  *
  * CODING STANDARD:
- *    Pure C, 79 columns, Stroustrup do/if braces
- *******************************************************************************
+ *    Pure C, 88 columns, Stroustrup do/if braces
+ ***************************************************************************************
  */
-
+#include "logger.h"
 #include "libxtg.h"
 
 int grd3d_translate (
@@ -46,18 +46,15 @@ int grd3d_translate (
 		     double yshift,
 		     double zshift,
 		     double *p_coord_v,
-		     double *p_zcorn_v,
-		     int   debug
+		     double *p_zcorn_v
 		     )
 
 {
     /* locals */
     int    i, j, ic, ib, nzcorn, iok=0;
-    char sub[24]="grd3d_translate";
-
-    xtgverbose(debug);
 
-    xtg_speak(sub,2,"Entering routine ...%s", sub);
+    logger_init(__FILE__, __FUNCTION__);
+    logger_info(__LINE__, "Do translation or pure flipping");
 
     if ((xflip == 1 || xflip == -1) &&
 	(yflip == 1 || yflip == -1) &&
@@ -67,7 +64,7 @@ int grd3d_translate (
 
     if (iok == 0) {
 	/* flip out of range */
-	xtg_warn(sub,1,"Error in flips ...%d %d %d", xflip, yflip, zflip);
+	logger_warn(__LINE__,"Error in flips ...%d %d %d", xflip, yflip, zflip);
 	return (-1);
     }
 
@@ -91,6 +88,6 @@ int grd3d_translate (
 	p_zcorn_v[ic]=zflip*(p_zcorn_v[ic] + zshift);
     }
 
-    xtg_speak(sub,2,"Exit from routine %s", sub);
+    logger_info(__LINE__,"Exit from routine");
     return(0);
 }

src/xtgeo/clib/xtg/libxtg.h

@@ -1770,10 +1770,17 @@ int grd3d_translate (
     double yshift,
     double zshift,
     double *p_coord_v,
-    double *p_zcorn_v,
-    int   debug
+    double *p_zcorn_v
     );
 
+int grd3d_reverse_jrows (
+    int   nx,
+    int   ny,
+    int   nz,
+    double *p_coord_v,
+    double *p_zcorn_v,
+    int *p_actnum_v
+    );
 
 int grd3d_point_val_crange(
     double x,

src/xtgeo/grid3d/_grid_etc1.py

@@ -903,14 +903,41 @@ def translate_coordinates(self, translate=(0, 0, 0), flip=(1, 1, 1)):
         tz,
         self._p_coord_v,
         self._p_zcorn_v,
-        XTGDEBUG,
     )
     if ier != 0:
         raise RuntimeError("Something went wrong in translate, code: {}".format(ier))
 
     logger.info("Translation of coords done")
 
 
+def reverse_row_axis(self, ijk_handedness=None):
+    """Reverse rows (aka flip) for geometry and assosiated properties"""
+
+    if ijk_handedness == self.ijk_handedness:
+        return
+
+    ier = _cxtgeo.grd3d_reverse_jrows(
+        self._ncol,
+        self._nrow,
+        self._nlay,
+        self._p_coord_v,
+        self._p_zcorn_v,
+        self._p_actnum_v,
+    )
+
+    if ier != 0:
+        raise RuntimeError("Something went wrong in jswapping, code: {}".format(ier))
+
+    if self._props is None:
+        return
+
+    # do it for properties
+    for prp in self._props.props:
+        prp.values = prp.values[:, ::-1, :]
+
+    logger.info("Reversing of rows done")
+
+
 def report_zone_mismatch(  # pylint: disable=too-many-statements
     self,
     well=None,

src/xtgeo/grid3d/grid.py

@@ -171,6 +171,7 @@ def __init__(self, *args, **kwargs):
 
         self._props = None  # None or a GridProperties instance
         self._subgrids = None  # A python dict if subgrids are given
+        self._ijk_handedness = None
 
         # Simulators like Eclipse may have a dual poro/perm model
         self._dualporo = False
@@ -256,6 +257,32 @@ def dimensions(self):
         """3-tuple: The grid dimensions as a tuple of 3 integers (read only)"""
         return (self._ncol, self._nrow, self._nlay)
 
+    @property
+    def ijk_handedness(self):
+        """IJK handedness for grids, "right" or "left"
+
+        For a non-rotated grid with K increasing with depth, 'left' is corner in
+        lower-left, while 'right' is origin in upper-left corner.
+
+        """
+        nflip = _grid_etc1.estimate_flip(self)
+        if nflip == 1:
+            self._ijk_handedness = "left"
+        elif nflip == -1:
+            self._ijk_handedness = "right"
+        else:
+            self._ijk_handedness = None  # cannot determine
+
+        return self._ijk_handedness
+
+    @ijk_handedness.setter
+    def ijk_handedness(self, value):
+        if value in ("right", "left"):
+            self.reverse_row_axis(ijk_handedness=value)
+        else:
+            raise ValueError("The value must be 'right' or 'left'")
+        self._ijk_handedness = value
+
     @property
     def subgrids(self):
         """:obj:`list` of :obj:`int`: A dictionary with subgrid name and
@@ -461,7 +488,7 @@ def create_box(
             increment (tuple of float): Grid increments (xinc, yinc, zinc)
             rotation (float): Roations in degrees, anticlock from X axis.
             flip (int): If 1, grid origin is lower left and left-handed; if 1, origin
-            is upper left and right-handed.
+            is upper left and right-handed (row flip).
 
         Returns:
             Instance is updated (previous instance content will be erased)
@@ -823,7 +850,11 @@ def estimate_design(self, nsub=None):
         return res
 
     def estimate_flip(self):
-        """Determine flip (handedness) of grid"""
+        """Determine flip (handedness) of grid returns; 1 for left-handed,
+        -1 for right-handed.
+
+        .. seealso:: :py:attr:`~ijk_handedness`
+        """
 
         return _grid_etc1.estimate_flip(self)
 
@@ -980,12 +1011,13 @@ def get_dz(self, name="dZ", flip=True, asmasked=True, mask=None):
 
         Args:
             name (str): name of property
-            flip (bool): Use False for Petrel grids (experimental)
+            flip (bool): Use False for Petrel grids were Z is negative down
+                (experimental)
             asmasked (bool): True if only for active cells, False for all cells
             mask (bool): Deprecated, use asmasked instead!
 
         Returns:
-            A XTGeo GridProperty object
+            A XTGeo GridProperty object for flip
         """
         if mask is not None:
             asmasked = self._evaluate_mask(mask)
@@ -1368,8 +1400,11 @@ def reduce_to_one_layer(self):
     def translate_coordinates(self, translate=(0, 0, 0), flip=(1, 1, 1)):
         """Translate (move) and/or flip grid coordinates in 3D.
 
+        By 'flip' here, it means that the full coordinate array are multiplied
+        with -1.
+
         Args:
-            translate (tuple): Tranlattion distance in X, Y, Z coordinates
+            translate (tuple): Translation distance in X, Y, Z coordinates
             flip (tuple): Flip array. The flip values must be 1 or -1.
 
         Raises:
@@ -1378,6 +1413,37 @@ def translate_coordinates(self, translate=(0, 0, 0), flip=(1, 1, 1)):
 
         _grid_etc1.translate_coordinates(self, translate=translate, flip=flip)
 
+    def reverse_row_axis(self, ijk_handedness=None):
+        """Reverse the row axis (J indices).
+
+        This means that IJK system will switched between a left vs right handed system.
+        It is here (by using ijk_handedness key), possible to set a wanted stated.
+
+        Note that properties that are assosiated with the grid (through the
+        :py:attr:`~gridprops` or :py:attr:`~props` attribute) will also be
+        reversed (which is desirable).
+
+        Args:
+            ijk_handedness (str): If set to "right" or "left", do only reverse rows if
+                handedness is not already achieved.
+
+        Example::
+
+            grd = xtgeo.Grid("somefile.roff")
+            prop1 = xtgeo.GridProperty("somepropfile1.roff")
+            prop2 = xtgeo.GridProperty("somepropfile2.roff")
+
+            grd.props = [prop1, prop2]
+
+            # secure that the grid geometry is IJK right-handed
+            grd.reverse_row_axis(ijk_handedness="right")
+
+        .. versionadded:: 2.4.0
+
+        """
+
+        _grid_etc1.reverse_row_axis(self, ijk_handedness=ijk_handedness)
+
     def make_zconsistent(self, zsep=1e-5):
         """Make the 3D grid consistent in Z, by a minimal gap (zsep).