dev

Author: davidhassell

cf/data/dask_utils.py

@@ -571,25 +571,25 @@
                 return the `esmpy.Regrid` instance that defines the
                 regridding operation.""",
     # HEALPix indexing schemes
-    "{{HEALPix indexing schemes}}": """The "nested" scheme indexes the pixels inside a single
+    "{{HEALPix indexing schemes}}": """The nested scheme indexes the pixels inside a single
                 coarser refinement level cell with consecutive
-                indices. The "ring" scheme indexes the pixels moving
+                indices. The ring scheme indexes the pixels moving
                 down from the north to the south pole along each
                 isolatitude ring. When the HEALPix axis is ordered
                 with monotonically increasing indices, each type of
                 indexing scheme is optimised for different types of
-                operation. For instance, the "ring" scheme is
-                optimised for Fourier transforms with spherical
-                harmonics; and the "nested" scheme is optimised for
-                geographical nearest-neighbour operations such as
-                decreasing the refinement level.
+                operation. For instance, the ring scheme is optimised
+                for Fourier transforxms with spherical harmonics; and
+                the nested scheme is optimised for geographical
+                nearest-neighbour operations such as decreasing the
+                refinement level.
 
                 A Multi-Order Coverage (MOC) has pixels with different
                 refinement levels stored in the same array. An
                 indexing scheme for an MOC has a unique index for each
                 cell at each refinement level.
 
-                The "nuniq" scheme defines MOC indices such that all
+                The nuniq scheme defines MOC indices such that all
                 cells within a particular refinement level form a
                 sequence of consecutive integers. E.g. for refinement
                 level 0 the indices are 4, ..., 15, for refinement
@@ -599,7 +599,7 @@
                 when the indices are sorted monotonically, for
                 within-refinement level data retrievals.
 
-                The "zuniq" scheme defines MOC indices such that, for
+                The zuniq scheme defines MOC indices such that, for
                 similar refinement levels, cells in the proximity of a
                 particular geographical location have similar index
                 values. This means that, unlike the nuniq case, the

cf/docstring/docstring.py

@@ -282,9 +282,10 @@ def create_healpix(
         resulting from an attempt to create an excessive amount of
         chunks for the healpix_index coordinates. For instance,
         healpix_index coordinates at refinement level 29 would need
-        ~206 billion chunks with the default Dask chunksize of 128
-        MiB, but with a chunksize of 1 pebibyte only 24576 chunks are
-        required::
+        ~206 billion chunks with the default Dask chunksize of 128 MiB
+        - almost certainly more than enough to cause a crash - but
+        with a chunksize of 1 pebibyte only 24576 chunks are
+        required, a much more manageable amount::
 
            >>> cf.chunksize()
            >>> 134217728
@@ -409,6 +410,11 @@ def create_healpix(
                 f"of {healpix_indexing_schemes!r}. Got {indexing_scheme!r}"
             )
 
+        indexing_scheme0 = indexing_scheme
+
+        if indexing_scheme == "zuniq":
+            indexing_scheme = "nuniq"
+
         domain = Domain()
         ncells = 12 * (4**refinement_level)
 
@@ -465,9 +471,9 @@ def create_healpix(
 
         domain.set_construct(cr)
 
-        if indexing_scheme == "zuniq":
+        if indexing_scheme0 == "zuniq":
             domain = domain.healpix_indexing_scheme("zuniq", sort=False)
-            
+
         return domain
 
     @_inplace_enabled(default=False)

cf/domain.py

@@ -4824,7 +4824,7 @@ def healpix_decrease_refinement_level(
         refinement level either
 
         * contains no original cells, in which case that larger cell
-          is not included in the output,
+          is not included in the output;
 
         or
 
@@ -4868,20 +4868,24 @@ def healpix_decrease_refinement_level(
                 ``'root_mean_square'``, ``'standard_deviation'``,
                 ``'sum'``, ``'sum_of_squares'``, ``'variance'``.
 
-                The method should be appropriate to nature of the
-                Field quantity, which is either intensive (i.e. that
-                does not depend on the size of the cells, such as
-                "sea_ice_amount" with units of kg m-2), or extensive
-                (i.e. that depends on the size of the cells, such as
-                "sea_ice_mass" with units of kg).
+                .. note:: The method should be appropriate to nature
+                          of the Field quantity, which is either
+                          intensive (i.e. that does not depend on the
+                          size of the cells, such as "sea_ice_amount"
+                          with units of kg m-2), or extensive
+                          (i.e. that depends on the size of the cells,
+                          such as "sea_ice_mass" with units of kg).
 
             reduction: function or `None`, optional
                 The function used to calculate the values in the new
                 larger cells, from the data on the original cells. The
-                function must i) calculate the quantity defined by the
-                *method* parameter, ii) take an array of values as its
-                first argument, and iii) have an *axis* keyword that
-                specifies which of the array axes is the HEALPix axis.
+                function must:
+
+                * calculate the quantity defined by the *method*
+                  parameter,
+                * take an array of values as its first argument,
+                * have an *axis* keyword that specifies which of the
+                  array axes is the HEALPix axis.
 
                 For some methods there are default *reduction*
                 functions, which are only used when *reduction* is
@@ -4907,11 +4911,11 @@ def healpix_decrease_refinement_level(
                 If True (the default) then the HEALPix grid is
                 automatically converted to a form suitable for having
                 its refinement level changed, i.e. the indexing scheme
-                is changed to 'nested' and the HEALPix axis is sorted
-                so that the nested HEALPix indices are monotonically
-                increasing. If False then either an exception is
-                raised if the HEALPix indexing scheme is not already
-                'nested', or else the HEALPix axis is not sorted.
+                is changed to nested and the HEALPix axis is sorted so
+                that the nested HEALPix indices are monotonically
+                increasing. If False then an exception is raised if
+                the HEALPix indexing scheme is not already nested and
+                the HEALPix axis is not sorted.
 
                 .. note:: Setting to False will speed up the operation
                           when the HEALPix indexing scheme is already
@@ -4924,12 +4928,12 @@ def healpix_decrease_refinement_level(
                 (but after the HEALPix grid has been conformed, when
                 *conform* is True):
 
-                1. The nested HEALPix indices are strictly
-                   monotonically increasing.
+                * The nested HEALPix indices are strictly
+                  monotonically increasing.
 
-                2. Every cell at the new lower refinement level
-                   contains zero or the maximum possible number of
-                   cells at the original refinement level.
+                * Every cell at the new lower refinement level
+                  contains zero or the maximum possible number of
+                  cells at the original refinement level.
 
                 If False then these checks are not carried out.
 
@@ -4938,14 +4942,16 @@ def healpix_decrease_refinement_level(
                              advance that these conditions are
                              satisfied. If set to False and any of the
                              conditions are not met then either an
-                             exception may be raised or, **much
+                             exception will be raised or, **much
                              worse**, the operation could complete and
                              return incorrect data values.
 
         :Returns:
 
             `Field`
-                A new Field with the new HEALPix grid.
+                A new Field with the new HEALPix grid. The HEALPix
+                indices of this field will follow the nested indexing
+                scheme.
 
         **Examples**
 
@@ -5084,6 +5090,7 @@ def healpix_decrease_refinement_level(
 
             # Re-get the HEALPix info
             hp = f.healpix_info()
+
         elif indexing_scheme != "nested":
             raise ValueError(
                 f"Can't decrease HEALPix refinement level of {f!r}: "
@@ -5149,7 +5156,13 @@ def healpix_decrease_refinement_level(
 
         # Get the HEALPix axis
         axis = hp["domain_axis_key"]
-        iaxis = f.get_data_axes().index(axis)
+        try:
+            iaxis = f.get_data_axes().index(axis)
+        except ValueError:
+            # The Field data doesn't span the size 1 HEALPix axis, so
+            # insert it on the left hand side.
+            f.insert_dimmension(axis, -1, inplace=True)
+            iaxis = f.ndim - 1
 
         # Whether or not to create lat/lon coordinates for the new
         # refinement level. Only do so if the original grid has
@@ -5169,7 +5182,7 @@ def healpix_decrease_refinement_level(
         # ------------------------------------------------------------
 
         # Note: Using 'Data.coarsen' works because a) the HEALPix
-        #       indexing scheme is now "nested", and b) we know that
+        #       indexing scheme is now nested, and b) we know that
         #       each new coarser cell contains the maximum possible
         #       number (i.e. 'ncells') of original cells.
         f.data.coarsen(
@@ -5263,7 +5276,7 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
         extensive or intensive quantity. 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),
-        and an extensive quantity depends on the size of the cells
+        and an extensive quantity does depend on the size of the cells
         (such as "sea_ice_mass" with units of kg, or "cell_area" with
         units of m2). For an extensive quantity only, the broadcast
         values are reduced to be consistent with the new smaller cell
@@ -5294,7 +5307,9 @@ def healpix_increase_refinement_level(self, refinement_level, quantity):
         :Returns:
 
             `Field`
-                A new Field with the new HEALPix grid.
+                A new Field with the new HEALPix grid. The HEALPix
+                indices of this field will follow the nested indexing
+                scheme.
 
         **Examples**
 

cf/field.py

@@ -316,7 +316,7 @@ def _healpix_increase_refinement_level_indices(
 
         healpix_index: `Coordinate`
             The HEALPix indices to be changed. They must use the
-            "nested" indexing scheme.
+            nested indexing scheme, but this is not checked.
 
         ncells: `int`
             The number of cells at the new higher refinement level
@@ -398,9 +398,13 @@ def _healpix_increase_refinement_level_indices(
         data._set_cached_elements({-1: x})
 
 
-def _healpix_indexing_scheme(f, hp, new_indexing_scheme):
+def _healpix_indexing_scheme(
+    f, hp, new_indexing_scheme, moc_refinement_level=None
+):
     """Change the indexing scheme of HEALPix indices in-place.
 
