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...")