Merge pull request #142 from jcrivenaes/xsect-seis

simplified cube xsections

Author: jcrivenaes

src/xtgeo/cube/_cube_utils.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 
+import xtgeo
 import xtgeo.cxtgeo.cxtgeo as _cxtgeo
 from xtgeo.common import XTGeoDialog
 
@@ -261,12 +262,22 @@ def get_xy_value_from_ij(self, iloc, jloc, ixline=False, zerobased=False):
 
 
 def get_randomline(
-    self, fencespec, zmin=None, zmax=None, zincrement=None, sampling="nearest"
+    self,
+    fencespec,
+    zmin=None,
+    zmax=None,
+    zincrement=None,
+    hincrement=None,
+    atleast=5,
+    extend=2,
+    sampling="nearest",
 ):
     """Get a random line from a fence spesification"""
 
-    if not isinstance(fencespec, np.ndarray):
-        raise ValueError("Fence is not a numpy array")
+    if isinstance(fencespec, xtgeo.Polygons):
+        logger.info("Estimate hincrement from Polygons instance...")
+        fencespec = _get_randomline_fence(self, fencespec, hincrement, atleast, extend)
+        logger.info("Estimate hincrement from Polygons instance... DONE")
 
     if not len(fencespec.shape) == 2:
         raise ValueError("Fence is not a 2D numpy")