+    The
+
     See CF Appendix F: Grid Mappings.
     https://doi.org/10.5281/zenodo.14274886
 
@@ -419,7 +423,11 @@ def _healpix_indexing_scheme(f, hp, new_indexing_scheme):
             The HEALPix info dictionary for *f*.
 
         new_indexing_scheme: `str`
-            The new indexing scheme.
+            The new indexing scheme. It is assumed to be different to
+            the original indexing scheme of *f*.
+
+        moc_refinement_level: `int` or `None`, optional
+            The unique refinement level of MOC indices of *f*.
 
     :Returns:
 
@@ -435,7 +443,7 @@ def _healpix_indexing_scheme(f, hp, new_indexing_scheme):
     refinement_level = hp.get("refinement_level")
 
     old_cache = healpix_index.data._get_cached_elements()
-    
+
     # Change the HEALPix indices
     #
     # `cf_healpix_indexing_scheme` has its own call to
@@ -447,21 +455,31 @@ def _healpix_indexing_scheme(f, hp, new_indexing_scheme):
     # Find the datatype for the largest possible index at this
     # refinement level
     if new_indexing_scheme == "nuniq":
-        # The largest possible "nuniq" index at refinement level N is
-        # 16*(4**N) - 1 = 4*(4**N) + 12*(4**N) - 1
-        dtype = integer_dtype(16 * (4**refinement_level) - 1)
+        if indexing_scheme in ("nested", "ring"):
+            # The largest possible nested or ring index at refinement
+            # level N is 12*(4**N) - 1
+            dtype = integer_dtype(16 * (4**refinement_level) - 1)
+        elif indexing_scheme == "zuniq":
+            dtype = np.dtype("int64")
+
     elif new_indexing_scheme in ("nested", "ring"):
