Merge pull request #140 from jcrivenaes/dev-xsection

Dev xsection

Author: jcrivenaes

docs/usextgeoroxar.rst

@@ -69,6 +69,72 @@ Take a surface in RMS and multiply values with 2:
     surf.to_roxar(project, 'TopReek', 'DS_tmp')
 
 
+Do operations on surfaces, also inside polygons:
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Find the diff maps in time domain, of the main surfaces. Also make a
+a version where cut by polygons where surfaces has interp (minimum
+common multiplum)
+
+.. code-block:: python
+
+   import xtgeo
+   from fmu.config import utilities as ut
+
+   CFG = ut.yaml_load("../../fmuconfig/output/global_variables.yml")["rms"]
+
+   # ========= SETTINGS ===================================================================
+
+   PRJ = project  # noqa
+
+   # input
+   TSCAT1 = "TS_interp_raw_ow"
+   PCAT = "TL_interp_raw_approx_outline"
+
+
+   # output
+   ISCAT1 = "IS_twt_main_interp_raw_ow"
+   ISCAT2 = "IS_twt_main_interp_raw_ow_cut"
+
+   # ========= END SETTINGS ===============================================================
+
+
+   def main():
+
+       topmainzones = CFG["horizons"]["TOP_MAINRES"]
+       mainzones = CFG["zones"]["MAIN_ZONES"]
+       for znum, mzone in enumerate(mainzones):
+
+           surf1 = xtgeo.surface_from_roxar(PRJ, topmainzones[znum], TSCAT1)
+           surf2 = xtgeo.surface_from_roxar(PRJ, topmainzones[znum + 1], TSCAT1)
+
+           diff = surf2.copy()
+           diff.values -= surf1.values
+           diff.to_roxar(PRJ, mzone, ISCAT1, stype="zones")
+           print("Store {} at {}".format(mzone, ISCAT1))
+
+           # extract differences inside a polygon and compute min/max values:
+
+           poly = xtgeo.polygons_from_roxar(PRJ, topmainzones[znum], PCAT)
+           surf1.eli_outside(poly)
+           surf2.eli_outside(poly)
+           diff2 = surf2.copy()
+           diff2.values -= surf1.values
+           print(
+              "Min and max values inside polygons {} : {} (negative OK) for {}".format(
+                    diff2.values.min(), diff2.values.max(), mzone
+                    )
+                )
+           diff2.to_roxar(PRJ, mzone, ISCAT2, stype="zones")
+           print("Store cut surface {} at {}".format(mzone, ISCAT2))
+
+
+    if __name__ == "__main__":
+        main()
+        print("Done, see <{}> and <{}>".format(ISCAT1, ISCAT2))
+
+
+
 3D grid data
 ------------
 

src/xtgeo/surface/_regsurf_export.py

@@ -5,6 +5,7 @@
 from __future__ import division, absolute_import
 from __future__ import print_function
 
+import xtgeo
 import xtgeo.cxtgeo.cxtgeo as _cxtgeo  # pylint: disable=import-error
 from xtgeo.common import XTGeoDialog
 
@@ -25,7 +26,7 @@ def export_irap_ascii(self, mfile):
     zmin = self.values.min()
     zmax = self.values.max()
 