@@ -325,6 +336,20 @@ def get_randomline(
     return (hcoords[0], hcoords[-1], zmin, zmax, arr)
 
 
+def _get_randomline_fence(self, fencespec, hincrement, atleast, extend):
+    """Compute a resampled fence from a Polygons instance"""
+
+    if hincrement is None:
+        avgdxdy = 0.5 * (self.xinc + self.yinc)
+        distance = 0.5 * avgdxdy
+
+    logger.info("Getting fence from a Polygons instance...")
+    return fencespec.get_fence(
+        distance=distance, atleast=atleast, extend=extend, asnumpy=True
+    )
+    logger.info("Getting fence from a Polygons instance... DONE")
+
+
 # copy (update) values from SWIG carray to numpy, 3D array, Fortran order
 # to be DEPRECATED
 def update_values(cube):

src/xtgeo/cube/cube1.py

@@ -8,17 +8,15 @@
 
 import numpy as np
 
+import xtgeo
 from xtgeo.common import XTGeoDialog
 from xtgeo.common import XTGDescription
-from xtgeo.common import constants
 
 from xtgeo.cube import _cube_import
 from xtgeo.cube import _cube_export
 from xtgeo.cube import _cube_utils
 from xtgeo.cube import _cube_roxapi
 
-UNDEF = constants.UNDEF
-UNDEF_LIMIT = constants.UNDEF_LIMIT
 
 xtg = XTGeoDialog()
 logger = xtg.functionlogger(__name__)
@@ -95,9 +93,6 @@ class Cube(object):  # pylint: disable=too-many-public-methods
     def __init__(self, *args, **kwargs):
         """Initiate a Cube instance."""
 
-        clsname = "{}.{}".format(type(self).__module__, type(self).__name__)
-        logger.info(clsname)
-
         self._filesrc = None
         self._segyfile = None
         self._ilines = None
@@ -117,8 +112,6 @@ def __init__(self, *args, **kwargs):
         self._xlines = np.array(range(1, 3 + 1), dtype=np.int32)
         self._rotation = 0.0
         self._traceidcodes = np.ones((5, 3), dtype=np.int32)
-        self._undef = UNDEF
-        self._undef_limit = UNDEF_LIMIT
 
         if len(args) >= 1:
             fformat = kwargs.get("fformat", "guess")
@@ -565,7 +558,15 @@ def get_xy_value_from_ij(self, iloc, jloc, ixline=False, zerobased=False):
     # =========================================================================
 
     def get_randomline(
-        self, fencespec, zmin=None, zmax=None, zincrement=None, sampling="nearest"
+        self,
+        fencespec,
+        zmin=None,
+        zmax=None,
+        zincrement=None,
+        hincrement=None,
+        atleast=5,
+        extend=2,
+        sampling="nearest",
     ):
         """Get a randomline from a fence spesification.
 
@@ -576,13 +577,24 @@ def get_randomline(
 
         The input fencespec is a 2D numpy where each row is X, Y, Z, HLEN,
         where X, Y are UTM coordinates, Z is depth/time, and HLEN is a
-        length along the fence.
+        length along the fence, or a Polygons instance with HLEN present.
+
+        It is important that the HLEN array has a constant increment and ideally
+        a sampling that is less than the Cube resolution.
 
         Args:
-            fencespec (np): 2D numpy with X, Y, Z, HLEN as rows.
+            fencespec (~np.ndarray or :class:`~xtgeo.xyz.polygons.Polygons`):
+                2D numpy with X, Y, Z, HLEN as rows or a xtgeo Polygons() object.
             zmin (float): Minimum Z (default is Cube Z minima/origin)
             zmax (float): Maximum Z (default is Cube Z maximum)
-            zincrement (float): Sampling, default is Cube ZINC/2
+            zincrement (float): Sampling vertically, default is Cube ZINC/2
+            hincrement (float or bool): Resampling horizontally. This applies only
+                if the fencespec is a Polygons() instance. If None (default),
+                the distance will be deduced automatically.
+            atleast (int): Minimum number of horizontal samples (only if
+                fencespec is a Polygons instance)
+            extend (int): Extend with extend*hincrement in both ends (only if
+                fencespec is a Polygons instance)
             sampling (str): Algorithm, 'nearest' or 'trilinear' (first is
                 faster, second is more precise for continuous fields)
 
@@ -592,16 +604,25 @@ def get_randomline(
         Raises:
             ValueError: Input fence is not according to spec.
 
-        """
+        .. versionadded:: 2.1.0 support for Polygons() as fencespec, and keywords
+           hincrement, atleast and sampling
 
+        """
+        if not isinstance(fencespec, (np.ndarray, xtgeo.Polygons)):
+            raise ValueError("fencespec must be a numpy or a Polygons() object")
+        logger.info("Getting randomline...")
         res = _cube_utils.get_randomline(
             self,
             fencespec,
             zmin=zmin,
             zmax=zmax,
             zincrement=zincrement,
+            hincrement=hincrement,
+            atleast=atleast,
+            extend=extend,
             sampling=sampling,
         )
+        logger.info("Getting randomline... DONE")
         return res
 
     # =========================================================================

src/xtgeo/cxtgeo/clib/src/pol_refine.c

@@ -69,7 +69,7 @@ int pol_refine (
 
 	if (n>1) {
 	    len=usedist/n;
-
+            xtg_speak(s, 2, "LEN is %d", len);
 	    for (ii=1;ii<n; ii++) {
 
 		frac = (len/usedist)*ii;

src/xtgeo/cxtgeo/clib/src/pol_resampling.c

@@ -112,7 +112,7 @@ int pol_resampling(
     thlenx = calloc(naddr, sizeof(double));
 
     /* ========================================================================
-     * Part one, look at current input, re-estimate smapling distance, and
+     * Part one, look at current input, re-estimate sampling distance, and
      * find extend angles using certain criteria
      * ========================================================================
      */

src/xtgeo/xyz/_xyz_oper.py

@@ -137,7 +137,7 @@ def get_fence(self, distance=20, atleast=5, extend=2, name=None, asnumpy=True):
     The atleast parameter will win over the distance, meaning that if total length
     horizontally is 50, and distance is set to 20, the actual length will be 50/5=10
 
-"""
+    """
 
     if len(self._df) < 2:
         xtg.warn("Well does not enough points in interval, outside range?")
@@ -150,11 +150,13 @@ def get_fence(self, distance=20, atleast=5, extend=2, name=None, asnumpy=True):
 
     nbuf = 1000000
 
+    logger.info("%s %s %s %s %s ", distance, atleast, extend, name, asnumpy)
     npxarr = np.zeros(nbuf, dtype=np.float64)
     npyarr = np.zeros(nbuf, dtype=np.float64)
     npzarr = np.zeros(nbuf, dtype=np.float64)
     npharr = np.zeros(nbuf, dtype=np.float64)
 
+    logger.info("Calling C routine...")
     # C function:
     ier, npxarr, npyarr, npzarr, npharr, nlen = _cxtgeo.pol_resampling(
         self._df[self.xname].values,
@@ -169,6 +171,7 @@ def get_fence(self, distance=20, atleast=5, extend=2, name=None, asnumpy=True):
         0,
         XTGDEBUG
     )
+    logger.info("Calling C routine... DONE")
 
     if ier != 0:
         raise RuntimeError("Nonzero code from_cxtgeo.pol_resampling")

src/xtgeo/xyz/polygons.py

@@ -463,18 +463,19 @@ def get_fence(self, distance=20, atleast=5, extend=2, name=None, asnumpy=True):
             extend (int): Number of samples to extend at each end
             name (str): Name of polygon (if asnumpy=False)
             asnumpy (bool): Return a [:, 5] numpy array with
-                columns X, Y, Z, H, dH
+                columns X.., Y.., Z.., HLEN, dH
 
         Returns:
             A numpy array (if asnumpy=True) or a new Polygons() object
 
         .. versionadded:: 2.1.0
         """
-
+        logger.info("Getting fence within a Polygons instance...")
         return _xyz_oper.get_fence(
             self, distance=distance, atleast=atleast, extend=extend,
             name=name, asnumpy=asnumpy
         )
+        logger.info("Getting fence within a Polygons instance... DONE")
 
     # =========================================================================
     # Operations restricted to inside/outside polygons

tests/test_cube/test_cube.py

@@ -2,8 +2,11 @@
 import os
 from os.path import join
 
+import numpy as np
+import pandas as pd
 import pytest
 
+import xtgeo
 from xtgeo.cube import Cube
 from xtgeo.common import XTGeoDialog
 
@@ -311,3 +314,32 @@ def test_cube_swapaxes():
     tsetup.assert_almostequal(diff.mean(), 0.0, 0.000001)
     tsetup.assert_almostequal(diff.std(), 0.0, 0.000001)
     assert incube.ilines.size == incube.ncol
+
+
+def test_cube_randomline():
+    """Import a cube, and compute a randomline given a simple Polygon"""
+
+    # import matplotlib.pyplot as plt
+
+    incube = Cube(SFILE4)
+
+    # make a polyline with two points
+    dfr = pd.DataFrame(
+        np.array([[778133, 6737650, 2000, 1], [776880, 6738820, 2000, 1]]),
+        columns=["X_UTME", "Y_UTMN", "Z_TVDSS", "POLY_ID"],
+    )
+    poly = xtgeo.Polygons()
+    poly.dataframe = dfr
+
+    logger.info("Generate random line...")
+    hmin, hmax, vmin, vmax, random = incube.get_randomline(poly)
+
+    tsetup.assert_almostequal(hmin, -15.655932861802714, 0.001)
+    tsetup.assert_almostequal(random.mean(), -11.8755, 0.001)
+
+    # plt.figure()
+    # plt.imshow(random, cmap='seismic', interpolation='sinc',
+    #            extent=(hmin, hmax, vmax, vmin))
+    # plt.axis('tight')
+    # plt.colorbar()
+    # plt.show()