-        # The largest possible "nested" or "ring" index at refinement
-        # level N is 12*(4**N) - 1
-        dtype = integer_dtype(12 * (4**refinement_level) - 1)
+        if indexing_scheme == "nuniq":
+            # The largest possible nuniq index at refinement level N
+            # is 16*(4**N) - 1 = 4*(4**N) + 12*(4**N) - 1
+            dtype = integer_dtype(16 * (4**moc_refinement_level) - 1)
+        elif indexing_scheme in ("nested", "ring"):
+            dtype = healpix_index.dtype
+        elif indexing_scheme == "zuniq":
+            dtype = np.dtype("int64")
+
     elif new_indexing_scheme == "zuniq":
-        dtype = np.dtype('int64')
+        dtype = np.dtype("int64")
+
     else:
         raise NotImplementedError(
-            "Can't yet change HEALPix indexing scheme from "
+            "Can't change HEALPix indexing scheme from "
             f"{indexing_scheme!r} to {new_indexing_scheme!r}"
         )  # pragma: no cover
-    
 
     dx = dx.map_blocks(
         cf_healpix_indexing_scheme,
@@ -471,30 +489,34 @@ def _healpix_indexing_scheme(f, hp, new_indexing_scheme):
         new_indexing_scheme=new_indexing_scheme,
         healpix_index_dtype=dtype,
         refinement_level=refinement_level,
+        moc_refinement_level=moc_refinement_level,
     )
     healpix_index.set_data(dx, copy=False)
 
     # If a Dimension Coordinate is now not monotonically ordered,
     # convert it to an Auxiliary Coordinate
-    if healpix_index.construct_type == "dimension_coordinate" and set(
-        (indexing_scheme, new_indexing_scheme)
-    ) != set(("nested", "nuniq")):
+    if (
+        healpix_index.construct_type == "dimension_coordinate"
+        and "ring" in set((indexing_scheme, new_indexing_scheme))
+    ):
         f.dimension_to_auxiliary(hp["coordinate_key"], inplace=True)
 
     # Create new cached elements
-    if old_cache:
+    if old_cache and moc_refinement_level is None:
         new_cache = {}
         for i, value in old_cache.items():
             new_cache[i] = cf_healpix_indexing_scheme(
                 np.array(value),
-                indexing_scheme,
-                new_indexing_scheme,
-                healpix_index_dtype,
-                refinement_level=refinement_level
+                indexing_scheme=indexing_scheme,
+                new_indexing_scheme=new_indexing_scheme,
+                healpix_index_dtype=dtype,
+                refinement_level=refinement_level,
+                moc_refinement_level=moc_refinement_level,
             )
-            
+
         healpix_index.data._set_cached_elements(new_cache)
-        
+
+
 def _healpix_locate(lat, lon, f):
     """Locate HEALPix cells containing latitude-longitude locations.
 
@@ -784,11 +806,11 @@ def healpix_max_refinement_level():
 
     """
     try:
-        from healpix import nside2order
+        import healpix
     except ImportError as e:
         raise ImportError(
             f"{e}. Must install healpix (https://pypi.org/project/healpix) "
             "to find the maximum HEALPix refinement level"
         )
 
-    return nside2order(healpix._chp.NSIDE_MAX)
+    return healpix.nside2order(healpix._chp.NSIDE_MAX)