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

@@ -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()

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":

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,19 +3071,23 @@ 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.
 
         This is an improved but private version over 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)

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

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