From c4cfe178efde45b1020ad5869dffaec27a2022c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Mon, 11 Sep 2023 15:30:19 +0200 Subject: [PATCH] ENH: allow RegularSurface() store np.float32 values As default, the numpy (masked) array is stored as 64 bit numpy, in line with Roxar API storage. Here it is allowed to cast to 32 bit in import. The methods surface_from_grid() and surface_from_cube() Still lack dtype, but it is easy to set with the property srf.dtype = ... --- src/xtgeo/surface/_regsurf_export.py | 14 ++-- src/xtgeo/surface/_regsurf_roxapi.py | 8 +- src/xtgeo/surface/regular_surface.py | 91 ++++++++++++++++++---- tests/test_roxarapi/test_roxarapi_reek.py | 16 ++++ tests/test_surface/test_regular_surface.py | 25 ++++++ 5 files changed, 129 insertions(+), 25 deletions(-) diff --git a/src/xtgeo/surface/_regsurf_export.py b/src/xtgeo/surface/_regsurf_export.py index 1091ca9b0..d7abd43e8 100644 --- a/src/xtgeo/surface/_regsurf_export.py +++ b/src/xtgeo/surface/_regsurf_export.py @@ -101,7 +101,7 @@ def _export_irap_ascii(self, mfile): self._xinc, self._yflip * self._yinc, self._rotation, - vals, + vals.astype(np.float64), 0, ) if ier != 0: @@ -199,7 +199,7 @@ def _export_irap_binary_cxtgeo(self, mfile): self._xinc, self._yflip * self._yinc, self._rotation, - vals, + vals.astype(np.float64), 0, ) @@ -225,7 +225,7 @@ def export_ijxyz_ascii(self, mfile): self._yflip, self._ilines, self._xlines, - vals, + vals.astype(np.float64), 0, ) @@ -324,7 +324,7 @@ def _export_zmap_ascii(self, mfile): scopy._yori, scopy._xinc, yinc, - vals, + vals.astype(np.float64), zmin, zmax, 0, @@ -361,7 +361,9 @@ def export_storm_binary(self, mfile): scopy._yori, scopy._xinc, yinc, - scopy.get_values1d(order="F", asmasked=False, fill_value=self.undef), + scopy.get_values1d(order="F", asmasked=False, fill_value=self.undef).astype( + np.float64 + ), zmin, zmax, 0, @@ -416,7 +418,7 @@ def export_petromod_binary(self, mfile, pmd_dataunits): _cxtgeo.surf_export_petromod_bin( mfile.get_cfhandle(), dsc, - values, + values.astype(np.float64), ) mfile.cfclose() diff --git a/src/xtgeo/surface/_regsurf_roxapi.py b/src/xtgeo/surface/_regsurf_roxapi.py index d0c6665f4..2adf06bd2 100644 --- a/src/xtgeo/surface/_regsurf_roxapi.py +++ b/src/xtgeo/surface/_regsurf_roxapi.py @@ -3,6 +3,8 @@ import os import tempfile +import numpy as np + from xtgeo import RoxUtils from xtgeo.common import XTGeoDialog @@ -152,7 +154,7 @@ def _roxapi_export_surface( try: roxroot = proj.horizons[name][category] roxg = _xtgeo_to_roxapi_grid(self) - roxg.set_values(self.values) + roxg.set_values(self.values.astype(np.float64)) roxroot.set_grid(roxg, realisation=realisation) except KeyError as kwe: logger.error(kwe) @@ -165,7 +167,7 @@ def _roxapi_export_surface( try: roxroot = proj.zones[name][category] roxg = _xtgeo_to_roxapi_grid(self) - roxg.set_values(self.values) + roxg.set_values(self.values.astype(np.float64)) roxroot.set_grid(roxg) except KeyError as kwe: logger.error(kwe) @@ -183,7 +185,7 @@ def _roxapi_export_surface( roxroot = styperef.create_surface(name, folders) roxg = _xtgeo_to_roxapi_grid(self) - roxg.set_values(self.values) + roxg.set_values(self.values.astype(np.float64)) roxroot.set_grid(roxg) elif stype == "trends": diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index 01dfc680b..b8a6792e2 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/src/xtgeo/surface/regular_surface.py @@ -74,7 +74,14 @@ # METHODS as wrappers to class init + import -def surface_from_file(mfile, fformat=None, template=None, values=True, engine="cxtgeo"): +def surface_from_file( + mfile, + fformat=None, + template=None, + values=True, + engine: Optional[str] = "cxtgeo", + dtype: Optional[Union[np.float64, np.float32]] = np.float64, +): """Make an instance of a RegularSurface directly from file import. Args: @@ -87,6 +94,7 @@ def surface_from_file(mfile, fformat=None, template=None, values=True, engine="c values (bool): If True (default), surface values will be read (Irap binary only) engine (str): Some import methods are implemnted in both C and Python. The C method ``cxtgeo`` is default. Alternative use ``python`` + dtype: Requsted numpy dtype of values; default is float64, alternatively float32 Example:: @@ -100,11 +108,23 @@ def surface_from_file(mfile, fformat=None, template=None, values=True, engine="c """ return RegularSurface._read_file( - mfile, fformat=fformat, load_values=values, engine=engine, template=template + mfile, + fformat=fformat, + load_values=values, + engine=engine, + template=template, + dtype=dtype, ) -def surface_from_roxar(project, name, category, stype="horizons", realisation=0): +def surface_from_roxar( + project, + name, + category, + stype="horizons", + realisation=0, + dtype: Optional[Union[np.float64, np.float32]] = np.float64, +): """This makes an instance of a RegularSurface directly from roxar input. Args: @@ -118,6 +138,7 @@ def surface_from_roxar(project, name, category, stype="horizons", realisation=0) stype (str): RMS folder type, 'horizons' (default), 'zones', 'clipboard', 'general2d_data' or 'trends' realisation (int): Realisation number, default is 0 + dtype: Requested numpy dtype for array; default is 64 bit Example:: @@ -136,7 +157,7 @@ def surface_from_roxar(project, name, category, stype="horizons", realisation=0) """ return RegularSurface._read_roxar( - project, name, category, stype=stype, realisation=realisation + project, name, category, stype=stype, realisation=realisation, dtype=dtype ) @@ -326,6 +347,7 @@ def __init__( filesrc: Optional[str] = None, fformat: Optional[str] = None, undef: Optional[float] = xtgeo.UNDEF, + dtype: Optional[Union[np.float64, np.float32]] = np.float64, ): """Instantiating a RegularSurface:: @@ -388,12 +410,13 @@ def __init__( self._metadata = xtgeo.MetaDataRegularSurface() self._values = None + if values is None: values = np.ma.zeros((self._ncol, self._nrow)) self._isloaded = False else: self._isloaded = True - self.values = values + self._ensure_correct_values(values, force_dtype=dtype) if ilines is None: self._ilines = np.array(range(1, self._ncol + 1), dtype=np.int32) @@ -685,9 +708,28 @@ def ymax(self): ymax = corner[1] return ymax + @property + def dtype(self) -> np.dtype: + """Getting the dtype of the values (np.array); float64 or float32""" + # this is not stored as state varible, but derived from the actual values + try: + infer_dtype = self._values.dtype + except AttributeError: + infer_dtype = np.float64 + return infer_dtype + + @dtype.setter + def dtype(self, wanted_dtype: Union[np.float64, np.float32]): + """Setting the dtype of the values (np.array); float64 or float32""" + try: + self._values = self._values.astype(wanted_dtype) + except TypeError as msg: + warnings.warn(f"Cannot change dtype: {msg}. Will keep current", UserWarning) + return + @property def values(self): - """The map values, as 2D masked numpy (float64), shape (ncol, nrow). + """The map values, as 2D masked numpy (float64/32), shape (ncol, nrow). When setting values as a scalar, the current mask will be preserved. @@ -985,10 +1027,11 @@ def _read_file( fformat = mfile.detect_fformat() else: fformat = mfile.generic_format_by_proposal(fformat) # default - kwargs = _data_reader_factory(fformat)(mfile, values=load_values, **kwargs) - kwargs["filesrc"] = mfile.file - kwargs["fformat"] = fformat - return cls(**kwargs) + new_kwargs = _data_reader_factory(fformat)(mfile, values=load_values, **kwargs) + new_kwargs["filesrc"] = mfile.file + new_kwargs["fformat"] = fformat + new_kwargs["dtype"] = kwargs.get("dtype", np.float64) + return cls(**new_kwargs) def load_values(self): """Import surface values in cases where metadata only is loaded. @@ -1200,7 +1243,13 @@ def to_hdf( @classmethod def _read_roxar( - cls, project, name, category, stype="horizons", realisation=0 + cls, + project, + name, + category, + stype="horizons", + realisation=0, + dtype=np.float64, ): # pragma: no cover """Load a surface from a Roxar RMS project. @@ -1223,11 +1272,14 @@ def _read_roxar( example 'DS_extracted' stype (str): RMS folder type, 'horizons' (default), 'zones' or 'clipboard' realisation (int): Realisation number, default is 0 + dtype: For supporting conversion to 32 bit float for the numpy; default + is 64 bit """ kwargs = _regsurf_roxapi.import_horizon_roxapi( project, name, category, stype, realisation ) + kwargs["dtype"] = dtype # eventual dtype change will be done in __init__ return cls(**kwargs) @@ -3019,7 +3071,7 @@ def translate_coordinates(self, translate=(0, 0, 0)): # ================================================================================== def _ensure_correct_values( - self, values + self, values, force_dtype=None ): # pylint: disable=too-many-branches, too-many-statements """Ensures that values is a 2D masked numpy (ncol, nrow), C order. @@ -3027,11 +3079,15 @@ def _ensure_correct_values( Args: values (array-like or scalar): Values to process. + force_dtype (numpy dtype or None): If not None, try to derive dtype from + current values Return: Nothing, self._values will be updated inplace """ + apply_dtype = force_dtype if force_dtype else self.dtype + currentmask = None if isinstance(self.values, ma.MaskedArray): currentmask = ma.getmaskarray(self.values) @@ -3060,7 +3116,7 @@ def _ensure_correct_values( elif isinstance(values, numbers.Number): if currentmask is not None: - vals = np.ones(self.dimensions, dtype=np.float64) * values + vals = np.ones(self.dimensions, dtype=apply_dtype) * values vals = np.ma.array(vals, mask=currentmask) # there maybe cases where values scalar input is some kind of UNDEF @@ -3070,11 +3126,11 @@ def _ensure_correct_values( self._values *= 0 self._values += vals else: - vals = ma.zeros((self.ncol, self.nrow), order="C", dtype=np.float64) + vals = ma.zeros((self.ncol, self.nrow), order="C", dtype=apply_dtype) self._values = vals + float(values) elif isinstance(values, (list, tuple, np.ndarray)): # ie values ~ list-like - vals = ma.array(values, order="C", dtype=np.float64) + vals = ma.array(values, order="C", dtype=apply_dtype) vals = ma.masked_greater(vals, self.undef_limit, copy=True) vals = ma.masked_invalid(vals, copy=True) @@ -3087,7 +3143,7 @@ def _ensure_correct_values( self._values = vals elif isinstance(values, (list, tuple)): # ie values ~ list-like - vals = ma.array(values, order="C", dtype=np.float64) + vals = ma.array(values, order="C", dtype=apply_dtype) vals = ma.masked_greater(vals, self.undef_limit, copy=True) vals = ma.masked_invalid(vals, copy=True) @@ -3104,3 +3160,6 @@ def _ensure_correct_values( if self._values.mask is ma.nomask: self._values = ma.array(self._values, mask=ma.getmaskarray(self._values)) + + # ensure dtype; avoid allocation and ID change if possible by setting copy=False + self._values = self._values.astype(apply_dtype, copy=False) diff --git a/tests/test_roxarapi/test_roxarapi_reek.py b/tests/test_roxarapi/test_roxarapi_reek.py index 701157681..cc5f29d11 100644 --- a/tests/test_roxarapi/test_roxarapi_reek.py +++ b/tests/test_roxarapi/test_roxarapi_reek.py @@ -171,6 +171,22 @@ def test_rox_surfaces(roxar_project): prj.close() +@pytest.mark.requires_roxar +def test_rox_surfaces_dtype_switching(roxar_project): + """Test dtype switching for from_roxar""" + srf = xtgeo.surface_from_roxar(roxar_project, "TopReek", SURFCAT1, dtype="float32") + assert srf.ncol == 554 + assert srf.values.mean() == pytest.approx(1698.648, abs=0.01) + assert srf.dtype == np.float32 + srf.to_roxar(roxar_project, "TopReek_copy", "SomeFolder", stype="clipboard") + + srf2 = srf.copy() + assert srf2.dtype == np.float32 + + srf2.dtype = np.float64 + np.testing.assert_allclose(srf.values, srf2.values) + + @pytest.mark.requires_roxar def test_rox_surfaces_alternative_open(roxar_project): """Based on previous but instead use a project ref as first argument""" diff --git a/tests/test_surface/test_regular_surface.py b/tests/test_surface/test_regular_surface.py index a381ff5d0..09f887552 100644 --- a/tests/test_surface/test_regular_surface.py +++ b/tests/test_surface/test_regular_surface.py @@ -181,6 +181,31 @@ def test_irapbin_import1(): xsurf.describe() +def test_irapbin_import_32bit(): + """Import Reek Irap binary, force 32 bit storage in import.""" + # setting dtype in importing surface + xsurf = xtgeo.surface_from_file(TESTSET2, dtype=np.float32) + assert xsurf.values.dtype == np.float32 + + # setting dtype explicit as 64 bit + xsurf.dtype = np.float64 + assert xsurf.values.dtype == np.float64 + + # setting dtype explicit as 32 bit + xsurf.dtype = np.float32 + assert xsurf.values.dtype == np.float32 + + xsurf.dtype = "float64" + assert xsurf.values.dtype == np.float64 + + +def test_irapbin_import_invalid_dtype(): + """Import Reek Irap binary, invalid dtype.""" + # setting dtype in importing surface + with pytest.raises(AttributeError): + xtgeo.surface_from_file(TESTSET2, dtype=np.float33) + + def test_irapbin_import_use_pathib(): """Import Reek Irap binary.""" logger.info("Import and export...")