Repository: ecmwf/earthkit-data (python)

Issue:
Provide access to single precision GRIB values, latitudes and longitudes
Version 2.30.0 of ecCodes and version 1.6.0 of its Python binding (this latter is yet to be released) support data values array decoded as "floats" (single-precision) to reduce memory consumption. The C function looks like this:

```C
int codes_get_float_array(const codes_handle* h, const char* key, float* vals, size_t* length);
```

The properties/methods to consider on a GRIBField are as follows:

```python
m.values
m.to_numpy()
m.to_numpy(dtype=np.float32)
m.data()
m.to_latlon()
m.to_points()
```

At the moment all these methods return double precision values. The only exception is `to_numpy(dtype=np.float32)`, but it still gets double precision values from ecCodes and performs the cast on it.

The simplest solution is to add the `dtype` option to all the methods and ensure when it is set to `np.float32` they call `codes_get_float_array` and do not perform a cast on the result.

Of course, this would not work for the `values` property.

Pull request:
Feature/single precision grib values
Fixes #140 

Patch:
diff --git a/docs/examples/grib_lat_lon_value.ipynb b/docs/examples/grib_lat_lon_value.ipynb
index 560b573..eb79d52 100644
--- a/docs/examples/grib_lat_lon_value.ipynb
+++ b/docs/examples/grib_lat_lon_value.ipynb
@@ -13,7 +13,7 @@
    "id": "c74a86a6-79ae-44bc-a9c9-8cc8a3aaa484",
    "metadata": {},
    "source": [
-    "We read a GRIB file containing pressure level data. First we ensure the example file is available."
+    "In this example we will work with pressure level GRIB data. First we ensure the example file is available."
    ]
   },
   {
@@ -32,7 +32,7 @@
    "id": "ae4ff8c4-d87c-4d57-83a5-38e1e2545556",
    "metadata": {},
    "source": [
-    "We create a fieldlist with the temperature fields only."
+    "We will only use the temperature fields, so we extract them from the data."
    ]
   },
   {
@@ -198,7 +198,7 @@
    "id": "c7ec94f3-4e68-4380-a88d-128b2fb46cd8",
    "metadata": {},
    "source": [
-    "The simplest way to access the latitudes, longitudes and values is to use the *data()* method. It works both on fields and fieldlists."
+    "The simplest way to access the latitudes, longitudes and values is using the *data()* method. It works both on fields and fieldlists."
    ]
   },
   {
@@ -311,7 +311,7 @@
    "id": "a21f46b4-6c79-4306-a005-390ef2a02173",
    "metadata": {},
    "source": [
-    "*data()* only works on a fieldlist if all the fields has the same grid. The first two elements of the resulting ndarray are the latitude and longitude arrays (shared between fields), while the rest of the elements are the value arrays per field. Since we have six fields in our data the size of the first axis of the resulting ndarray is 2+6=8."
+    "*data()* only works on a fieldlist if all the fields have the same grid. The first two elements of the resulting ndarray are the latitude and longitude arrays (shared between fields), while the rest of the elements are the value arrays per field. Since we have six fields in our data the size of the first axis of the resulting ndarray is 2+6=8."
    ]
   },
   {
@@ -582,7 +582,7 @@
    "id": "392bdc90-0e15-415f-892b-7ebbf86cbc8c",
    "metadata": {},
    "source": [
-    "*to_latlon()* only works on a fieldlist if all the fields has the same grid. In this case it returns the same dict that we would get for any of the fields:"
+    "*to_latlon()* only works on a fieldlist if all the fields have the same grid. In this case it returns the same dict that we would get for any of the fields:"
    ]
   },
   {
@@ -679,7 +679,7 @@
    "id": "3bbbaad8-f9b7-4faa-b9ad-7a696f391b40",
    "metadata": {},
    "source": [
-    "For *to_numpy()* we can set the array type with the *dtype* keyword both for fields and fieldlists:"
+    "For all the methods above we can set the array type with the *dtype* keyword both for fields and fieldlists:"
    ]
   },
   {
@@ -708,11 +708,35 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 18,
    "id": "c89a473d-c6fd-4400-bfbc-da35bdcf8dc3",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "llv = ds.data(dtype=np.float32)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "id": "1e888bd9-48ec-4d18-8ef3-ee6d4fd107c0",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 90.     ,   0.     , 272.56418, 272.53918, 271.26532, 255.84306,\n",
+       "       244.00323, 226.65315], dtype=float32)"
+      ]
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "llv[:,0,0]"
+   ]
   }
  ],
  "metadata": {
diff --git a/earthkit/data/core/fieldlist.py b/earthkit/data/core/fieldlist.py
index e1140d3..f01dd9d 100644
--- a/earthkit/data/core/fieldlist.py
+++ b/earthkit/data/core/fieldlist.py
@@ -22,11 +22,29 @@ class Field(Base):
     def __init__(self, metadata=None):
         self.__metadata = metadata
 
-    @property
     @abstractmethod
+    def _values(self, dtype=None):
+        r"""Return the values stored in the field as a 1D ndarray.
+
+        Parameters
+        ----------
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
+
+        Returns
+        -------
+        ndarray
+            Field values
+
+        """
+        self._not_implemented()
+
+    @property
     def values(self):
         r"""ndarray: Get the values stored in the field as a 1D ndarray."""
-        self._not_implemented()
+        return self._values()
 
     def _make_metadata(self):
         r"""Create a field metadata object."""
@@ -47,8 +65,9 @@ class Field(Base):
         flatten: bool
             When it is True a flat ndarray is returned. Otherwise an ndarray with the field's
             :obj:`shape` is returned.
-        dtype: str or dtype
-            Typecode or data-type to which the array is cast.
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is ``np.float64``.
 
         Returns
         -------
@@ -58,12 +77,13 @@ class Field(Base):
         """
         values = self.values
         if not flatten:
-            values = self.values.reshape(self.shape)
+            # values = self.values.reshape(self.shape)
+            values = self._values(dtype=dtype).reshape(self.shape)
         if dtype is not None:
             values = values.astype(dtype)
         return values
 
-    def data(self, keys=("lat", "lon", "value"), flatten=False):
+    def data(self, keys=("lat", "lon", "value"), flatten=False, dtype=None):
         r"""Return the values and/or the geographical coordinates for each grid point.
 
         Parameters
@@ -74,6 +94,10 @@ class Field(Base):
         flatten: bool
             When it is True a flat ndarray per key is returned. Otherwise an ndarray with the field's
             :obj:`shape` is returned for each key.
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the arrays. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is ``np.float64``.
+
 
         Returns
         -------
@@ -115,7 +139,7 @@ class Field(Base):
         _keys = dict(
             lat=self._metadata.geography.latitudes,
             lon=self._metadata.geography.longitudes,
-            value=lambda: self.values,
+            value=self._values,
         )
 
         if isinstance(keys, str):
@@ -125,7 +149,7 @@ class Field(Base):
             if k not in _keys:
                 raise ValueError(f"data: invalid argument: {k}")
 
-        r = [_keys[k]() for k in keys]
+        r = [_keys[k](dtype=dtype) for k in keys]
         if not flatten:
             shape = self.shape
             r = [x.reshape(shape) for x in r]
@@ -136,7 +160,7 @@ class Field(Base):
 
             return np.array(r)
 
-    def to_points(self, flatten=False):
+    def to_points(self, flatten=False, dtype=None):
         r"""Return the geographical coordinates in the data's original
         Coordinate Reference System (CRS).
 
@@ -145,6 +169,10 @@ class Field(Base):
         flatten: bool
             When it is True 1D ndarrays are returned. Otherwise ndarrays with the field's
             :obj:`shape` are returned.
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the arrays. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
 
         Returns
         -------
@@ -162,8 +190,8 @@ class Field(Base):
         to_latlon
 
         """
-        x = self._metadata.geography.x()
-        y = self._metadata.geography.y()
+        x = self._metadata.geography.x(dtype=dtype)
+        y = self._metadata.geography.y(dtype=dtype)
         if x is not None and y is not None:
             if not flatten:
                 shape = self.shape
@@ -171,14 +199,14 @@ class Field(Base):
                 y = y.reshape(shape)
             return dict(x=x, y=y)
         elif self.projection().CARTOPY_CRS == "PlateCarree":
-            lon, lat = self.data(("lon", "lat"), flatten=flatten)
+            lon, lat = self.data(("lon", "lat"), flatten=flatten, dtype=dtype)
             return dict(x=lon, y=lat)
         else:
             raise ValueError(
                 "to_points(): geographical coordinates in original CRS are not available"
             )
 
-    def to_latlon(self, flatten=False):
+    def to_latlon(self, flatten=False, dtype=None):
         r"""Return the latitudes/longitudes of all the gridpoints in the field.
 
         Parameters
@@ -186,6 +214,10 @@ class Field(Base):
         flatten: bool
             When it is True 1D ndarrays are returned. Otherwise ndarrays with the field's
             :obj:`shape` are returned.
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the arrays. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
 
         Returns
         -------
@@ -198,7 +230,7 @@ class Field(Base):
         to_points
 
         """
-        lon, lat = self.data(("lon", "lat"), flatten=flatten)
+        lon, lat = self.data(("lon", "lat"), flatten=flatten, dtype=dtype)
         return dict(lat=lat, lon=lon)
 
     @property
@@ -544,7 +576,7 @@ class FieldList(Index):
         -------
         dict
             Unique, sorted metadata values from all the
-            :obj:`GribField <data.readers.grib.codes.GribField>`\ s.
+            :obj:`Field`\ s.
 
         See Also
         --------
@@ -590,7 +622,7 @@ class FieldList(Index):
         -------
         list
             Unique, sorted values of ``key`` from all the
-            :obj:`GribField <data.readers.grib.codes.GribField>`\ s.
+            :obj:`Field`\ s.
 
         See Also
         --------
@@ -660,7 +692,7 @@ class FieldList(Index):
 
         return np.array([f.values for f in self])
 
-    def data(self, keys=("lat", "lon", "value"), flatten=False):
+    def data(self, keys=("lat", "lon", "value"), flatten=False, dtype=None):
         r"""Return the values and/or the geographical coordinates.
 
         Only works when all the fields have the same grid geometry.
@@ -673,6 +705,10 @@ class FieldList(Index):
         flatten: bool
             When it is True the "lat", "lon" arrays and the "value" arrays per field
             will all be flattened. Otherwise they will preserve the field's :obj:`shape`.
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the arrays. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
 
         Returns
         -------
@@ -729,7 +765,7 @@ class FieldList(Index):
                 keys = [keys]
 
             if "lat" in keys or "lon" in keys:
-                latlon = self[0].to_latlon(flatten=flatten)
+                latlon = self[0].to_latlon(flatten=flatten, dtype=dtype)
 
             r = []
             for k in keys:
@@ -738,7 +774,7 @@ class FieldList(Index):
                 elif k == "lon":
                     r.append(latlon["lon"])
                 elif k == "value":
-                    r.extend([f.to_numpy(flatten=flatten) for f in self])
+                    r.extend([f.to_numpy(flatten=flatten, dtype=dtype) for f in self])
                 else:
                     raise ValueError(f"data: invalid argument: {k}")
 
@@ -977,7 +1013,7 @@ class FieldList(Index):
         ----------
         **kwargs: dict, optional
             Keyword arguments passed to
-            :obj:`GribField.to_points() <data.readers.grib.codes.GribField.to_points>`
+            :obj:`Field.to_points() <data.core.fieldlist.Field.to_points>`
 
         Returns
         -------
@@ -1004,7 +1040,7 @@ class FieldList(Index):
         ----------
         **kwargs: dict, optional
             Keyword arguments passed to
-            :obj:`GribField.to_latlon() <data.readers.grib.codes.GribField.to_latlon>`
+            :meth:`Field.to_latlon() <data.core.fieldlist.Field.to_latlon>`
 
         Returns
         -------
@@ -1094,7 +1130,7 @@ class FieldList(Index):
         -------
         list
             List with one :obj:`BoundingBox <data.utils.bbox.BoundingBox>` per
-            :obj:`GribField <data.readers.grib.codes.GribField>`
+            :obj:`Field`
         """
         return [s.bounding_box() for s in self]
 
diff --git a/earthkit/data/core/geography.py b/earthkit/data/core/geography.py
index 75788ee..7e9db2b 100644
--- a/earthkit/data/core/geography.py
+++ b/earthkit/data/core/geography.py
@@ -15,9 +15,16 @@ class Geography(metaclass=ABCMeta):
     """Geographical information about a field or data unit"""
 
     @abstractmethod
-    def latitudes(self):
+    def latitudes(self, dtype=None):
         r"""Return the latitudes.
 
+        Parameters
+        ----------
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
+
         Returns
         -------
         ndarray
@@ -25,9 +32,16 @@ class Geography(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def longitudes(self):
+    def longitudes(self, dtype=None):
         r"""Return the longitudes.
 
+        Parameters
+        ----------
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
+
         Returns
         -------
         ndarray
@@ -35,9 +49,16 @@ class Geography(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def x(self):
+    def x(self, dtype=None):
         r"""Return the x coordinates in the original CRS.
 
+        Parameters
+        ----------
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
+
         Returns
         -------
         ndarray
@@ -45,9 +66,16 @@ class Geography(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def y(self):
+    def y(self, dtype=None):
         r"""Return the y coordinates in the original CRS.
 
+        Parameters
+        ----------
+        dtype: str, numpy.dtype or None
+            Typecode or data-type of the array. When it is :obj:`None` the default
+            type used by the underlying data accessor is used. For GRIB it is
+            ``np.float64``.
+
         Returns
         -------
         ndarray
diff --git a/earthkit/data/readers/grib/codes.py b/earthkit/data/readers/grib/codes.py
index eb876fe..13e5cc9 100644
--- a/earthkit/data/readers/grib/codes.py
+++ b/earthkit/data/readers/grib/codes.py
@@ -33,6 +33,60 @@ def missing_is_none(x):
     return None if x == 2147483647 else x
 
 
+class GribCodesFloatArrayAccessor:
+    HAS_FLOAT_SUPPORT = None
+    KEY = None
+
+    def __init__(self):
+        if GribCodesFloatArrayAccessor.HAS_FLOAT_SUPPORT is None:
+            GribCodesFloatArrayAccessor.HAS_FLOAT_SUPPORT = hasattr(
+                eccodes, "codes_get_float_array"
+            )
+
+    def get(self, handle, dtype=None):
+        if dtype is np.float32 and self.HAS_FLOAT_SUPPORT:
+            return eccodes.codes_get_array(handle, self.KEY, ktype=dtype)
+        else:
+            return eccodes.codes_get_array(handle, self.KEY)
+
+
+class GribCodesValueAccessor(GribCodesFloatArrayAccessor):
+    def __init__(self):
+        super().__init__()
+        self.KEY = "values"
+
+
+class GribCodesLatitudeAccessor(GribCodesFloatArrayAccessor):
+    def __init__(self):
+        super().__init__()
+        self.KEY = "latitudes"
+
+    def get(self, handle, dtype=None):
+        v = eccodes.codes_get_array(handle, self.KEY)
+        if dtype is not None:
+            return v.astype(dtype)
+        else:
+            return v
+
+
+class GribCodesLongitudeAccessor(GribCodesFloatArrayAccessor):
+    def __init__(self):
+        super().__init__()
+        self.KEY = "longitudes"
+
+    def get(self, handle, dtype=None):
+        v = eccodes.codes_get_array(handle, self.KEY)
+        if dtype is not None:
+            return v.astype(dtype)
+        else:
+            return v
+
+
+VALUE_ACCESSOR = GribCodesValueAccessor()
+LATITUDE_ACCESSOR = GribCodesLatitudeAccessor()
+LONGITUDE_ACCESSOR = GribCodesLongitudeAccessor()
+
+
 class GribCodesMessagePositionIndex(CodesMessagePositionIndex):
     # This does not belong here, should be in the C library
     def _get_message_positions(self, path):
@@ -153,18 +207,18 @@ class GribCodesHandle(CodesHandle):
 
     # TODO: once missing value handling is implemented in the base class this method
     # can be removed
-    def get_values(self):
+    def get_values(self, dtype=None):
         eccodes.codes_set(self._handle, "missingValue", CodesHandle.MISSING_VALUE)
-        vals = eccodes.codes_get_values(self._handle)
+        vals = VALUE_ACCESSOR.get(self._handle, dtype=dtype)
         if self.get_long("bitmapPresent"):
             vals[vals == CodesHandle.MISSING_VALUE] = np.nan
         return vals
 
-    def get_latitudes(self):
-        return self.get("latitudes")
+    def get_latitudes(self, dtype=None):
+        return LATITUDE_ACCESSOR.get(self._handle, dtype=dtype)
 
-    def get_longitudes(self):
-        return self.get("longitudes")
+    def get_longitudes(self, dtype=None):
+        return LONGITUDE_ACCESSOR.get(self._handle, dtype=dtype)
 
     def get_data_points(self):
         return eccodes.codes_grib_get_data(self._handle)
@@ -214,10 +268,13 @@ class GribField(Field):
             self._handle = GribCodesReader.from_cache(self.path).at_offset(self._offset)
         return self._handle
 
-    @property
-    def values(self):
-        r"""ndarray: Gets the values stored in the GRIB field as a 1D ndarray."""
-        return self.handle.get_values()
+    def _values(self, dtype=None):
+        return self.handle.get_values(dtype=dtype)
+
+    # @property
+    # def values(self):
+    #     r"""ndarray: Gets the values stored in the GRIB field as a 1D ndarray."""
+    #     return self.handle.get_values()
 
     @property
     def offset(self):
diff --git a/earthkit/data/readers/grib/metadata.py b/earthkit/data/readers/grib/metadata.py
index 1b5520a..1441150 100644
--- a/earthkit/data/readers/grib/metadata.py
+++ b/earthkit/data/readers/grib/metadata.py
@@ -24,25 +24,25 @@ class GribFieldGeography(Geography):
     def __init__(self, metadata):
         self.metadata = metadata
 
-    def latitudes(self):
+    def latitudes(self, dtype=None):
         r"""Return the latitudes of the field.
 
         Returns
         -------
         ndarray
         """
-        return self.metadata._handle.get_latitudes()
+        return self.metadata._handle.get_latitudes(dtype=dtype)
 
-    def longitudes(self):
+    def longitudes(self, dtype=None):
         r"""Return the longitudes of the field.
 
         Returns
         -------
         ndarray
         """
-        return self.metadata._handle.get_longitudes()
+        return self.metadata._handle.get_longitudes(dtype=dtype)
 
-    def x(self):
+    def x(self, dtype=None):
         r"""Return the x coordinates in the field's original CRS.
 
         Returns
@@ -51,9 +51,9 @@ class GribFieldGeography(Geography):
         """
         grid_type = self.metadata.get("gridType", None)
         if grid_type in ["regular_ll", "reduced_gg", "regular_gg"]:
-            return self.longitudes()
+            return self.longitudes(dtype=dtype)
 
-    def y(self):
+    def y(self, dtype=None):
         r"""Return the y coordinates in the field's original CRS.
 
         Returns
@@ -62,7 +62,7 @@ class GribFieldGeography(Geography):
         """
         grid_type = self.metadata.get("gridType", None)
         if grid_type in ["regular_ll", "reduced_gg", "regular_gg"]:
-            return self.latitudes()
+            return self.latitudes(dtype=dtype)
 
     def shape(self):
         r"""Get the shape of the field.
diff --git a/earthkit/data/readers/netcdf.py b/earthkit/data/readers/netcdf.py
index 20c7cac..0074ff2 100644
--- a/earthkit/data/readers/netcdf.py
+++ b/earthkit/data/readers/netcdf.py
@@ -151,13 +151,13 @@ class XArrayFieldGeography(Geography):
         self._ds = ds
         self.north, self.west, self.south, self.east = self._ds.bbox(variable)
 
-    def latitudes(self):
-        return self.x()
+    def latitudes(self, dtype=None):
+        return self.x(dtype=dtype)
 
-    def longitudes(self):
-        return self.y()
+    def longitudes(self, dtype=None):
+        return self.y(dtype=dtype)
 
-    def _get_xy(self, axis, flatten=False):
+    def _get_xy(self, axis, flatten=False, dtype=None):
         if axis not in ("x", "y"):
             raise ValueError(f"Invalid axis={axis}")
 
@@ -177,13 +177,16 @@ class XArrayFieldGeography(Geography):
         points["x"], points["y"] = np.meshgrid(points["x"], points["y"])
         if flatten:
             points[axis] = points[axis].flatten()
-        return points[axis]
+        if dtype is not None:
+            return points[axis].astype(dtype)
+        else:
+            return points[axis]
 
-    def x(self):
-        return self._get_xy("x", flatten=True)
+    def x(self, dtype=None):
+        return self._get_xy("x", flatten=True, dtype=dtype)
 
-    def y(self):
-        return self._get_xy("y", flatten=True)
+    def y(self, dtype=None):
+        return self._get_xy("y", flatten=True, dtype=dtype)
 
     def shape(self):
         return self._da.shape[-2:]
@@ -311,8 +314,7 @@ class XArrayField(Field):
     def _to_numpy(self):
         return self.to_xarray().to_numpy()
 
-    @property
-    def values(self):
+    def _values(self, dtype=None):
         return self._to_numpy().flatten()
 
     def to_numpy(self, flatten=False, dtype=None):
diff --git a/earthkit/data/sources/list_of_dicts.py b/earthkit/data/sources/list_of_dicts.py
index b0b2be2..3417d6c 100644
--- a/earthkit/data/sources/list_of_dicts.py
+++ b/earthkit/data/sources/list_of_dicts.py
@@ -98,21 +98,29 @@ class VirtualGribMetadata(RawMetadata):
     def namespaces(self):
         return []
 
-    def latitudes(self):
-        return self.get("latitudes")
-
-    def longitudes(self):
-        return self.get("longitudes")
-
-    def x(self):
+    def latitudes(self, dtype=None):
+        v = self.get("latitudes")
+        if dtype is None:
+            return v
+        else:
+            return v.astype(dtype)
+
+    def longitudes(self, dtype=None):
+        v = self.get("longitudes")
+        if dtype is None:
+            return v
+        else:
+            return v.astype(dtype)
+
+    def x(self, dtype=None):
         grid_type = self.get("gridType", None)
         if grid_type in ["regular_ll", "reduced_gg", "regular_gg"]:
-            return self.longitudes()
+            return self.longitudes(dtype=dtype)
 
-    def y(self):
+    def y(self, dtype=None):
         grid_type = self.get("gridType", None)
         if grid_type in ["regular_ll", "reduced_gg", "regular_gg"]:
-            return self.latitudes()
+            return self.latitudes(dtype=dtype)
 
     def _unique_grid_id(self):
         return self.get("md5GridSection", None)
@@ -156,9 +164,12 @@ class VirtualGribField(Field):
     def __init__(self, d):
         super().__init__(metadata=VirtualGribMetadata(d))
 
-    @property
-    def values(self):
-        return self._metadata["values"]
+    def _values(self, dtype=None):
+        v = self._metadata["values"]
+        if dtype is None:
+            return v
+        else:
+            return v.astype(dtype)
 
     def _make_metadata(self):
         pass


Test patch:
diff --git a/tests/grib/test_grib_geography.py b/tests/grib/test_grib_geography.py
index 7d735f0..8e6795f 100644
--- a/tests/grib/test_grib_geography.py
+++ b/tests/grib/test_grib_geography.py
@@ -34,6 +34,8 @@ def test_grib_to_latlon_single(index):
     assert isinstance(v, dict)
     assert isinstance(v["lon"], np.ndarray)
     assert isinstance(v["lat"], np.ndarray)
+    assert v["lon"].dtype == np.float64
+    assert v["lat"].dtype == np.float64
     check_array(
         v["lon"],
         (84,),
@@ -64,26 +66,31 @@ def test_grib_to_latlon_single_shape(index):
 
     # x
     assert v["lon"].shape == (7, 12)
+    assert v["lon"].dtype == np.float64
     for x in v["lon"]:
         assert np.allclose(x, np.linspace(0, 330, 12))
 
     # y
     assert v["lat"].shape == (7, 12)
+    assert v["lon"].dtype == np.float64
     for i, y in enumerate(v["lat"]):
         assert np.allclose(y, np.ones(12) * (90 - i * 30))
 
 
-def test_grib_to_latlon_multi():
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_grib_to_latlon_multi(dtype):
     f = from_source("file", earthkit_examples_file("test.grib"))
 
-    v_ref = f[0].to_latlon(flatten=True)
-    v = f.to_latlon(flatten=True)
+    v_ref = f[0].to_latlon(flatten=True, dtype=dtype)
+    v = f.to_latlon(flatten=True, dtype=dtype)
     assert isinstance(v, dict)
     assert v.keys() == v_ref.keys()
 
     assert isinstance(v, dict)
     assert np.allclose(v["lat"], v_ref["lat"])
     assert np.allclose(v["lon"], v_ref["lon"])
+    assert v["lat"].dtype == dtype
+    assert v["lon"].dtype == dtype
 
 
 def test_grib_to_latlon_multi_non_shared_grid():
@@ -105,6 +112,8 @@ def test_grib_to_points_single(index):
     assert isinstance(v, dict)
     assert isinstance(v["x"], np.ndarray)
     assert isinstance(v["y"], np.ndarray)
+    assert v["x"].dtype == np.float64
+    assert v["y"].dtype == np.float64
     check_array(
         v["x"],
         (84,),
@@ -129,17 +138,20 @@ def test_grib_to_points_unsupported_grid():
         f[0].to_points()
 
 
-def test_grib_to_points_multi():
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_grib_to_points_multi(dtype):
     f = from_source("file", earthkit_examples_file("test.grib"))
 
-    v_ref = f[0].to_points(flatten=True)
-    v = f.to_points(flatten=True)
+    v_ref = f[0].to_points(flatten=True, dtype=dtype)
+    v = f.to_points(flatten=True, dtype=dtype)
     assert isinstance(v, dict)
     assert v.keys() == v_ref.keys()
 
     assert isinstance(v, dict)
     assert np.allclose(v["x"], v_ref["x"])
     assert np.allclose(v["y"], v_ref["y"])
+    assert v["x"].dtype == dtype
+    assert v["y"].dtype == dtype
 
 
 def test_grib_to_points_multi_non_shared_grid():
diff --git a/tests/grib/test_grib_values.py b/tests/grib/test_grib_values.py
index 3da1749..57ef274 100644
--- a/tests/grib/test_grib_values.py
+++ b/tests/grib/test_grib_values.py
@@ -31,6 +31,7 @@ def test_grib_values_1():
     # whole file
     v = f.values
     assert isinstance(v, np.ndarray)
+    assert v.dtype == np.float64
     assert v.shape == (1, 84)
     v = v[0].flatten()
     check_array(
@@ -57,6 +58,7 @@ def test_grib_values_18():
     # whole file
     v = f.values
     assert isinstance(v, np.ndarray)
+    assert v.dtype == np.float64
     assert v.shape == (18, 84)
     vf = v[0].flatten()
     check_array(
@@ -85,6 +87,7 @@ def test_grib_to_numpy_1():
     eps = 1e-5
     v = f.to_numpy()
     assert isinstance(v, np.ndarray)
+    assert v.dtype == np.float64
     v = v[0].flatten()
     check_array(
         v,
@@ -116,6 +119,7 @@ def test_grib_to_numpy_1_shape(first, options, expected_shape):
     data = f[0] if first else f
     v1 = data.to_numpy(**options)
     assert isinstance(v1, np.ndarray)
+    assert v1.dtype == np.float64
     assert v1.shape == expected_shape
     v1 = v1.flatten()
     assert np.allclose(v_ref, v1, eps)
@@ -129,6 +133,7 @@ def test_grib_to_numpy_18():
     # whole file
     v = f.to_numpy(flatten=True)
     assert isinstance(v, np.ndarray)
+    assert v.dtype == np.float64
     assert v.shape == (18, 84)
     vf0 = v[0].flatten()
     check_array(
@@ -180,6 +185,7 @@ def test_grib_to_numpy_18_shape(options, expected_shape):
     # whole file
     v = f.to_numpy()
     assert isinstance(v, np.ndarray)
+    assert v.dtype == np.float64
     assert v.shape == (18, 7, 12)
     vf0 = f[0].to_numpy().flatten()
     assert vf0.shape == (84,)
@@ -188,6 +194,7 @@ def test_grib_to_numpy_18_shape(options, expected_shape):
 
     v1 = f.to_numpy(**options)
     assert isinstance(v1, np.ndarray)
+    assert v1.dtype == np.float64
     assert v1.shape == expected_shape
     vr = v1[0].flatten()
     assert np.allclose(vf0, vr, eps)
@@ -218,10 +225,18 @@ def test_grib_to_numpy_18_dtype(dtype):
 
 
 @pytest.mark.parametrize(
-    "kwarg,expected_shape",
-    [({}, (11, 19)), ({"flatten": True}, (209,)), ({"flatten": False}, (11, 19))],
+    "kwarg,expected_shape,expected_dtype",
+    [
+        ({}, (11, 19), np.float64),
+        ({"flatten": True}, (209,), np.float64),
+        ({"flatten": True, "dtype": np.float32}, (209,), np.float32),
+        ({"flatten": True, "dtype": np.float64}, (209,), np.float64),
+        ({"flatten": False}, (11, 19), np.float64),
+        ({"flatten": False, "dtype": np.float32}, (11, 19), np.float32),
+        ({"flatten": False, "dtype": np.float64}, (11, 19), np.float64),
+    ],
 )
-def test_grib_field_data(kwarg, expected_shape):
+def test_grib_field_data(kwarg, expected_shape, expected_dtype):
     ds = from_source("file", earthkit_examples_file("test.grib"))
 
     latlon = ds[0].to_latlon(**kwarg)
@@ -229,6 +244,7 @@ def test_grib_field_data(kwarg, expected_shape):
 
     d = ds[0].data(**kwarg)
     assert isinstance(d, np.ndarray)
+    assert d.dtype == expected_dtype
     assert len(d) == 3
     assert d[0].shape == expected_shape
     assert np.allclose(d[0], latlon["lat"])
@@ -237,35 +253,40 @@ def test_grib_field_data(kwarg, expected_shape):
 
     d = ds[0].data(keys="lat", **kwarg)
     assert d.shape == expected_shape
+    assert d.dtype == expected_dtype
     assert np.allclose(d, latlon["lat"])
 
     d = ds[0].data(keys="lon", **kwarg)
     assert d.shape == expected_shape
+    assert d.dtype == expected_dtype
     assert np.allclose(d, latlon["lon"])
 
     d = ds[0].data(keys="value", **kwarg)
     assert d.shape == expected_shape
+    assert d.dtype == expected_dtype
     assert np.allclose(d, v)
 
     d = ds[0].data(keys=("value", "lon"), **kwarg)
     assert isinstance(d, np.ndarray)
+    assert d.dtype == expected_dtype
     assert len(d) == 2
     assert np.allclose(d[0], v)
     assert np.allclose(d[1], latlon["lon"])
 
 
 @pytest.mark.parametrize(
-    "kwarg,expected_shape",
+    "kwarg,expected_shape,expected_dtype",
     [
-        ({}, (11, 19)),
-        (
-            {"flatten": True},
-            (209,),
-        ),
-        ({"flatten": False}, (11, 19)),
+        ({}, (11, 19), np.float64),
+        ({"flatten": True}, (209,), np.float64),
+        ({"flatten": True, "dtype": np.float32}, (209,), np.float32),
+        ({"flatten": True, "dtype": np.float64}, (209,), np.float64),
+        ({"flatten": False}, (11, 19), np.float64),
+        ({"flatten": False, "dtype": np.float32}, (11, 19), np.float32),
+        ({"flatten": False, "dtype": np.float64}, (11, 19), np.float64),
     ],
 )
-def test_grib_fieldlist_data(kwarg, expected_shape):
+def test_grib_fieldlist_data(kwarg, expected_shape, expected_dtype):
     ds = from_source("file", earthkit_examples_file("test.grib"))
 
     latlon = ds.to_latlon(**kwarg)
@@ -274,6 +295,7 @@ def test_grib_fieldlist_data(kwarg, expected_shape):
     d = ds.data(**kwarg)
     assert isinstance(d, np.ndarray)
     assert d.shape == tuple([4, *expected_shape])
+    assert d.dtype == expected_dtype
     assert np.allclose(d[0], latlon["lat"])
     assert np.allclose(d[1], latlon["lon"])
     assert np.allclose(d[2], v[0])
@@ -281,19 +303,23 @@ def test_grib_fieldlist_data(kwarg, expected_shape):
 
     d = ds.data(keys="lat", **kwarg)
     assert d.shape == tuple([1, *expected_shape])
+    assert d.dtype == expected_dtype
     assert np.allclose(d[0], latlon["lat"])
 
     d = ds.data(keys="lon", **kwarg)
     assert d.shape == tuple([1, *expected_shape])
+    assert d.dtype == expected_dtype
     assert np.allclose(d[0], latlon["lon"])
 
     d = ds.data(keys="value", **kwarg)
     assert d.shape == tuple([2, *expected_shape])
+    assert d.dtype == expected_dtype
     assert np.allclose(d, v)
 
     d = ds.data(keys=("value", "lon"), **kwarg)
     assert isinstance(d, np.ndarray)
     assert d.shape == tuple([3, *expected_shape])
+    assert d.dtype == expected_dtype
     assert np.allclose(d[0], v[0])
     assert np.allclose(d[1], v[1])
     assert np.allclose(d[2], latlon["lon"])
diff --git a/tests/netcdf/test_netcdf_geography.py b/tests/netcdf/test_netcdf_geography.py
index 3ad1636..ee79ea9 100644
--- a/tests/netcdf/test_netcdf_geography.py
+++ b/tests/netcdf/test_netcdf_geography.py
@@ -28,14 +28,20 @@ def check_array(v, shape=None, first=None, last=None, meanv=None, eps=1e-3):
     assert np.isclose(v.mean(), meanv, eps)
 
 
-def test_netcdf_to_points_1():
+@pytest.mark.parametrize(
+    "dtype,expected_dtype",
+    [(None, np.float64), (np.float32, np.float32), (np.float64, np.float64)],
+)
+def test_netcdf_to_points_1(dtype, expected_dtype):
     ds = from_source("file", earthkit_test_data_file("test_single.nc"))
 
     eps = 1e-5
-    v = ds[0].to_points(flatten=True)
+    v = ds[0].to_points(flatten=True, dtype=dtype)
     assert isinstance(v, dict)
     assert isinstance(v["x"], np.ndarray)
     assert isinstance(v["y"], np.ndarray)
+    assert v["x"].dtype == expected_dtype
+    assert v["y"].dtype == expected_dtype
     check_array(
         v["x"],
         (84,),