dev

davidhassell · davidhassell · commit 63ca9cec2db9 · 2025-08-05T13:25:57.000+01:00
diff --git a/cf/data/dask_utils.py b/cf/data/dask_utils.py
@@ -754,15 +754,18 @@ def cf_healpix_coordinates(
 
 
 def cf_healpix_increase_refinement(a, ncells, iaxis, quantity):
-    """Increase the refinement level of a HEALPix array.
+    """Increase the refinement level of a HEALPix data.
 
     Data are broadcast to cells of at the higher refinement level. For
-    an extensive field quantity (that depends on the size of the
-    cells, such as "sea_ice_mass" with units of kg), the new values
-    are reduced so that they are consistent with the new smaller cell
-    areas. For an intensive field quantity (that does not depend on
-    the size of the cells, such as "sea_ice_amount" with units of kg
-    m-2), the broadcast values are not changed.
+    an "extensive" quantity only, the broadcast values are reduced to
+    be consistent with the new smaller cell areas.
+
+    .. warning:: The returned `numpy` array will take up at least
+                 *ncells* times more memory than the input array
+                 *a*. For instance, if *a* has shape ``(3600, 48)``
+                 (1.3 MiB), the HEALPix axis is ``1``, and *ncells* is
+                 ``4096``, then the returned array will have shape
+                 ``(3600, 196608)`` (5.3 GiB), where 196608=48*4096.
 
     K. Gorski, Eric Hivon, A. Banday, B. Wandelt, M. Bartelmann, et
     al.. HEALPix: A Framework for High-Resolution Discretization and
@@ -781,7 +784,7 @@ def cf_healpix_increase_refinement(a, ncells, iaxis, quantity):
 
         ncells: `int`
             The number of cells at the new refinement level which are
-            contained in one cell at the original refinement level
+            contained in one cell at the original refinement level.
 
         iaxis: `int`
             The position of the HEALPix axis in the array dimensions.
@@ -819,7 +822,7 @@ def cf_healpix_increase_refinement(a, ncells, iaxis, quantity):
 
 
 def cf_healpix_increase_refinement_indices(a, ncells):
-    """Increase the refinement level of HEALPix indices
+    """Increase the refinement level of HEALPix indices.
 
     For instance, when going from refinement level 1 to refinement
     level 2, if *a* is ``(2, 23, 17)`` then it will be transformed to
@@ -828,6 +831,13 @@ def cf_healpix_increase_refinement_indices(a, ncells):
     ``ncells`` is the number of cells at refinement level 2 that lie
     inside one cell at refinement level 1, i.e. ``ncells=4**(2-1)=4``.
 
+    .. warning:: The returned `numpy` array will take up at least
+                 *ncells* times more memory than the input array
+                 *a*. For instance, if *a* has shape ``(48,)`` (384
+                 B), and *ncells* is ``262144``, then the returned
+                 array will have shape ``(12582912,)`` (96 MiB), where
+                 12582912=48*262144.
+
     K. Gorski, Eric Hivon, A. Banday, B. Wandelt, M. Bartelmann, et
     al.. HEALPix: A Framework for High-Resolution Discretization and
     Fast Analysis of Data Distributed on the Sphere. The Astrophysical
@@ -853,9 +863,6 @@ def cf_healpix_increase_refinement_indices(a, ncells):
             The array at the new refinement level.
 
     """
-    # PERFORMANCE: This function can use a lot of memory when 'a'
-    #              and/or 'ncells' are large.
-
     a = cfdm_to_memory(a)
 
     a = a * ncells
diff --git a/cf/domain.py b/cf/domain.py
@@ -415,10 +415,11 @@ def create_healpix(
             start = 0
 
         stop = start + ncells
-        data = Data(da.arange(start, stop), units="1")
+        dtype = cfdm.integer_dtype(stop - 1)
+        data = Data(da.arange(start, stop, dtype=dtype), units="1")
 
         # Set cached data elements
-        data._set_cached_elements({0: start, -1: stop - 1})
+        data._set_cached_elements({0: start, 1: start + 1, -1: stop - 1})
 
         c.set_data(data, copy=False)
         key = domain.set_construct(c, axes=axis, copy=False)
diff --git a/cf/field.py b/cf/field.py
@@ -4994,7 +4994,7 @@ def healpix_decrease_refinement_level(
 
         Set the refinement level to 0 using the ``'`range'`` *method*,
         which requires a new *reduction* function to be defined:
-        
+
         >>> import numpy as np
         >>> def range_func(a, axis=None):
         ...     return np.max(a, axis=axis) - np.min(a, axis=axis)
@@ -5114,16 +5114,17 @@ def healpix_decrease_refinement_level(
                 f"Can't decrease HEALPix refinement level of {f!r}: "
                 "'level' must be a non-negative integer less than "
                 "or equal to the current refinement level of "
-                f"{old_refinement_level}. Got {refinement_level!r}"
+                f"{old_refinement_level}. Got: {refinement_level!r}"
             )
 
+        # Parse 'level'
         if refinement_level == old_refinement_level:
             # No change in refinement level
             return f
 
         # Get the number of cells at the original refinement level
-        # which are contained in one cell at the lower refinement
-        # level
+        # which are contained in one larger cell at the new lower
+        # refinement level
         ncells = 4 ** (old_refinement_level - refinement_level)
 
         # Get the healpix_index coordinates
@@ -5260,18 +5261,20 @@ def healpix_decrease_refinement_level(
     def healpix_increase_refinement_level(self, refinement_level, quantity):
         """Increase the refinement level of a HEALPix grid.
 
-        Increasing the refinement level increases the resolution of
-        the HEALPix grid by broadcasting the Field data from each
-        original cell to all of the new smaller cells at the new
-        higher refinement level that lie inside it.
-
-        For an extensive field quantity (i.e. one that depends on the
-        size of the cells, such as "sea_ice_mass" with units of kg),
-        the broadcast values are also reduced to be consistent with
-        the new smaller cell areas. For an intensive field quantity
-        (i.e. one that does not depend on the size of the cells, such
-        as "sea_ice_amount" with units of kg m-2), the broadcast
-        values are not changed.
+        Increasing the HEALPix refinement level increases the
+        resolution of the HEALPix grid by broadcasting the Field data
+        from each original cell to all of the new smaller cells at the
+        new higher refinement level that lie inside it.
+
+        It must be specified whether the field data contains an
+        extensive or intensive quantity. An extensive quantity depends
+        on the size of the cells (such as "sea_ice_mass" with units of
+        kg, or "cell_area" with units of m2), and an intensive
+        quantity does not depend on the size of the cells (such as
+        "sea_ice_amount" with units of kg m-2, or "air_temperature"
+        with units of K). For an extensive quantity only, the
+        broadcast values are reduced to be consistent with the new
+        smaller cell areas.
 
         **References**
 
@@ -5332,12 +5335,12 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
         >>> print(g[0, :8] .array)
         [[291.5 291.5 291.5 291.5 293.5 293.5 293.5 293.5]]
 
-        For an extensive quantity (which ``f`` is not in this example,
-        but we can assume that it is for demonstration purposes), each
-        cell at the higher refinement level has the value of a cell at
-        the original refinement level after dividing it by the number
-        of cells at the higher refinement level that lie in one cell
-        of the original refinement level (4 in this case):
+        For an extensive quantity (which the ``f`` is this example is
+        not, but we can assume that it is for demonstration purposes),
+        each cell at the higher refinement level has the value of a
+        cell at the original refinement level after dividing it by the
+        number of cells at the higher refinement level that lie in one
+        cell of the original refinement level (4 in this case):
 
         >>> g = f.healpix_increase_refinement_level(2, 'extensive')
         >>> print(f[0, :2] .array)
@@ -5346,11 +5349,11 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
         [[72.875 72.875 72.875 72.875 73.375 73.375 73.375 73.375]]
 
         """
-        from .data.dask_utils import (
-            cf_healpix_increase_refinement,
-            cf_healpix_increase_refinement_indices,
+        from .data.dask_utils import cf_healpix_increase_refinement_indices
+        from .healpix import (
+            _healpix_increase_refinement_level,
+            healpix_max_refinement_level,
         )
-        from .healpix import healpix_max_refinement_level
 
         # Increasing the refinement level requires the nested indexing
         # scheme
@@ -5380,7 +5383,7 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
                 "'level' must be an integer greater than or equal to the "
                 f"current refinement level of {old_refinement_level}, and "
                 f"less than or equal to {healpix_max_refinement_level()}. "
-                f"Got {refinement_level!r}"
+                f"Got: {refinement_level!r}"
             )
 
         if refinement_level == old_refinement_level:
@@ -5393,7 +5396,7 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
             raise ValueError(
                 f"Can't increase HEALPix refinement level of {f!r}: "
                 f"'quantity' keyword must be one of {valid_quantities}. "
-                f"Got {quantity!r}"
+                f"Got: {quantity!r}"
             )
 
         # Get the number of cells at the higher refinement level which
@@ -5424,97 +5427,35 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
             )
         )
 
-        # Increase the refinement of the Field data
-        dx = f.data.to_dask_array(
-            _force_mask_hardness=False, _force_to_memory=False
-        )
-
-        if quantity == "extensive":
-            # Extensive data get divided by 'ncells' in
-            # `cf_healpix_increase_refinement`, so they end up with a
-            # data type of float64.
-            dtype = np.dtype("float64")
-        else:
-            dtype = dx.dtype
-
-        # Each chunk is going to get larger by a factor of 'ncells'
-        chunks = list(dx.chunks)
-        chunks[iaxis] = (np.array(chunks[iaxis]) * ncells).tolist()
-
-        dx = dx.map_blocks(
-            cf_healpix_increase_refinement,
-            chunks=tuple(chunks),
-            dtype=dtype,
-            meta=np.array((), dtype=dtype),
-            ncells=ncells,
-            iaxis=iaxis,
-            quantity=quantity,
-        )
-
         # Re-size the HEALPix axis
         domain_axis = f.domain_axis(axis)
-        domain_axis.set_size(dx.shape[iaxis])
+        domain_axis.set_size(f.shape[iaxis] * ncells)
 
-        f.set_data(dx, copy=False)
+        # Increase the refinement of the Field data
+        _healpix_increase_refinement_level(f, ncells, iaxis, quantity)
 
         # Increase the refinement level of domain ancillary constructs
         # that span the HEALPix axis. We're assuming that domain
         # ancillary data are intensive (i.e. do not depend on the size
         # of the cell).
-        meta = np.array((), dtype=dtype)
         for key, domain_ancillary in f.domain_ancillaries(
             filter_by_axis=(axis,), axis_mode="and", todict=True
         ).items():
             iaxis = f.get_data_axes(key).index(axis)
-            dx = domain_ancillary.data.to_dask_array(
-                _force_mask_hardness=False, _force_to_memory=False
-            )
-            dtype = dx.dtype
-
-            # Each chunk is going to get larger by a factor of
-            # 'ncells'
-            chunks = list(dx.chunks)
-            chunks[iaxis] = (np.array(chunks[iaxis]) * ncells).tolist()
-
-            dx = dx.map_blocks(
-                cf_healpix_increase_refinement,
-                chunks=tuple(chunks),
-                dtype=dtype,
-                meta=meta,
-                ncells=ncells,
-                iaxis=iaxis,
-                quantity="intensive",
+            _healpix_increase_refinement_level(
+                domain_ancillary, ncells, iaxis, "intensive"
             )
-            domain_ancillary.set_data(dx, copy=False)
 
         # Increase the refinement level of cell measure constructs
         # that span the HEALPix axis. Cell measure data are extensive
         # (i.e. depend on the size of the cell).
-        dtype = np.dtype("float64")
-        meta = np.array((), dtype=dtype)
         for key, cell_measure in f.cell_measures(
             filter_by_axis=(axis,), axis_mode="and", todict=True
         ).items():
             iaxis = f.get_data_axes(key).index(axis)
-            dx = cell_measure.data.to_dask_array(
-                _force_mask_hardness=False, _force_to_memory=False
-            )
-
-            # Each chunk is going to get larger by a factor of
-            # 'ncells'
-            chunks = list(dx.chunks)
-            chunks[iaxis] = (np.array(chunks[iaxis]) * ncells).tolist()
-
-            dx = dx.map_blocks(
-                cf_healpix_increase_refinement,
-                chunks=tuple(chunks),
-                dtype=dtype,
-                meta=meta,
-                ncells=ncells,
-                iaxis=iaxis,
-                quantity="extensive",
+            _healpix_increase_refinement_level(
+                cell_measure, ncells, iaxis, "extensive"
             )
-            cell_measure.set_data(dx, copy=False)
 
         # Remove all other metadata constructs that span the HEALPix
         # axis (including the original healpix_index coordinate
diff --git a/cf/healpix.py b/cf/healpix.py
@@ -157,8 +157,78 @@ def _healpix_create_latlon_coordinates(f, pole_longitude):
     return lat_key, lon_key
 
 
+def _healpix_increase_refinement_level(x, ncells, iaxis, quantity):
+    """Increase the HEALPix refinement level in-place.
+
+    K. Gorski, Eric Hivon, A. Banday, B. Wandelt, M. Bartelmann, et
+    al.. HEALPix: A Framework for High-Resolution Discretization and
+    Fast Analysis of Data Distributed on the Sphere. The Astrophysical
+    Journal, 2005, 622 (2), pp.759-771.
+    https://dx.doi.org/10.1086/427976
+
+    .. versionadded:: NEXTVERSION
+
+    .. seealso:: `cf.Field.healpix_increase_refinement_level`
+
+    :Parameters:
+
+        x: construct
+            The construct containing data that is to be changed.
+
+        ncells: `int`
+            The number of cells at the new refinement level which are
+            contained in one cell at the original refinement level.
+
+        iaxis: `int`
+            The position of the HEALPix axis in the construct's data
+            dimensions.
+
+        quantity: `str`
+            Whether the data represent intensive or extensive
+            quantities, specified with ``'intensive'`` and
+            ``'extensive'`` respectively.
+
+    :Returns:
+
+        `None`
+
+    """
+    from .data.dask_utils import cf_healpix_increase_refinement
+
+    # Get the Dask array.
+    #
+    # `cf_healpix_increase_refinement` has its own call to
+    # `cfdm_to_memory`, so we can set _force_to_memory=False.
+    dx = x.data.to_dask_array(
+        _force_mask_hardness=False, _force_to_memory=False
+    )
+
+    if quantity == "extensive":
+        # Extensive data get divided by 'ncells' in
+        # `cf_healpix_increase_refinement`, so they end up with a
+        # data type of float64.
+        dtype = np.dtype("float64")
+    else:
+        dtype = dx.dtype
+
+    # Each chunk is going to get larger by a factor of 'ncells'
+    chunks = list(dx.chunks)
+    chunks[iaxis] = (np.array(chunks[iaxis]) * ncells).tolist()
+
+    dx = dx.map_blocks(
+        cf_healpix_increase_refinement,
+        chunks=tuple(chunks),
+        dtype=dtype,
+        meta=np.array((), dtype=dtype),
+        ncells=ncells,
+        iaxis=iaxis,
+        quantity=quantity,
+    )
+    x.set_data(dx, copy=False)
+
+
 def _healpix_indexing_scheme(healpix_index, hp, new_indexing_scheme):
-    """Change the indexing scheme of HEALPix indices.
+    """Change the indexing scheme of HEALPix indices in-place.
 
     K. Gorski, Eric Hivon, A. Banday, B. Wandelt, M. Bartelmann, et
     al.. HEALPix: A Framework for High-Resolution Discretization and
diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py