-    vals = self.get_values1d(fill_value=self.undef)
+    vals = self.get_values1d(fill_value=xtgeo.UNDEF)
     logger.debug("SHAPE %s %s", vals.shape, vals.dtype)
 
     ier = _cxtgeo.surf_export_irap_ascii(
@@ -52,7 +53,7 @@ def export_irap_ascii(self, mfile):
 def export_irap_binary(self, mfile):
     """Export to Irap RMS binary format."""
 
-    vals = self.get_values1d(fill_value=self.undef)
+    vals = self.get_values1d(fill_value=xtgeo.UNDEF)
     ier = _cxtgeo.surf_export_irap_bin(
         mfile,
         self._ncol,
@@ -76,7 +77,7 @@ def export_irap_binary(self, mfile):
 def export_ijxyz_ascii(self, mfile):
     """Export to DSG IJXYZ ascii format."""
 
-    vals = self.get_values1d(fill_value=self.undef)
+    vals = self.get_values1d(fill_value=xtgeo.UNDEF)
     ier = _cxtgeo.surf_export_ijxyz(
         mfile,
         self._ncol,
@@ -117,7 +118,7 @@ def export_zmap_ascii(self, mfile):
 
     yinc = scopy._yinc * scopy._yflip
 
-    vals = scopy.get_values1d(order="F", asmasked=False, fill_value=self.undef)
+    vals = scopy.get_values1d(order="F", asmasked=False, fill_value=xtgeo.UNDEF)
 
     ier = _cxtgeo.surf_export_zmap_ascii(
         mfile,

src/xtgeo/surface/_regsurf_oper.py

@@ -5,6 +5,7 @@
 import numpy as np
 import numpy.ma as ma
 
+import xtgeo
 from xtgeo.xyz import Polygons
 import xtgeo.cxtgeo.cxtgeo as _cxtgeo
 from xtgeo.common import XTGeoDialog
@@ -129,7 +130,7 @@ def get_value_from_xy(self, point=(0.0, 0.0)):
         XTGDEBUG,
     )
 
-    if zcoord > self._undef_limit:
+    if zcoord > xtgeo.UNDEF_LIMIT:
         return None
 
     return zcoord
@@ -165,7 +166,7 @@ def get_xy_value_from_ij(self, iloc, jloc, zvalues=None):
     else:
         raise ValueError("Index i and/or j out of bounds")
 
-    if value > self.undef_limit:
+    if value > xtgeo.UNDEF_LIMIT:
         value = None
 
     return xval, yval, value
@@ -308,7 +309,7 @@ def get_fence(self, xyfence):
         logger.warning("Seem to be rotten")
 
     xyfence[:, 2] = czarr
-    xyfence = ma.masked_greater(xyfence, self._undef_limit)
+    xyfence = ma.masked_greater(xyfence, xtgeo.UNDEF_LIMIT)
     xyfence = ma.mask_rows(xyfence)
 
     return xyfence
@@ -325,7 +326,7 @@ def operation_polygons(self, poly, value, opname="add", inside=True):
 
     proxy = self.copy()
     proxy.values *= 0.0
-    vals = proxy.get_values1d(fill_value=self.undef)
+    vals = proxy.get_values1d(fill_value=xtgeo.UNDEF)
 
     # value could be a scalar or another surface; if another surface,
     # must ensure same topology
@@ -397,8 +398,8 @@ def operation_polygons(self, poly, value, opname="add", inside=True):
     elif opname == "set":
         tmp = value
     elif opname == "eli":
-        tmp = value * 0 + self.undef
-        tmp = ma.masked_greater(tmp, self.undef_limit)
+        tmp = value * 0 + xtgeo.UNDEF
+        tmp = ma.masked_greater(tmp, xtgeo.UNDEF_LIMIT)
 
     self.values[proxyv == proxytarget] = tmp[proxyv == proxytarget]
     del tmp

src/xtgeo/surface/_regsurf_roxapi.py

@@ -1,26 +1,23 @@
 # coding: utf-8
 """Roxar API functions for XTGeo RegularSurface"""
-import numpy as np
-
 from xtgeo.common import XTGeoDialog
+from xtgeo import RoxUtils
 
 xtg = XTGeoDialog()
 
 logger = xtg.functionlogger(__name__)
 
-# pylint: disable=protected-access
-
 
 def import_horizon_roxapi(self, project, name, category, stype, realisation):
     """Import a Horizon surface via ROXAR API spec."""
-    import roxar  # pylint: disable=import-error
 
-    if project is not None and isinstance(project, str):
-        projectname = project
-        with roxar.Project.open_import(projectname) as proj:
-            _roxapi_import_surface(self, proj, name, category, stype, realisation)
-    else:
-        _roxapi_import_surface(self, project, name, category, stype, realisation)
+    rox = RoxUtils(project, readonly=True)
+
+    proj = rox.project
+
+    _roxapi_import_surface(self, proj, name, category, stype, realisation)
+
+    rox.safe_close()
 
 
 def _roxapi_import_surface(self, proj, name, category, stype, realisation):
@@ -39,6 +36,7 @@ def _roxapi_import_surface(self, proj, name, category, stype, realisation):
             _roxapi_horizon_to_xtgeo(self, rox)
         except KeyError as kwe:
             logger.error(kwe)
+
     elif stype == "zones":
         if name not in proj.zones:
             raise ValueError("Name {} is not within Zones".format(name))
@@ -51,6 +49,19 @@ def _roxapi_import_surface(self, proj, name, category, stype, realisation):
             _roxapi_horizon_to_xtgeo(self, rox)
         except KeyError as kwe:
             logger.error(kwe)
+
+    elif stype == "clipboard":
+        if category:
+            if "|" in category:
+                folders = category.split("|")
+            else:
+                folders = category.split("/")
+            rox = proj.clipboard.folders[folders]
+        else:
+            rox = proj.clipboard
+        roxsurf = rox[name].get_grid(realisation)
+        _roxapi_horizon_to_xtgeo(self, roxsurf)
+
     else:
         raise ValueError("Invalid stype")
 
@@ -64,20 +75,19 @@ def _roxapi_horizon_to_xtgeo(self, rox):
     self._xinc, self._yinc = rox.increment
     self._rotation = rox.rotation
 
-    # since XTGeo is F order, while RMS is C order...
-    self._values = np.asanyarray(rox.get_values(), order="C")
+    self._values = rox.get_values()
+    logger.info("Surface from roxapi to xtgeo... DONE")
 
 
 def export_horizon_roxapi(self, project, name, category, stype, realisation):
     """Export (store) a Horizon surface to RMS via ROXAR API spec."""
-    import roxar  # pylint: disable=import-error
+    rox = RoxUtils(project, readonly=False)
 
-    if project is not None and isinstance(project, str):
-        projectname = project
-        with roxar.Project.open_import(projectname) as proj:
-            _roxapi_export_surface(self, proj, name, category, stype, realisation)
-    else:
-        _roxapi_export_surface(self, project, name, category, stype, realisation)
+    logger.info("Surface from xtgeo to roxapi...")
+    _roxapi_export_surface(self, rox.project, name, category, stype, realisation)
+
+    logger.info("Surface from xtgeo to roxapi... DONE")
+    rox.safe_close()
 
 
 def _roxapi_export_surface(self, proj, name, category, stype, realisation):
@@ -90,9 +100,9 @@ def _roxapi_export_surface(self, proj, name, category, stype, realisation):
             )
         try:
             roxroot = proj.horizons[name][category]
-            rox = _xtgeo_to_roxapi_grid(self)
-            rox.set_values(np.asanyarray(self.values, order="C"))
-            roxroot.set_grid(rox, realisation=realisation)
+            roxg = _xtgeo_to_roxapi_grid(self)
+            roxg.set_values(self.values)
+            roxroot.set_grid(roxg, realisation=realisation)
         except KeyError as kwe:
             logger.error(kwe)
 
@@ -105,12 +115,26 @@ def _roxapi_export_surface(self, proj, name, category, stype, realisation):
             )
         try:
             roxroot = proj.zones[name][category]
-            rox = _xtgeo_to_roxapi_grid(self)
-            rox.set_values(np.asanyarray(self.values, order="C"))
-            roxroot.set_grid(rox)
+            roxg = _xtgeo_to_roxapi_grid(self)
+            roxg.set_values(self.values)
+            roxroot.set_grid(roxg)
         except KeyError as kwe:
             logger.error(kwe)
 
+    elif stype == "clipboard":
+        folders = []
+        if category:
+            if "|" in category:
+                folders = category.split("|")
+            else:
+                folders = category.split("/")
+        if folders:
+            proj.clipboard.folders.create(folders)
+
+        roxroot = proj.clipboard.create_surface(name, folders)
+        roxg = _xtgeo_to_roxapi_grid(self)
+        roxg.set_values(self.values)
+        roxroot.set_grid(roxg)
     else:
         raise ValueError("Invalid stype")
 

src/xtgeo/surface/_regsurf_utils.py

@@ -1,5 +1,6 @@
 """RegularSurface utilities (basic low level)"""
 
+import xtgeo
 import xtgeo.cxtgeo.cxtgeo as _cxtgeo
 from xtgeo.common import XTGeoDialog
 
@@ -31,7 +32,7 @@ def swapaxes(self):
     _cxtgeo.doublepointer_assign(yinc, self._yinc)
     _cxtgeo.doublepointer_assign(rota, self._rotation)
 
-    val = self.get_values1d(fill_value=self.undef)
+    val = self.get_values1d(fill_value=xtgeo.UNDEF)
 
     ier = _cxtgeo.surf_swapaxes(
         ncol, nrow, yflip, self.xori, xinc, self.yori, yinc, rota, val, 0, XTGDEBUG
@@ -81,7 +82,7 @@ def swapaxes(self):
 #     xvv = np.reshape(xvv, (self._ncol, self._nrow), order='F')
 
 #     # make it masked
-#     xvv = ma.masked_greater(xvv, self._undef_limit)
+#     xvv = ma.masked_greater(xvv, xtgeo.UNDEF_LIMIT)
 
 #     self._values = xvv
 
@@ -109,7 +110,7 @@ def swapaxes(self):
 #         sys.exit(9)
 
 #     # make a 1D F order numpy array, and update C array
-#     xvv = ma.filled(self._values, self._undef)
+#     xvv = ma.filled(self._values, xtgeo.UNDEF)
 #     xvv = np.reshape(xvv, -1, order='F')
 
 #     self._cvalues = _cxtgeo.new_doublearray(nnum)