diff --git a/Changelog.rst b/Changelog.rst index 1bbaf8583d..ebd09a5130 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -1,10 +1,26 @@ Version NEXTVERSION -------------------- - -**2025-12-??** +-------------- +**2025-??-??** + +* Support for HEALPix grids + (https://github.com/NCAS-CMS/cf-python/issues/909) +* New HEALPix methods: `cf.Field.healpix_info`, + `cf.Field.healpix_decrease_refinement_level`, + `cf.Field.healpix_increase_refinement_level`, + `cf.Field.healpix_indexing_scheme`, `cf.Field.healpix_to_ugrid`, + `cf.Domain.create_healpix` + (https://github.com/NCAS-CMS/cf-python/issues/909) +* New method: `cf.Field.create_latlon_coordinates` + (https://github.com/NCAS-CMS/cf-python/issues/909) +* New method: `cf.Data.coarsen` + (https://github.com/NCAS-CMS/cf-python/issues/909) +* New function: `cf.locate` + (https://github.com/NCAS-CMS/cf-python/issues/909) * Reduce the time taken to import `cf` (https://github.com/NCAS-CMS/cf-python/issues/902) +* New optional dependency: ``healpix>=2025.1`` +* Changed dependency: ``cfdm>=1.13.0.0, <1.13.1.0`` ---- @@ -34,8 +50,8 @@ Version 3.18.1 * Allow multiple conditions for the same axis in `cf.Field.subspace` and `cf.Field.indices` (https://github.com/NCAS-CMS/cf-python/issues/881) -* Fix bug in `cf.Field.collapse` that causes a ``ValueError`` to be raised - for missing external cell measures data +* Fix bug in `cf.Field.collapse` that causes a ``ValueError`` to be + raised for missing external cell measures data (https://github.com/NCAS-CMS/cf-python/issues/885) * New dependency: ``distributed>=2025.5.1`` * Changed dependency: ``cfdm>=1.12.3.0, <1.12.4.0`` diff --git a/README.md b/README.md index 8829d96f6a..97743a0257 100644 --- a/README.md +++ b/README.md @@ -89,6 +89,7 @@ of its array manipulation and can: * create new field constructs in memory, * write and append field and domain constructs to netCDF datasets on disk, * read, create, and manipulate UGRID mesh topologies, +* read, write, and manipulate HEALPix grids, * read, write, and create coordinates defined by geometry cells, * read netCDF and CDL datasets containing hierarchical groups, * inspect field constructs, @@ -105,11 +106,12 @@ of its array manipulation and can: * manipulate field construct data by arithmetical and trigonometrical operations, * perform weighted statistical collapses on field constructs, - including those with geometry cells and UGRID mesh topologies, + including those with geometry cells, UGRID mesh topologies, and + HEALPix grids, * perform histogram, percentile and binning operations on field constructs, -* regrid structured grid, mesh and DSG field constructs with - (multi-)linear, nearest neighbour, first- and second-order +* regrid structured grid, UGRID, HEALPix, and DSG field constructs + with (multi-)linear, nearest neighbour, first- and second-order conservative and higher order patch recovery methods, including 3-d regridding, and large-grid support, * apply convolution filters to field constructs, diff --git a/cf/__init__.py b/cf/__init__.py index deaac4a813..b4e6c715c5 100644 --- a/cf/__init__.py +++ b/cf/__init__.py @@ -119,7 +119,6 @@ from .tiepointindex import TiePointIndex from .bounds import Bounds -from .domain import Domain from .datum import Datum from .coordinateconversion import CoordinateConversion @@ -138,6 +137,7 @@ from .cellconnectivity import CellConnectivity from .cellmethod import CellMethod from .cellmeasure import CellMeasure +from .domain import Domain from .domainancillary import DomainAncillary from .domainaxis import DomainAxis from .domaintopology import DomainTopology diff --git a/cf/constants.py b/cf/constants.py index 685bbddb83..0c5f9650ed 100644 --- a/cf/constants.py +++ b/cf/constants.py @@ -493,6 +493,37 @@ }, } +# -------------------------------------------------------------------- +# CF cell methods +# -------------------------------------------------------------------- +cell_methods = set( + ( + "point", + "sum", + "maximum", + "maximum_absolute_value", + "median", + "mid_range", + "minimum", + "minimum_absolute_value", + "mean", + "mean_absolute_value", + "mean_of_upper_decile", + "mode", + "range", + "root_mean_square", + "standard_deviation", + "sum_of_squares", + "variance", + ) +) + + +# -------------------------------------------------------------------- +# CF HEALPix indexing schemes supported by cf +# -------------------------------------------------------------------- +healpix_indexing_schemes = ("nested", "ring", "nuniq") + # -------------------------------------------------------------------- # Logging level setup diff --git a/cf/data/dask_utils.py b/cf/data/dask_utils.py index 20751482df..45d66d1746 100644 --- a/cf/data/dask_utils.py +++ b/cf/data/dask_utils.py @@ -464,3 +464,558 @@ def cf_filled(a, fill_value=None): """ a = cfdm_to_memory(a) return np.ma.filled(a, fill_value=fill_value) + + +def cf_healpix_bounds( + a, + indexing_scheme, + refinement_level=None, + latitude=False, + longitude=False, + pole_longitude=None, +): + """Calculate HEALPix cell bounds. + + Latitude or longitude locations of the cell vertices are derived + from HEALPix indices. Each cell has four bounds which are returned + in an anticlockwise direction, as seen from above, starting with + the northern-most vertex. + + 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 + + M. Reinecke and E. Hivon: Efficient data structures for masks on + 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549 + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.Field.create_latlon_coordinates` + + :Parameters: + + a: `numpy.ndarray` + The array of HEALPix indices. + + indexing_scheme: `str` + The HEALPix indexing scheme. One of ``'nested'``, + ``'ring'``, or ``'nuniq'``. + + refinement_level: `int` or `None`, optional + The refinement level of the grid within the HEALPix + hierarchy, starting at 0 for the base tessellation with 12 + cells. Must be an `int` for *indexing_scheme* ``'nested'`` + or ``'ring'``, but is ignored for *indexing_scheme* + ``'nuniq'``, in which case *refinement_level* may be + `None`. + + latitude: `bool`, optional + If True then return the bounds' latitudes. + + longitude: `bool`, optional + If True then return the bounds' longitudes. + + pole_longitude: `None` or number, optional + Define the longitudes of vertices that lie exactly on the + north or south pole. If `None` (the default) then the + longitude of such a vertex on the north (south) pole will + be the same as the longitude of the south (north) vertex + of the same cell. If set to a number, then the longitudes + of all vertices on the north or south pole will be given + the value *pole_longitude*. Ignored if *longitude* is + False. + + :Returns: + + `numpy.ndarray` + A 2-d array containing the HEALPix cell bounds. + + **Examples** + + >>> cf.data.dask_utils.cf_healpix_bounds( + ... np.array([0, 1, 2, 3]), 'nested', 1, latitude=True + ) + array([[41.8103149 , 19.47122063, 0. , 19.47122063], + [66.44353569, 41.8103149 , 19.47122063, 41.8103149 ], + [66.44353569, 41.8103149 , 19.47122063, 41.8103149 ], + [90. , 66.44353569, 41.8103149 , 66.44353569]]) + >>> cf.data.dask_utils.cf_healpix_bounds( + ... np.array([0, 1, 2, 3]), 'nested', 1, longitude=True + ) + array([[45. , 22.5, 45. , 67.5], + [90. , 45. , 67.5, 90. ], + [ 0. , 0. , 22.5, 45. ], + [45. , 0. , 45. , 90. ]]) + >>> cf.data.dask_utils.cf_healpix_bounds( + ... np.array([0, 1, 2, 3]), 'nested', 1, longitude=True, + ... pole_longitude=3.14159 + ) + array([[45. , 22.5 , 45. , 67.5 ], + [90. , 45. , 67.5 , 90. ], + [ 0. , 0. , 22.5 , 45. ], + [ 3.14159, 0. , 45. , 90. ]]) + + """ + try: + import healpix + except ImportError as e: + raise ImportError( + f"{e}. Must install healpix (https://pypi.org/project/healpix) " + "for the calculation of latitude/longitude coordinate bounds " + "of a HEALPix grid" + ) + + a = cfdm_to_memory(a) + + # Keep an eye on https://github.com/ntessore/healpix/issues/66 + if a.ndim > 1: + raise ValueError( + "Can only calculate HEALPix cell bounds when the " + f"healpix_index array has one dimension. Got shape: {a.shape}" + ) + + if latitude: + pos = 1 + elif longitude: + pos = 0 + + # Define the function that's going to calculate the bounds from + # the HEALPix indices + match indexing_scheme: + case "ring": + bounds_func = healpix._chp.ring2ang_uv + case "nested" | "nuniq": + bounds_func = healpix._chp.nest2ang_uv + case _: + raise ValueError( + "Can't calculate HEALPix cell bounds: Unknown " + f"'indexing_scheme' in cf_healpix_bounds: " + f"{indexing_scheme!r}" + ) + + # Define the cell vertices in an anticlockwise direction, as seen + # from above, starting with the northern-most vertex. Each vertex + # is defined in the form needed by `bounds_func`. + north = (1, 1) + west = (0, 1) + south = (0, 0) + east = (1, 0) + vertices = (north, west, south, east) + + # Initialise the output bounds array + b = np.empty((a.size, 4), dtype="float64") + + match indexing_scheme: + case "nuniq": + # Create bounds for 'nuniq' indices + orders, a = healpix.uniq2pix(a, nest=True) + for order in np.unique(orders): + nside = healpix.order2nside(order) + indices = np.where(orders == order)[0] + for j, (u, v) in enumerate(vertices): + thetaphi = bounds_func(nside, a[indices], u, v) + b[indices, j] = healpix.lonlat_from_thetaphi(*thetaphi)[ + pos + ] + + del orders, indices + + case "nested" | "ring": + # Create bounds for 'nested' or 'ring' indices + nside = healpix.order2nside(refinement_level) + for j, (u, v) in enumerate(vertices): + thetaphi = bounds_func(nside, a, u, v) + b[:, j] = healpix.lonlat_from_thetaphi(*thetaphi)[pos] + + del thetaphi, a + + if longitude: + # Ensure that longitude bounds are less than 360 + where_ge_360 = np.where(b >= 360) + if where_ge_360[0].size: + b[where_ge_360] -= 360.0 + del where_ge_360 + + # Vertices on the north or south pole come out with a + # longitude of NaN, so replace these with sensible values: + # Either the constant 'pole_longitude', or else the longitude + # of the cell vertex that is opposite the vertex on the pole. + north = 0 + south = 2 + for pole, replacement in ((north, south), (south, north)): + indices = np.argwhere(np.isnan(b[:, pole])).flatten() + if not indices.size: + continue + + if pole_longitude is None: + longitude = b[indices, replacement] + else: + longitude = pole_longitude + + b[indices, pole] = longitude + + return b + + +def cf_healpix_coordinates( + a, indexing_scheme, refinement_level=None, latitude=False, longitude=False +): + """Calculate HEALPix cell centre coordinates. + + 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 + + M. Reinecke and E. Hivon: Efficient data structures for masks on + 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549 + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.Field.create_latlon_coordinates` + + :Parameters: + + a: `numpy.ndarray` + The array of HEALPix indices. + + indexing_scheme: `str` + The HEALPix indexing scheme. One of ``'nested'``, + ``'ring'``, or ``'nuniq'``. + + refinement_level: `int` or `None`, optional + The refinement level of the grid within the HEALPix + hierarchy, starting at 0 for the base tessellation with 12 + cells. Must be an `int` for *indexing_scheme* ``'nested'`` + or ``'ring'``, but is ignored for *indexing_scheme* + ``'nuniq'``, in which case *refinement_level* may be + `None`. + + latitude: `bool`, optional + If True then return the coordinate latitudes. + + longitude: `bool`, optional + If True then return the coordinate longitudes. + + :Returns: + + `numpy.ndarray` + A 1-d array containing the HEALPix cell coordinates. + + **Examples** + + >>> cf.data.dask_utils.cf_healpix_coordinates( + ... np.array([0, 1, 2, 3]), 'nested', 1, latitude=True + ) + array([19.47122063, 41.8103149 , 41.8103149 , 66.44353569]) + >>> cf.data.dask_utils.cf_healpix_coordinates( + ... np.array([0, 1, 2, 3]), 'nested', 1, longitude=True + ) + array([45. , 67.5, 22.5, 45. ]) + + """ + try: + import healpix + except ImportError as e: + raise ImportError( + f"{e}. Must install healpix (https://pypi.org/project/healpix) " + "for the calculation of latitude/longitude coordinates of a " + "HEALPix grid" + ) + + a = cfdm_to_memory(a) + + if a.ndim > 1: + raise ValueError( + "Can only calculate HEALPix cell coordinates when the " + f"healpix_index array has one dimension. Got shape: {a.shape}" + ) + + if latitude == longitude: + raise ValueError( + "Can't calculate HEALPix cell coordinates: " + f"latitude={latitude!r} and longitude={longitude!r}" + ) + + if latitude: + pos = 1 + elif longitude: + pos = 0 + + match indexing_scheme: + case "nuniq": + # Create coordinates for 'nuniq' indices + c = np.empty(a.shape, dtype="float64") + + nest = True + orders, a = healpix.uniq2pix(a, nest=nest) + for order in np.unique(orders): + nside = healpix.order2nside(order) + indices = np.where(orders == order)[0] + c[indices] = healpix.pix2ang( + nside=nside, ipix=a[indices], nest=nest, lonlat=True + )[pos] + + case "nested" | "ring": + # Create coordinates for 'nested' or 'ring' indices + nest = indexing_scheme == "nested" + nside = healpix.order2nside(refinement_level) + c = healpix.pix2ang( + nside=nside, + ipix=a, + nest=nest, + lonlat=True, + )[pos] + + case _: + raise ValueError( + "Can't calculate HEALPix cell coordinates: Unknown " + f"'indexing_scheme': {indexing_scheme!r}" + ) + + return c + + +def cf_healpix_indexing_scheme( + a, + indexing_scheme, + new_indexing_scheme, + healpix_index_dtype, + refinement_level=None, +): + """Change the indexing scheme of HEALPix indices. + + 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 + + M. Reinecke and E. Hivon: Efficient data structures for masks on + 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549 + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.Field.healpix_indexing_scheme` + + :Parameters: + + a: `numpy.ndarray` + The array of HEALPix indices. + + indexing_scheme: `str` + The original HEALPix indexing scheme. One of ``'nested'``, + ``'ring'``, or ``'nuniq'``. + + new_indexing_scheme: `str` + The new HEALPix indexing scheme to change to. One of + ``'nested'``, ``'ring'``, or ``'nuniq'``. + + healpix_index_dtype: `str` or `numpy.dtype` + Typecode or data-type to which the new indices will be + cast. This should be the smallest data type needed for + storing the largest possible index value at the refinement + level. + + refinement_level: `int` or `None`, optional + The refinement level of the grid within the HEALPix + hierarchy, starting at 0 for the base tessellation with 12 + cells. Must be an `int` for *indexing_scheme* ``'nested'`` + or ``'ring'``, but is ignored for *indexing_scheme* + ``'nuniq'`` (in which case *refinement_level* may be + `None`). + + :Returns: + + `numpy.ndarray` + An array containing the new HEALPix indices. + + **Examples** + + >>> cf.data.dask_utils.cf_healpix_indexing_scheme( + ... [0, 1, 2, 3], 'nested', 'ring', 1 + ... ) + array([13, 5, 4, 0]) + >>> cf.data.dask_utils.cf_healpix_indexing_scheme( + ... [0, 1, 2, 3], 'nested', 'nuniq', 1 + ) + array([16, 17, 18, 19]) + >>> cf.data.dask_utils.cf_healpix_indexing_scheme( + ... [16, 17, 18, 19], 'nuniq', 'nested', None + ) + array([0, 1, 2, 3]) + + """ + if new_indexing_scheme == indexing_scheme: + # Null operation + return a + + from ..constants import healpix_indexing_schemes + + if new_indexing_scheme not in healpix_indexing_schemes: + raise ValueError( + "Can't change HEALPix indexing scheme: Unknown " + f"'new_indexing_scheme' in cf_healpix_indexing_scheme: " + f"{new_indexing_scheme!r}" + ) + + try: + import healpix + except ImportError as e: + raise ImportError( + f"{e}. Must install healpix (https://pypi.org/project/healpix) " + "for changing the HEALPix indexing scheme" + ) + + a = cfdm_to_memory(a) + + match indexing_scheme: + case "nested": + match new_indexing_scheme: + case "ring": + nside = healpix.order2nside(refinement_level) + a = healpix.nest2ring(nside, a) + + case "nuniq": + a = healpix.pix2uniq(refinement_level, a, nest=True) + + case "ring": + match new_indexing_scheme: + case "nested": + nside = healpix.order2nside(refinement_level) + a = healpix.ring2nest(nside, a) + + case "nuniq": + a = healpix.pix2uniq(refinement_level, a, nest=False) + + case "nuniq": + match new_indexing_scheme: + case "nested" | "ring": + nest = new_indexing_scheme == "nested" + order, a = healpix.uniq2pix(a, nest=nest) + + order = np.unique(order) + if order.size > 1: + raise ValueError( + "Can't change HEALPix indexing scheme from " + f"'nuniq' to {new_indexing_scheme!r}: " + "HEALPix indices span multiple refinement levels " + f"(at least levels {order.tolist()})" + ) + + case _: + raise ValueError( + "Can't change HEALPix indexing scheme: Unknown " + f"'indexing_scheme' in cf_healpix_indexing_scheme: " + f"{indexing_scheme!r}" + ) + + # Cast the new indices to the given data type + a = a.astype(healpix_index_dtype, copy=False) + return a + + +def cf_healpix_weights(a, indexing_scheme, measure=False, radius=None): + """Calculate HEALPix cell area weights. + + 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 + + M. Reinecke and E. Hivon: Efficient data structures for masks on + 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549 + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.weights.Weights.healpix_area` + + :Parameters: + + a: `numpy.ndarray` + The array of HEALPix 'nuniq' indices. + + indexing_scheme: `str` + The HEALPix indexing scheme. Must be ``'nuniq'``. + + measure: `bool`, optional + If True then create weights that are actual cell areas, in + units of the square of the radius units. + + radius: number, optional + The radius of the sphere. Must be set if *measure * is + True, otherwise ignored. + + :Returns: + + `numpy.ndarray` + An array containing the HEALPix cell weights. + + **Examples** + + >>> cf.data.dask_utils.cf_healpix_weights( + ... [76, 77, 78, 79, 20, 21], 'nuniq' + ) + array([0.0625, 0.0625, 0.0625, 0.0625, 0.25 , 0.25 ]) + >>> cf.data.dask_utils.cf_healpix_weights( + ... [76, 77, 78, 79, 20, 21], 'nuniq', + ... measure=True, radius=6371000 + ) + array([2.65658579e+12, 2.65658579e+12, 2.65658579e+12, 2.65658579e+12, + 1.06263432e+13, 1.06263432e+13]) + + """ + if indexing_scheme != "nuniq": + raise ValueError( + "cf_healpix_weights: Can only calulate weights for the " + "'nuniq' indexing scheme" + ) + + try: + import healpix + except ImportError as e: + raise ImportError( + f"{e}. Must install healpix (https://pypi.org/project/healpix) " + "for the calculation of cell area weights of a HEALPix grid" + ) + + a = cfdm_to_memory(a) + + if a.ndim != 1: + raise ValueError( + "Can only calculate HEALPix cell area weights when the " + f"healpix_index array has one dimension. Got shape: {a.shape}" + ) + + # Each cell at refinement level N has weight x/(4**N), where ... + if measure: + # Cell weights equal cell areas. Surface area of a sphere is + # 4*pi*(r**2), number of HEALPix cells at refinement level N + # is 12*(4**N) => Area of each cell is (pi*(r**2)/3)/(4**N) + x = np.pi * (radius**2) / 3.0 + else: + # Normalised weights + x = 1.0 + + orders = healpix.uniq2pix(a, nest=True)[0] + orders, index, inverse = np.unique( + orders, return_index=True, return_inverse=True + ) + + # Initialise the output weights array + w = np.empty(a.shape, dtype="float64") + + # For each refinement level N, put the weights (= x/4**N) into 'w' + # at the correct locations. + for N, i in zip(orders, index): + w = np.where(inverse == inverse[i], x / (4**N), w) + + return w diff --git a/cf/data/data.py b/cf/data/data.py index e906087c93..a89a41a7fd 100644 --- a/cf/data/data.py +++ b/cf/data/data.py @@ -1349,6 +1349,97 @@ def ceil(self, inplace=False, i=False): d._set_dask(dx) return d + @_inplace_enabled(default=False) + def coarsen( + self, + reduction, + axes, + trim_excess=False, + inplace=False, + ): + """Coarsen the data. + + Coarsen the data by applying the *reduction* function to + combine the elements within fixed-size neighborhoods. + + .. versionadded:: NEXTVERSION + + :Parameters: + + reduction: function + The function with which to coarsen the data. + + axes: `dict` + Define how to coarsening neighbourhood for each + axis. A dictionary key is an integer axis position, + with correponding value giving the integer size of the + coarsening neighbourhood for that axis. Unspecified + axes are not coarsened, which is equivalent to + providing a coarsening neighbourhood of ``1``. + + *Example:* + Coarsen the axis in position 1 by combining every 4 + elements: ``{1: 4}`` + + *Example:* + Coarsen the first axis by combining every 3 + elements, and the last axis by combining every 4 + elements: ``{0: 3, -1: 4}`` + + trim_excess: `bool`, optional + If True then omit a partially-full neighbourhood at + the end of a coarsened axis. If False (the default) + then an exception is raised if there are any + partially-filled neighbourhoods. + + {{inplace: `bool`, optional}} + + :Returns: + + `Data` or `None` + The coarsened data, or `None` if the operation was + in-place. + + **Examples** + + >>> import numpy as np + >>> d = cf.Data(np.arange(24).reshape((4, 6))) + >>> print(d.array) + [[ 0 1 2 3 4 5] + [ 6 7 8 9 10 11] + [12 13 14 15 16 17] + [18 19 20 21 22 23]] + >>> e = d.coarsen(np.min, {0: 2, 1: 3}) + >>> print(e.array) + [[ 0 3] + [12 15]] + >>> e = d.coarsen(np.max, {-1: 5}, trim_excess=True) + >>> print(e.array) + [[ 4] + [10] + [16] + [22]] + >>> e = d.coarsen(np.max, {-1: 5}, trim_excess=False) + ValueError: Coarsening factors {1: 5} do not align with array shape (4, 6). + + """ + import dask.array as da + + d = _inplace_enabled_define_and_cleanup(self) + + # Parse axes, making sure that all axes are non-negative. + ndim = self.ndim + for k in axes: + if k < -ndim or k > ndim: + raise ValueError("axis {k} is out of bounds for {ndim}-d data") + + axes = {(k + ndim if k < 0 else k): v for k, v in axes.items()} + + dx = d.to_dask_array() + dx = da.coarsen(reduction, dx, axes, trim_excess=trim_excess) + d._set_dask(dx) + return d + @_inplace_enabled(default=False) def convolution_filter( self, @@ -6258,8 +6349,8 @@ def reshape(self, *shape, merge_chunks=True, limit=None, inplace=False): """Change the shape of the data without changing its values. It assumes that the array is stored in row-major order, and - only allows for reshapings that collapse or merge dimensions - like ``(1, 2, 3, 4) -> (1, 6, 4)`` or ``(64,) -> (4, 4, 4)``. + only allows for reshapings that collapse or merge dimensions, + e.g. ``(1, 2, 3, 4) -> (1, 6, 4)`` or ``(64,) -> (4, 4, 4)``. :Parameters: diff --git a/cf/dimensioncoordinate.py b/cf/dimensioncoordinate.py index 644709bbda..6b0752a00a 100644 --- a/cf/dimensioncoordinate.py +++ b/cf/dimensioncoordinate.py @@ -440,8 +440,7 @@ def direction(self): @classmethod def create_regular(cls, args, units=None, standard_name=None, bounds=True): - """ - Create a new `DimensionCoordinate` with the given range and cellsize. + """Create a new `DimensionCoordinate` with the given range and cellsize. .. versionadded:: 3.15.0 @@ -455,9 +454,9 @@ def create_regular(cls, args, units=None, standard_name=None, bounds=True): bounds: `bool`, optional If True (the default) then the given range represents - the bounds, and the coordinate points will be the midpoints of - the bounds. If False, the range represents the coordinate points - directly. + the bounds, and the coordinate points will be the + midpoints of the bounds. If False, the range + represents the coordinate points directly. units: str or `Units`, optional The units of the new `DimensionCoordinate` object. @@ -492,31 +491,34 @@ def create_regular(cls, args, units=None, standard_name=None, bounds=True): f"Expected a sequence of three numbers, got {args}." ) - range = (args[0], args[1]) + r = (args[0], args[1]) cellsize = args[2] - range_diff = range[1] - range[0] + range_diff = r[1] - r[0] if cellsize > 0 and range_diff <= 0: raise ValueError( - f"Range ({range[0], range[1]}) must be increasing for a " + f"Range {r} must be increasing for a " f"positive cellsize ({cellsize})" ) elif cellsize < 0 and range_diff >= 0: raise ValueError( - f"Range ({range[0], range[1]}) must be decreasing for a " + f"Range {r} must be decreasing for a " f"negative cellsize ({cellsize})" ) + elif cellsize > abs(range_diff): + raise ValueError( + f"cellsize ({cellsize}) can not exceed the range {r}" + ) if standard_name is not None and not isinstance(standard_name, str): raise ValueError("standard_name must be either None or a string.") + start = r[0] + end = r[1] if bounds: cellsize2 = cellsize / 2 - start = range[0] + cellsize2 - end = range[1] - cellsize2 - else: - start = range[0] - end = range[1] + start += cellsize2 + end -= cellsize2 points = np.arange(start, end + cellsize, cellsize) @@ -527,8 +529,11 @@ def create_regular(cls, args, units=None, standard_name=None, bounds=True): ) if bounds: - b = coordinate.create_bounds() - coordinate.set_bounds(b, copy=False) + if points.size > 1: + coordinate.create_bounds(inplace=True) + else: + b = Bounds(data=Data([r])) + coordinate.set_bounds(b, copy=False) return coordinate diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index dfd65e28d1..bf9b9e035f 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -92,6 +92,16 @@ elements will be automatically reduced if including the full amount defined by the halo would extend the subspace beyond the axis limits.""", + # HEALPix references + "{{HEALPix references}}": """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 + + M. Reinecke and E. Hivon: Efficient data structures for masks + on 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549""", # ---------------------------------------------------------------- # Method description substitutions (3 levels of indentation) # ---------------------------------------------------------------- @@ -569,6 +579,30 @@ If True then do not perform the regridding, rather 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 + coarser refinement level cell with consecutive + 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 optmised 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. + + A Multi-Order Coverage (MOC) has pixels with different + refinement levels stored in the same array. The + "nuniq" scheme for an MOC has a unique index for each + cell at each refinement level, such that within each + refinement level a nested-type scheme is employed. In + the "nuniq" scheme, pixels at different refinement + levels inside a single coarser refinement level cell + can have widely different indices. For each refinement + level *n*, the "nuniq" indices are in the range + :math:`4^{(n+1)}, ..., 4^{(n+2)}-1`.""", # dst_grid_partitions "{{dst_grid_partitions: `int` or `str`, optional}}": """dst_grid_partitions: `int` or `str`, optional Calculating the weights matrix for grids with a very large diff --git a/cf/domain.py b/cf/domain.py index 98ca769cdb..7a73c6928a 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -1,15 +1,19 @@ from math import prod +from numbers import Integral import cfdm import numpy as np from . import mixin from .auxiliarycoordinate import AuxiliaryCoordinate +from .bounds import Bounds from .constructs import Constructs +from .coordinatereference import CoordinateReference from .data import Data from .decorators import _inplace_enabled, _inplace_enabled_define_and_cleanup from .dimensioncoordinate import DimensionCoordinate from .domainaxis import DomainAxis +from .domaintopology import DomainTopology from .functions import ( _DEPRECATION_ERROR_ARG, _DEPRECATION_ERROR_METHOD, @@ -78,11 +82,14 @@ class Domain(mixin.FieldDomain, mixin.Properties, cfdm.Domain): def __new__(cls, *args, **kwargs): """Creates a new Domain instance.""" instance = super().__new__(cls) + instance._Bounds = Bounds instance._Constructs = Constructs + instance._CoordinateReference = CoordinateReference instance._Data = Data instance._DomainAxis = DomainAxis instance._DimensionCoordinate = DimensionCoordinate instance._AuxiliaryCoordinate = AuxiliaryCoordinate + instance._DomainTopology = DomainTopology return instance @property @@ -151,12 +158,12 @@ def close(self): @classmethod def create_regular(cls, x_args, y_args, bounds=True): - """ - Create a new domain with the regular longitudes and latitudes. + """Create a new domain with the regular longitudes and latitudes. .. versionadded:: 3.15.1 - .. seealso:: `cf.DimensionCoordinate.create_regular` + .. seealso:: `cf.DimensionCoordinate.create_regular`, + `cf.Domain.create_healpix` :Parameters: @@ -180,7 +187,6 @@ def create_regular(cls, x_args, y_args, bounds=True): **Examples** - >>> import cf >>> domain = cf.Domain.create_regular((-180, 180, 1), (-90, 90, 1)) >>> domain.dump() -------- @@ -257,6 +263,208 @@ def create_regular(cls, x_args, y_args, bounds=True): return domain + @classmethod + def create_healpix( + cls, refinement_level, indexing_scheme="nested", radius=None + ): + r"""Create a new global HEALPix grid. + + The HEALPix axis of the new Domain is ordered so that the + HEALPix indices are monotonically increasing. + + **Performance** + + Very high refinement levels may require the setting of a very + large Dask chunksize, to prevent a possible run-time failure + 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 chunksize of 128 MiB, but + with a chunksize of 1 pebibyte only 24576 chunks are + required:: + + >>> cf.chunksize() + >>> 134217728 + >>> d = cf.Domain.create_healpix(10) + >>> assert d.coord('healpix_index').data.npartitions == 1 + >>> d = cf.Domain.create_healpix(15) + >>> assert d.coord('healpix_index').data.npartitions == 768 + >>> with cf.chunksize('1 PiB'): + ... d = cf.Domain.create_healpix(29) + ... assert d.coord('healpix_index').data.npartitions == 24576 + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.Domain.create_regular`, + `cf.Domain.create_latlon_coordinates` + + :Parameters: + + refinement_level: `int` + The refinement level of the grid within the HEALPix + hierarchy, starting at 0 for the base tessellation + with 12 cells. The number of cells in the global + HEALPix grid for refinement level *n* is + :math:`12\times 4^n`. + + indexing_scheme: `str` + The HEALPix indexing scheme. One of ``'nested'`` (the + default), ``'ring'``, or ``'nuniq'``. + + {{HEALPix indexing schemes}} + + radius: optional + Specify the radius of the latitude-longitude plane + defined in spherical polar coordinates. May be set to + any numeric scalar object, including `numpy` and + `Data` objects. The units of the radius are assumed to + be metres, unless specified by a `Data` object. If the + special value ``'earth'`` is given then the radius + taken as 6371229 metres. If `None` (the default) then + no radius is set. + + *Example:* + Equivalent ways to set a radius of 6371229 metres: + ``6371229``, ``numpy.array(6371229)``, + ``cf.Data(6371229)``, ``cf.Data(6371229, 'm')``, + ``cf.Data(6371.229, 'km')``, ``'earth'``. + + :Returns: + + `Domain` + The newly created HEALPix domain. + + **Examples** + + .. code-block:: python + + >>> d = cf.Domain.create_healpix(4) + >>> d.dump() + -------- + Domain: + -------- + Domain Axis: healpix_index(3072) + + Dimension coordinate: healpix_index + standard_name = 'healpix_index' + Data(healpix_index(3072)) = [0, ..., 3071] + + Coordinate reference: grid_mapping_name:healpix + Coordinate conversion:grid_mapping_name = healpix + Coordinate conversion:indexing_scheme = nested + Coordinate conversion:refinement_level = 4 + Dimension Coordinate: healpix_index + + .. code-block:: python + + >>> d = cf.Domain.create_healpix(4, "nuniq", radius=6371000) + >>> d.dump() + -------- + Domain: + -------- + Domain Axis: healpix_index(3072) + + Dimension coordinate: healpix_index + standard_name = 'healpix_index' + Data(healpix_index(3072)) = [1024, ..., 4095] + + Coordinate reference: grid_mapping_name:healpix + Coordinate conversion:grid_mapping_name = healpix + Coordinate conversion:indexing_scheme = nuniq + Datum:earth_radius = 6371000.0 + Dimension Coordinate: healpix_index + + .. code-block:: python + + >>> d.create_latlon_coordinates(inplace=True) + >>> print(d) + Dimension coords: healpix_index(3072) = [1024, ..., 4095] + Auxiliary coords: latitude(healpix_index(3072)) = [2.388015463268772, ..., -2.388015463268786] degrees_north + : longitude(healpix_index(3072)) = [45.0, ..., 315.0] degrees_east + Coord references: grid_mapping_name:healpix + + """ + import dask.array as da + + from .constants import healpix_indexing_schemes + from .healpix import healpix_max_refinement_level + + if ( + not isinstance(refinement_level, Integral) + or refinement_level < 0 + or refinement_level > healpix_max_refinement_level() + ): + raise ValueError( + "Can't create HEALPix Domain: 'refinement_level' must be a " + "non-negative integer less than or equal to " + f"{healpix_max_refinement_level()}. Got {refinement_level!r}" + ) + + if indexing_scheme not in healpix_indexing_schemes: + raise ValueError( + "Can't create HEALPix Domain: 'indexing_scheme' must be one " + f"of {healpix_indexing_schemes!r}. Got {indexing_scheme!r}" + ) + + nuniq = indexing_scheme == "nuniq" + + domain = Domain() + ncells = 12 * (4**refinement_level) + + # domain_axis: ncdim%cell + d = domain._DomainAxis(ncells) + d.nc_set_dimension("cell") + axis = domain.set_construct(d, copy=False) + + # dimension_coordinate: healpix_index + c = domain._DimensionCoordinate() + c.set_properties({"standard_name": "healpix_index"}) + c.nc_set_variable("healpix_index") + + # Create the healpix_index data + if nuniq: + start = 4 ** (refinement_level + 1) + else: + start = 0 + + stop = start + ncells + dtype = cfdm.integer_dtype(stop - 1) + data = Data(da.arange(start, stop, dtype=dtype)) + + # Set cached data elements + 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) + + # coordinate_reference: grid_mapping_name:healpix + cr = domain._CoordinateReference() + cr.nc_set_variable("healpix") + cr.set_coordinates({key}) + + if radius is not None: + radius = domain.radius(default=radius) + cr.datum.set_parameter("earth_radius", radius.datum()) + + cr.coordinate_conversion.set_parameters( + { + "grid_mapping_name": "healpix", + "indexing_scheme": indexing_scheme, + } + ) + if not nuniq: + cr.coordinate_conversion.set_parameter( + "refinement_level", refinement_level + ) + + domain.set_construct(cr) + + return domain + @_inplace_enabled(default=False) def flip(self, axes=None, inplace=False): """Flip (reverse the direction of) domain axes. @@ -546,7 +754,10 @@ def indices(self, *config, **kwargs): Metadata constructs are selected conditions are specified on their data. Indices for subspacing are then automatically - inferred from where the conditions are met. + inferred from where the conditions are met. If a condition is + a callable function then if is automateically replaced with + the result of calling that function with the Domain as its + only argument. Metadata constructs and the conditions on their data are defined by keyword parameters. @@ -901,7 +1112,10 @@ def subspace(self, *config, **kwargs): Subspacing by metadata selects metadata constructs and specifies conditions on their data. Indices for subspacing are - then automatically inferred from where the conditions are met. + then automatically inferred from where the conditions are + met. If a condition is a callable function then if is + automateically replaced with the result of calling that + function with the Domain as its only argument. Metadata constructs and the conditions on their data are defined by keyword parameters. diff --git a/cf/field.py b/cf/field.py index 53b889259d..2eb1177733 100644 --- a/cf/field.py +++ b/cf/field.py @@ -1,6 +1,7 @@ import logging from dataclasses import dataclass from functools import reduce +from numbers import Integral from operator import mul as operator_mul import cfdm @@ -18,6 +19,7 @@ Domain, DomainAncillary, DomainAxis, + DomainTopology, FieldList, Flags, Index, @@ -178,8 +180,6 @@ # -------------------------------------------------------------------- _collapse_ddof_methods = set(("sd", "var")) -_earth_radius = 6371229.0 - _relational_methods = ( "__eq__", "__ne__", @@ -281,6 +281,7 @@ def __new__(cls, *args, **kwargs): instance._Domain = Domain instance._DomainAncillary = DomainAncillary instance._DomainAxis = DomainAxis + instance._DomainTopology = DomainTopology instance._Quantization = Quantization instance._RaggedContiguousArray = RaggedContiguousArray instance._RaggedIndexedArray = RaggedIndexedArray @@ -2431,8 +2432,8 @@ def cell_area( Specify the radius used for calculating the areas of cells defined in spherical polar coordinates. The radius is that which would be returned by this call of - the field construct's `~cf.Field.radius` method: - ``f.radius(radius)``. See `cf.Field.radius` for + the field construct's `radius` method: + ``f.radius(default=radius)``. See `radius` for details. By default *radius* is ``'earth'`` which means that if @@ -2608,112 +2609,6 @@ def get_domain(self): return domain - def radius(self, default=None): - """Return the radius of a latitude-longitude plane defined in - spherical polar coordinates. - - The radius is taken from the datums of any coordinate - reference constructs, but if and only if this is not possible - then a default value may be used instead. - - .. versionadded:: 3.0.2 - - .. seealso:: `bin`, `cell_area`, `collapse`, `weights` - - :Parameters: - - default: optional - The radius is taken from the datums of any coordinate - reference constructs, but if and only if this is not - possible then the value set by the *default* parameter - is used. May be set to any numeric scalar object, - including `numpy` and `Data` objects. The units of the - radius are assumed to be metres, unless specified by a - `Data` object. If the special value ``'earth'`` is - given then the default radius taken as 6371229 - metres. If *default* is `None` an exception will be - raised if no unique datum can be found in the - coordinate reference constructs. - - *Parameter example:* - Five equivalent ways to set a default radius of - 6371200 metres: ``6371200``, - ``numpy.array(6371200)``, ``cf.Data(6371200)``, - ``cf.Data(6371200, 'm')``, ``cf.Data(6371.2, - 'km')``. - - :Returns: - - `Data` - The radius of the sphere, in units of metres. - - **Examples** - - >>> f.radius() - - - >>> g.radius() - ValueError: No radius found in coordinate reference constructs and no default provided - >>> g.radius('earth') - - >>> g.radius(1234) - - - """ - radii = [] - for cr in self.coordinate_references(todict=True).values(): - r = cr.datum.get_parameter("earth_radius", None) - if r is not None: - r = Data.asdata(r) - if not r.Units: - r.override_units("m", inplace=True) - - if r.size != 1: - radii.append(r) - continue - - got = False - for _ in radii: - if r == _: - got = True - break - - if not got: - radii.append(r) - - if len(radii) > 1: - raise ValueError( - "Multiple radii found from coordinate reference " - f"constructs: {radii!r}" - ) - - if not radii: - if default is None: - raise ValueError( - "No radius found from coordinate reference constructs " - "and no default provided" - ) - - if isinstance(default, str): - if default != "earth": - raise ValueError( - "The default radius must be numeric, 'earth', " - "or None" - ) - - return Data(_earth_radius, "m") - - r = Data.asdata(default).squeeze() - else: - r = Data.asdata(radii[0]).squeeze() - - if r.size != 1: - raise ValueError(f"Multiple radii: {r!r}") - - r.Units = Units("m") - r.dtype = float - return r - def laplacian_xy( self, x_wrap=None, one_sided_at_boundary=False, radius=None ): @@ -3461,6 +3356,18 @@ def weights( ): # Found linear weights from dimension coordinates pass + elif Weights.healpix_area( + self, + da_key, + comp, + weights_axes, + measure=measure, + radius=radius, + methods=methods, + auto=True, + ): + # Found area weights from HEALPix cells + pass weights_axes = [] for key in comp: @@ -3628,6 +3535,18 @@ def weights( ): # Found area weights from UGRID/geometry cells area_weights = True + elif Weights.healpix_area( + self, + None, + comp, + weights_axes, + measure=measure, + radius=radius, + methods=methods, + auto=True, + ): + # Found area weights from HEALPix cells + area_weights = True if not area_weights: raise ValueError( @@ -3668,6 +3587,18 @@ def weights( ): # Found linear weights from line geometries pass + elif Weights.healpix_area( + self, + da_key, + comp, + weights_axes, + measure=measure, + radius=radius, + methods=methods, + auto=True, + ): + # Found area weights from HEALPix cells + pass else: Weights.linear( self, @@ -4871,6 +4802,696 @@ def bin( return out + def healpix_decrease_refinement_level( + self, + refinement_level, + method, + reduction=None, + conform=True, + check_healpix_index=True, + ): + r"""Decrease the refinement level of a HEALPix grid. + + Decreasing the HEALPix refinement level reduces the resolution + of the HEALPix grid by combining, using the given *method*, + the Field data from the original cells that lie inside each + larger cell at the new lower refinement level. + + A new cell method is added, when appropriate, to describe the + reduction. + + The operation requires that a new larger cell at the lower + refinement level either + + * contains no original cells, in which case that larger cell + is not included in the output, + + or + + * is completely covered by original cells. + + It is not allowed for a larger cell to be only partially + covered by original cells. For instance, if the original + refinement level is 10 and the new refinement level is 8, then + each output cell will be the combination of :math:`16\equiv + 4^2\equiv 4^{(10-8)}` original cells, and if a larger cell + contains at least one but fewer than 16 original cells then an + exception is raised (assuming that *check_healpix_index* is + True). + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + .. seealso:: `healpix_increase_refinement_level`, + `healpix_info`, `healpix_indexing_scheme`, + `healpix_to_ugrid`, `collapse` + + :Parameters: + + refinement_level: `int` or `None` + Specify the new lower refinement level as a + non-negative integer less than or equal to the current + refinement level, or if `None` then the refinement + level is not changed. + + method: `str` + The method used to calculate the values in the new + larger cells, from the data on the original + cells. Must be one of these CF standardised cell + method names: ``'maximum'``, + ``'maximum_absolute_value'``, ``'mean'``, + ``'mean_absolute_value'``, ``'mean_of_upper_decile'``, + ``'median'``, ``'mid_range'``, ``'minimum'``, + ``'minimum_absolute_value'``, ``'mode'``, ``'range'``, + ``'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). + + 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. + + For some methods there are default *reduction* + functions, which are only used when *reduction* is + `None` (the default): + + ======================== =================== + *method* Default *reduction* + ======================== =================== + ``'maximum'`` `np.max` + ``'mean'`` `np.mean` + ``'median'`` `np.median` + ``'minimum'`` `np.min` + ``'standard_deviation'`` `np.std` + ``'sum'`` `np.sum` + ``'variance'`` `np.var` + ======================== =================== + + Note that these methods may be also calculated by + another function provided by the *reduction* + parameter. + + conform: `bool`, optional + 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. + + .. note:: Setting to False will speed up the operation + when the HEALPix indexing scheme is already + nested and the HEALPix axis is already + sorted monotonically. + + check_healpix_index: `bool`, optional + If True (the default) then the following conditions + will be checked before the creation of the new Field + (and after the HEALPix grid has been conformed, when + *conform* is True): + + 1. 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. + + If False then these checks are not carried out. + + .. warning:: Only set to False, which will speed up + the operation, when it is known in + 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 + worse**, the operation could complete and + return incorrect data values. + + :Returns: + + `Field` + A new Field with the new HEALPix grid. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(48) = [0, ..., 47] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> f.healpix_info()['refinement_level'] + 1 + + Set the refinement level to 0, showing that every 4 cells + (i.e. the number of cells at the original refinement level + that lie in one cell of the lower refinement leve1) in the + orginal field correspond to one cell at the lower level: + + >>> g = f.healpix_decrease_refinement_level(0, 'maximum') + >>> print(g) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(12)) K + Cell methods : time(2): mean area: mean area: maximum + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(12) = [0, ..., 11] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> g.healpix_info()['refinement_level'] + 0 + >>> print(f[0, :4].array) + [[291.5 293.5 285.3 286.3]] + >>> print(g[0, 0].array) + [[293.5]] + + Set the refinement level to 0 using the ``'range'`` *method*, + which requires a *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) + ... + >>> g = f.healpix_decrease_refinement_level(0, 'range', range_func) + >>> print(g) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(12)) K + Cell methods : time(2): mean area: mean area: range + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(12) = [0, ..., 11] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> print(g[0, 0].array) + [[8.2]] + + """ + from .constants import cell_methods + + # "point" is not a reduction method + cell_methods = cell_methods.copy() + cell_methods.remove("point") + + f = self.copy() + + # Parse 'method' + if method not in cell_methods: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "'method' is not one of the standardised CF cell methods " + f"{sorted(cell_methods)}. Got: {method!r}" + ) + + # Parse 'reduction' + if reduction is None: + # Use a default reduction function for selected methods + match method: + case "maximum": + reduction = np.max + case "minimum": + reduction = np.min + case "mean": + reduction = np.mean + case "standard_deviation": + reduction = np.std + case "variance": + reduction = np.var + case "sum": + reduction = np.sum + case "median": + reduction = np.median + case _: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "Must provide a 'reduction' function for " + f"method={method!r}" + ) + + elif not callable(reduction): + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + f"'reduction' must be a callable. Got: {reduction!r} of " + f"type {type(reduction)}" + ) + + # Get the HEALPix info + hp = f.healpix_info() + + indexing_scheme = hp.get("indexing_scheme") + if indexing_scheme is None: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "indexing_scheme has not been set in the healpix grid " + "mapping coordinate reference" + ) + + old_refinement_level = hp.get("refinement_level") + if old_refinement_level is None: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "refinement_level has not been set in the healpix grid " + "mapping coordinate reference" + ) + + # Decreasing the refinement level requires the nested indexing + # scheme with ordered HEALPix indices + if conform: + try: + f = f.healpix_indexing_scheme("nested", sort=True) + except ValueError as e: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: {e}" + ) + + # 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}: " + "indexing_scheme in the healpix grid mapping coordinate " + f"reference is {indexing_scheme!r}, and not 'nested'. " + "Consider setting conform=True" + ) + + # Parse 'level' + if refinement_level is None: + # No change in refinement level + return f + + if ( + not isinstance(refinement_level, Integral) + or refinement_level < 0 + or refinement_level > old_refinement_level + ): + raise ValueError( + 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}" + ) + + # 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 larger cell at the new lower + # refinement level + ncells = 4 ** (old_refinement_level - refinement_level) + + # Get the healpix_index coordinates + healpix_index = hp.get("healpix_index") + if healpix_index is None: + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "There are no healpix_index coordinates" + ) + + if check_healpix_index: + d = healpix_index.data + if not (d.diff() > 0).all(): + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "Nested healpix_index coordinates are not strictly " + "monotonically increasing. Consider setting conform=True" + ) + + if (d[::ncells] % ncells).any() or ( + d[ncells - 1 :: ncells] - d[::ncells] != ncells - 1 + ).any(): + raise ValueError( + f"Can't decrease HEALPix refinement level of {f!r}: " + "At least one cell at the new lower refinement level " + f"({refinement_level}) contains fewer than {ncells} " + "cells at the original refinement level " + f"({old_refinement_level})" + ) + + # Get the HEALPix axis + axis = hp["domain_axis_key"] + iaxis = f.get_data_axes().index(axis) + + # Whether or not to create lat/lon coordinates for the new + # refinement level. Only do so if the original grid has + # lat/lon coordinates. + create_latlon = bool( + f.coordinates( + "latitude", + "longitude", + filter_by_axis=(axis,), + axis_mode="exact", + todict=True, + ) + ) + + # ------------------------------------------------------------ + # Chenge the refinement level of the Field's data + # ------------------------------------------------------------ + + # Note: Using 'Data.coarsen' works because a) the HEALPix + # indexing scheme is "nested", and b) each new coarser + # cell contains the maximum possible number + # (i.e. 'ncells') of original cells. + f.data.coarsen( + reduction, axes={iaxis: ncells}, trim_excess=False, inplace=True + ) + + # Re-size the HEALPix axis + domain_axis = f.domain_axis(axis) + domain_axis.set_size(f.shape[iaxis]) + + # ------------------------------------------------------------ + # Change the refinement level of the Field's metadata + # ------------------------------------------------------------ + + # Coarsen the 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), so + # we use np.mean for the reduction function. + 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) + domain_ancillary.data.coarsen( + np.mean, axes={iaxis: ncells}, trim_excess=False, inplace=True + ) + + # Coarsen the cell measure constructs that span the HEALPix + # axis. Cell measure data are extensive (i.e. depend on the + # size of the cell), so we use np.sum for the reduction + # function. + 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) + cell_measure.data.coarsen( + np.sum, axes={iaxis: ncells}, trim_excess=False, inplace=True + ) + + # Remove all other metadata constructs that span the HEALPix + # axis, including the original healpix_index coordinates and + # any lat/lon coordinates. + for key in ( + f.constructs.filter_by_axis(axis, axis_mode="and") + .filter_by_type("cell_measure", "domain_ancillary") + .inverse_filter(1) + .todict() + ): + f.del_construct(key) + + # ------------------------------------------------------------ + # Change the refinement level of the healpix_index coordinates + # ------------------------------------------------------------ + healpix_index = healpix_index[::ncells] // ncells + hp_key = f.set_construct(healpix_index, axes=axis, copy=False) + + # ------------------------------------------------------------ + # Update the healpix Coordinate Reference + # ------------------------------------------------------------ + cr = hp.get("grid_mapping_name:healpix") + cr.coordinate_conversion.set_parameter( + "refinement_level", refinement_level + ) + cr.set_coordinate(hp_key) + + # ------------------------------------------------------------ + # Update the cell methods + # ------------------------------------------------------------ + f._update_cell_methods( + method=method, + input_axes=("area",), + domain_axes=f.domain_axes(axis, todict=True), + ) + + # ------------------------------------------------------------ + # Create lat/lon coordinates for the new refinement level + # ------------------------------------------------------------ + if create_latlon: + f.create_latlon_coordinates(two_d=False, inplace=True) + + return f + + def healpix_increase_refinement_level(self, refinement_level, quantity): + """Increase the refinement level of a HEALPix grid. + + 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 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 + (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 + areas (by dividing them by the number of new cells per + original cell). + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + .. seealso:: `healpix_decrease_refinement_level`, + `healpix_info`, `healpix_indexing_scheme`, + `healpix_to_ugrid` + + :Parameters: + + refinement_level: `int` or `None` + Specify the new higher refinement level as an integer + greater than or equal to the current refinement level, + or if `None` then the refinement level is not changed. + + quantity: `str` + Whether the Field data represent intensive or + extensive quantities, specified with ``'intensive'`` + and ``'extensive'`` respectively. + + :Returns: + + `Field` + A new Field with the new HEALPix grid. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(48) = [0, ..., 47] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> f.healpix_info()['refinement_level'] + 1 + >>> g = f.healpix_increase_refinement_level(3, 'intensive') + >>> g + + >>> g.healpix_info()['refinement_level'] + 3 + >>> g = f.healpix_increase_refinement_level(10, 'intensive') + >>> g + + + Set the refinement level to 2, showing that, for an intensive + quantity, every 4 cells at the higher refinement level + (i.e. the number of cells at the higher refinement level that + lie in one cell of the original refinement leve1) have the + same value as a single cell at the original refinement level: + + >>> g = f.healpix_increase_refinement_level(2, 'intensive') + >>> g + + >>> g.healpix_info()['refinement_level'] + 2 + >>> print(f[0, :2] .array) + [[291.5 293.5]] + >>> 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 the ``f`` is this example is + not, but we can assume that it is for demonstration purposes), + each output cell has the value of the original cell in which + it lies, divided by the ratio of the cells' areas (4 in this + case): + + >>> g = f.healpix_increase_refinement_level(2, 'extensive') + >>> print(f[0, :2] .array) + [[291.5 293.5]] + >>> print(g[0, :8] .array) + [[72.875 72.875 72.875 72.875 73.375 73.375 73.375 73.375]] + + """ + from .healpix import ( + _healpix_increase_refinement_level, + _healpix_increase_refinement_level_indices, + healpix_max_refinement_level, + ) + + # Increasing the refinement level requires the nested indexing + # scheme + try: + f = self.healpix_indexing_scheme("nested", sort=False) + except ValueError as e: + raise ValueError( + f"Can't increase HEALPix refinement level of {self!r}: {e}" + ) + + # Get the HEALPix info + hp = f.healpix_info() + old_refinement_level = hp["refinement_level"] + + # Parse 'refinement_level' + if refinement_level is None: + # No change in refinement level + return f + + if ( + not isinstance(refinement_level, Integral) + or refinement_level < old_refinement_level + or refinement_level > healpix_max_refinement_level() + ): + raise ValueError( + f"Can't increase HEALPix refinement level of {f!r}: " + "'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}" + ) + + if refinement_level == old_refinement_level: + # No change in refinement level + return f + + # Parse 'quantity' + valid_quantities = ("intensive", "extensive") + if quantity not in valid_quantities: + 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}" + ) + + # Get the number of cells at the higher refinement level which + # are contained in one cell at the original refinement level + ncells = 4 ** (refinement_level - old_refinement_level) + + # Get the HEALPix axis + axis = hp["domain_axis_key"] + + # Whether or not to create lat/lon coordinates for the new + # refinement level. Only do so if the original grid has + # lat/lon coordinates. + create_latlon = bool( + f.coordinates( + "latitude", + "longitude", + filter_by_axis=(axis,), + axis_mode="exact", + todict=True, + ) + ) + + # Find the position of the HEALPix axis in the Field data + try: + iaxis = f.get_data_axes().index(axis) + except ValueError: + # The Field data doesn't span the HEALPix axis, so insert + # it. + f.insert_dimension(axis, -1, inplace=True) + iaxis = f.get_data_axes().index(axis) + + # Re-size the HEALPix axis + domain_axis = f.domain_axis(axis) + domain_axis.set_size(f.shape[iaxis] * ncells) + + # Increase the refinement level 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. they do not depend on the + # size of the cell). + 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) + _healpix_increase_refinement_level( + domain_ancillary, ncells, iaxis, "intensive" + ) + + # Increase the refinement level of cell measure constructs + # that span the HEALPix axis. Cell measure data are extensive + # (i.e. they depend on the size of the cell). + 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) + _healpix_increase_refinement_level( + cell_measure, ncells, iaxis, "extensive" + ) + + # Remove all other metadata constructs that span the HEALPix + # axis (including the original healpix_index coordinate + # construct, and any lat/lon coordinate constructs that span + # the HEALPix axis). + for key in ( + f.constructs.filter_by_axis(axis, axis_mode="and") + .filter_by_type("cell_measure", "domain_ancillary") + .inverse_filter(1) + ): + f.del_construct(key) + + # Create the healpix_index coordinates for the new refinement + # level + healpix_index = hp["healpix_index"] + _healpix_increase_refinement_level_indices( + healpix_index, ncells, refinement_level + ) + hp_key = f.set_construct(healpix_index, axes=axis, copy=False) + + # Update the healpix Coordinate Reference + cr = hp.get("grid_mapping_name:healpix") + cr.coordinate_conversion.set_parameter( + "refinement_level", refinement_level + ) + cr.set_coordinate(hp_key) + + if create_latlon: + # Create 1-d lat/lon coordinates for the new refinement + # level + f.create_latlon_coordinates(two_d=False, inplace=True) + + return f + def histogram(self, digitized): """Return a multi-dimensional histogram of the data. @@ -5302,7 +5923,9 @@ def collapse( .. versionadded:: 1.0 .. seealso:: `bin`, `cell_area`, `convolution_filter`, - `moving_window`, `radius`, `weights` + `moving_window`, + `healpix_decrease_refinement_level`, `radius`, + `weights` :Parameters: @@ -6516,6 +7139,7 @@ def collapse( input_axes = all_axes all_axes = [] + f_latlon = None for axes in input_axes: if axes is None: domain_axes = self.domain_axes( @@ -6524,6 +7148,10 @@ def collapse( all_axes.append(list(domain_axes)) continue + if f_latlon is None: + # Temporarily create any implied lat/lon coordinates + f_latlon = f.create_latlon_coordinates(cache=False) + axes2 = [] for axis in axes: msg = ( @@ -6541,7 +7169,7 @@ def collapse( msg = "Can't find the collapse axis identified by {!r}" for x in iterate_over: - a = self.domain_axis(x, key=True, default=None) + a = f_latlon.domain_axis(x, key=True, default=None) if a is None: raise ValueError(msg.format(x)) @@ -6549,6 +7177,8 @@ def collapse( all_axes.append(axes2) + del f_latlon + if debug: logger.debug( " all_methods, all_axes, all_within, all_over = " @@ -6579,12 +7209,8 @@ def collapse( if debug: logger.debug( - f" axes = {axes}" - ) # pragma: no cover - logger.debug( - f" method = {method}" - ) # pragma: no cover - logger.debug( + f" axes = {axes}\n" + f" method = {method}\n" f" collapse_axes_all_sizes = {collapse_axes_all_sizes}" ) # pragma: no cover @@ -6872,8 +7498,8 @@ def collapse( # axes. Also delete the corresponding domain ancillaries. # # This is because missing domain ancillaries in a - # coordinate refernce are assumed to have the value zero, - # which is most likely inapproriate. + # coordinate reference are assumed to have the value zero, + # which is most likely inappropriate. # -------------------------------------------------------- if remove_vertical_crs: for ref_key, ref in f.coordinate_references( @@ -6901,6 +7527,26 @@ def collapse( if set(ref_axes).intersection(flat(all_axes)): f.del_coordinate_reference(ref_key) + # -------------------------------------------------------- + # Remove a HEALPix Coordinate Reference and Dimension + # Coordinate. (A HEALPix Auxiliary Coordinate gets removed + # later with the other 1-d Auxiliary Coordinates.) + # -------------------------------------------------------- + healpix_axis = f.domain_axis( + "healpix_index", key=True, default=None + ) + domain_axis = collapse_axes.get(healpix_axis) + if domain_axis is not None and domain_axis.get_size() > 1: + from .healpix import del_healpix_coordinate_reference + + del_healpix_coordinate_reference(f) + + key = f.dimension_coordinate( + "healpix_index", key=True, default=None + ) + if key is not None: + f.del_construct(key) + # --------------------------------------------------------- # Update dimension coordinates, auxiliary coordinates, # cell measures and domain ancillaries @@ -6929,7 +7575,6 @@ def collapse( # REMOVE all 2+ dimensional auxiliary coordinates # which span this axis - # c = auxiliary_coordinates.filter_by_naxes(gt(1), view=True) c = f.auxiliary_coordinates( filter_by_naxes=(gt(1),), filter_by_axis=(axis,), @@ -6944,14 +7589,13 @@ def collapse( f.del_construct(key) - # REMOVE all 1 dimensional auxiliary coordinates which - # span this axis and have different values in their - # data array and bounds. + # REMOVE all 1-d auxiliary coordinates which span this + # axis and have different values in their data array + # and bounds. # - # KEEP, after changing their data arrays, all - # one-dimensional auxiliary coordinates which span - # this axis and have the same values in their data - # array and bounds. + # KEEP, after changing their data arrays, all 1-d + # auxiliary coordinates which span this axis and have + # the same values in their data array and bounds. c = f.auxiliary_coordinates( filter_by_axis=(axis,), axis_mode="exact", todict=True ) @@ -8550,7 +9194,10 @@ def indices(self, *config, **kwargs): Metadata constructs are selected by conditions specified on their data. Indices for subspacing are then automatically - inferred from where the conditions are met. + inferred from where the conditions are met. If a condition is + a callable function then if is automateically replaced with + the result of calling that function with the Field as its only + argument. The returned tuple of indices may be used to created a subspace by indexing the original field construct with them. @@ -12735,10 +13382,13 @@ def subspace(self): **Subspacing by metadata** - Subspacing by metadata, signified by the use of round brackets, - selects metadata constructs and specifies conditions on their - data. Indices for subspacing are then automatically inferred from - where the conditions are met. + Subspacing by metadata, signified by the use of round + brackets, selects metadata constructs and specifies conditions + on their data. Indices for subspacing are then automatically + inferred from where the conditions are met. If a condition is + a callable function then if is automateically replaced with + the result of calling that function with the Field as its only + argument. Metadata constructs and the conditions on their data are defined by keyword parameters. @@ -12767,6 +13417,12 @@ def subspace(self): acting along orthogonal dimensions, some missing data may still need to be inserted into the field construct's data. + * If a condition is a given as a callable function, then it is + replaced with the output of it being called with the field + as its ony argument. For instance, + ``f.subspace(X=cf.locate(30, 180))`` is equivalent to + ``f.subspace(X=cf.locate(30, 180)(f))``. + **Subspacing by index** Subspacing by indexing, signified by the use of square brackets, @@ -13013,7 +13669,22 @@ def regrids( Data defined on UGRID face or node cells may be regridded to any other latitude-longitude grid, including other UGRID - meshes and DSG feature types. + meshes and DSG feature types. Note that for conservative + regridding, a cell edge is assumed to be the great circle arc + that connects its two vertices. + + **HEALPix grids** + + Data on HEALPix grids may be regridded to any other + latitude-longitude grid, including other HEALPix grids, UGRID + meshes and DSG feature types. This is done internally by + converting HEALPix grids to UGRID meshes and carrying out a + UGRID regridding. This means that for conservative regridding, + a cell edge is assumed to be the great circle arc that + connects its two vertices, which is certainly not true for + HEALPix cells. However, the errors are likely to be small, + particularly at higher resolutions. See also + `healpix_decrease_refinement_level`. **DSG feature types** @@ -13037,7 +13708,7 @@ def regrids( .. versionadded:: 1.0.4 - .. seealso:: `regridc` + .. seealso:: `regridc`, `healpix_decrease_refinement_level` :Parameters: diff --git a/cf/functions.py b/cf/functions.py index f63b037305..768a3c45a4 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -3340,6 +3340,136 @@ def unique_constructs(constructs, ignore_properties=None, copy=True): ) +def locate(lat, lon, f=None): + """Locate cells containing latitude-longitude locations. + + The cells must be defined by a discrete axis that has 1-d latitude + and longitude coordinate, or for which it is possible to create + 1-d latitude and longitude coordinates from other metadata (as is + the case for a HEALPix axis). At present, only a HEALPix axis is + supported. + + Returns the indices of the discrete axis for the cells containing + the latitude-longitude locations. + + If a single latitude is given then it is paired with each + longitude, and if a single longitude is given then it is paired + with each latitude. If multiple latitudes and multiple longitudes + are provided then they are paired element-wise. + + If a cell contains more than one of the latitude-longitude + locations, then that cell's index appears only once in the output. + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.contains`, `cf.Field.subspace`, + `cf.Field.indices` + + :Parameters: + + lat: (sequence of) number + The latitude(s), in degrees north, for which to find the + cell indices. Must be in the range [-90, 90]. + + lon: (sequence of) number + The longitude(s), in degrees east, for which to find the + cell indices. + + f: `Field` or `Domain` or `None`, optional + The Field or Domain containing the UGRID or HEALPix + grid. + + If `None` (the default) then a callable function is + returned, which when called with A Field as its argument + returns the indices for that field, i.e ``cf.contains(lat, + lon, f)`` is equivalent to ``cf.contains(lat, + lon)(f)``. This returned function may be used as a + condition in a Field's `subspace` and `indices` methods, + since, for instance, ``f.subspace(X=cf.locate(0, 45))`` is + equivalent to ``f.subspace(X=cf.locate(0, 45, f))``. + + `numpy.ndarray` or function + Indices for the discrete axis that contain the + latitude-longitude locations, or if *f* is `None`, a + function that will return those indices. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : height(1) = [1.5] m + Auxiliary coords: healpix_index(healpix_index(48)) = [0, ..., 47] 1 + Coord references: grid_mapping_name:healpix + >>> cf.locate(20, 90, f) + array([23]) + >>> cf.locate(-70, 90, f) + array([36]) + >>> cf.locate([-70, 20], 90, f) + array([23, 36]) + >>> cf.locate([-70, 20], [90, 280], f) + array([31, 36]) + >>> cf.locate([-70, 20], [90, 280])(f) + array([31, 36]) + >>> cf.locate(20, [280, 280.001], f) + array([31]) + >>> func = cf.locate(20, [280, 280.001]) + >>> func(f) + array([31]) + + """ + + def _locate(lat, lon, f): + healpix = f.coordinate( + "healpix_index", filter_by_naxes=(1,), default=None + ) + if healpix: + # HEALPix + from .healpix import _healpix_locate + + return _healpix_locate(lat, lon, f) + + raise ValueError( + f"Can't find cell locations for {f!r}: Can only find locations " + "for HEALPix cells (at present)" + ) + + ugrid = f.domain_topologies(todict=True) + if ugrid: + # UGRID - not coded up, yet. + pass + + geometry = any( + aux.get_geometry(False) + for aux in f.auxiliary_coordinates( + filter_by_naxes=(1,), todict=True + ).values() + ) + if geometry: + # Geometries - not coded up, yet. + pass + + raise ValueError( + f"Can't find cell locations for {f!r}: Can only find locations " + "for UGRID, HEALPix, or geometry cells" + ) + + if np.abs(lat).max() > 90: + raise ValueError( + f"Can't find cell locations for {f!r}: Latitudes must be in " + f"the range [-90, 90]. Got: {lat}" + ) + + if f is None: + return partial(_locate, lat, lon) + + return _locate(lat, lon, f) + + def _DEPRECATION_ERROR(message="", version="3.0.0", removed_at="4.0.0"): if removed_at: removed_at = f" and will be removed at version {removed_at}" diff --git a/cf/healpix.py b/cf/healpix.py new file mode 100644 index 0000000000..963d5a76be --- /dev/null +++ b/cf/healpix.py @@ -0,0 +1,792 @@ +"""HEALPix functionality.""" + +import logging + +import numpy as np +from cfdm import is_log_level_info + +logger = logging.getLogger(__name__) + + +def _healpix_create_latlon_coordinates(f, pole_longitude, cache=True): + """Create HEALPix latitude and longitude coordinates and bounds. + + When it is not possible to create latitude and longitude + coordinates, the reason why will be reported if the log level is + at ``2``/``'INFO'`` or higher. + + 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.create_latlon_coordinates`, + `cf.Field.healpix_to_ugrid` + + :Parameters: + + f: `Field` or `Domain` + The Field or Domain containing the HEALPix grid, which + will be updated in-place. + + pole_longitude: `None` or number + The longitude of HEALPix coordinate bounds that lie + exactly on the north (south) pole. If `None` then the + longitude of such a vertex will be the same as the south + (north) vertex of the same cell. If set to a number, then + the longitudes of such vertices will all be given that + value. + + cache: `bool`, optional + If True (the default) then cache in memory the first and + last of any newly-created coordinates and bounds. This + may slightly slow down the coordinate creation process, + but may greatly speed up, and reduce the memory + requirement of, a future inspection of the coordinates and + bounds. Even when *cache* is True, new cached coordinate + values can only be created if the existing healpix_index + coordinates themselves have cached first and last values. + + :Returns: + + (`str`, `str`) or (`None`, `None`) + The keys of the new latitude and longitude coordinate + constructs, in that order, or two `None`s if the + coordinates could not be created. + + """ + import dask.array as da + + from .constants import healpix_indexing_schemes + from .data.dask_utils import cf_healpix_bounds, cf_healpix_coordinates + + hp = f.healpix_info() + + indexing_scheme = hp.get("indexing_scheme") + if indexing_scheme not in healpix_indexing_schemes: + if is_log_level_info(logger): + logger.info( + "Can't create 1-d latitude and longitude coordinates for " + f"{f!r}: indexing_scheme in the healpix grid mapping " + "coordinate reference must be one of " + f"{healpix_indexing_schemes!r}. Got {indexing_scheme!r}" + ) # pragma: no cover + + return (None, None) + + refinement_level = hp.get("refinement_level") + if refinement_level is None and indexing_scheme != "nuniq": + if is_log_level_info(logger): + logger.info( + "Can't create 1-d latitude and longitude coordinates for " + f"{f!r}: refinement_level has not been set in the healpix " + "grid mapping coordinate reference" + ) # pragma: no cover + + return (None, None) + + healpix_index = hp.get("healpix_index") + if healpix_index is None: + if is_log_level_info(logger): + logger.info( + "Can't create 1-d latitude and longitude coordinates for " + f"{f!r}: Missing healpix_index coordinates" + ) # pragma: no cover + + return (None, None) + + # Create latitude and longitude cached elements + if cache: + cache = healpix_index.data._get_cached_elements() + cached_lon_coords = {} + cached_lat_coords = {} + cached_lon_bounds = {} + cached_lat_bounds = {} + for i, value in cache.items(): + if i not in (0, -1): + continue + + cached_lon_coords[i] = cf_healpix_coordinates( + value, indexing_scheme, refinement_level, longitude=True + ) + cached_lat_coords[i] = cf_healpix_coordinates( + value, indexing_scheme, refinement_level, latitude=True + ) + cached_lon_bounds[i] = cf_healpix_bounds( + value, + indexing_scheme, + refinement_level, + longitude=True, + pole_longitude=pole_longitude, + )[0, i] + cached_lat_bounds[i] = cf_healpix_bounds( + value, indexing_scheme, refinement_level, latitude=True + )[0, i] + + # Get the Dask array of HEALPix indices. + # + # `cf_healpix_coordinates` and `cf_healpix_bounds` have their own + # calls to `cfdm_to_memory`, so we can set _force_to_memory=False. + dx = healpix_index.data.to_dask_array( + _force_mask_hardness=False, _force_to_memory=False + ) + meta = np.array((), dtype="float64") + + # Create latitude coordinates + dy = dx.map_blocks( + cf_healpix_coordinates, + meta=meta, + indexing_scheme=indexing_scheme, + refinement_level=refinement_level, + latitude=True, + ) + data = f._Data(dy, "degrees_north") + if cache: + data._set_cached_elements(cached_lat_coords) + + lat = f._AuxiliaryCoordinate( + data=data, + properties={"standard_name": "latitude"}, + copy=False, + ) + + # Create longitude coordinates + dy = dx.map_blocks( + cf_healpix_coordinates, + meta=meta, + indexing_scheme=indexing_scheme, + refinement_level=refinement_level, + longitude=True, + ) + data = f._Data(dy, "degrees_east") + if cache: + data._set_cached_elements(cached_lon_coords) + + lon = f._AuxiliaryCoordinate( + data=data, + properties={"standard_name": "longitude"}, + copy=False, + ) + + # Create latitude bounds + dy = da.blockwise( + cf_healpix_bounds, + "ij", + dx, + "i", + new_axes={"j": 4}, + meta=meta, + indexing_scheme=indexing_scheme, + refinement_level=refinement_level, + latitude=True, + ) + data = f._Data(dy, "degrees_north") + if cache: + data._set_cached_elements(cached_lat_bounds) + + bounds = f._Bounds(data=data) + lat.set_bounds(bounds, copy=False) + + # Create longitude bounds + dy = da.blockwise( + cf_healpix_bounds, + "ij", + dx, + "i", + new_axes={"j": 4}, + meta=meta, + indexing_scheme=indexing_scheme, + refinement_level=refinement_level, + longitude=True, + pole_longitude=pole_longitude, + ) + data = f._Data(dy, "degrees_east") + if cache: + data._set_cached_elements(cached_lon_bounds) + + bounds = f._Bounds(data=data) + lon.set_bounds(bounds, copy=False) + + # Set the new latitude and longitude coordinates + axis = hp["domain_axis_key"] + lat_key = f.set_construct(lat, axes=axis, copy=False) + lon_key = f.set_construct(lon, axes=axis, copy=False) + + 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 + in-place. + + 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` + + """ + import dask.array as da + from dask.array.core import normalize_chunks + + # Get any cached data values + cached = x.data._get_cached_elements().copy() + + # Get the Dask array (e.g. dx.shape is (12, 19, 48)) + dx = x.data.to_dask_array(_force_mask_hardness=False) + + if quantity == "extensive": + # Divide extensive data by the number of new cells + dx = dx / ncells + if cached: + cached = {i: value / ncells for i, value in cached.items()} + + # Add a new size dimension just after the HEALPix dimension + # (e.g. dx.shape becomes (12, 19, 48, 1)) + new_axis = iaxis + 1 + dx = da.expand_dims(dx, new_axis) + + # Modify the size of the new dimension to be the number of cells + # at the new refinement level which are contained in one cell at + # the original refinement level (e.g. size becomes 16) + shape = list(dx.shape) + shape[new_axis] = ncells + + # Work out what the chunks should be for the new dimension + chunks = list(dx.chunks) + chunks[new_axis] = "auto" + chunks = normalize_chunks(chunks, shape, dtype=dx.dtype) + + # Broadcast the data along the new dimension (e.g. dx.shape + # becomes (12, 19, 48, 16)) + dx = da.broadcast_to(dx, shape, chunks=chunks) + + # Reshape the array so that it has a single, larger HEALPix + # dimension (e.g. dx.shape becomes (12, 19, 768)) + dx = dx.reshape( + shape[:iaxis] + + [shape[iaxis] * shape[new_axis]] + + shape[new_axis + 1 :] + ) + + x.set_data(dx, copy=False) + + # Set cached data elements + if cached: + x.data._set_cached_elements(cached) + + +def _healpix_increase_refinement_level_indices( + healpix_index, ncells, refinement_level +): + """Increase the refinement level 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 + 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: + + healpix_index: `Coordinate` + The HEALPix indices to be changed. They must use the + "nested" indexing scheme. + + ncells: `int` + The number of cells at the new higher refinement level + which are contained in one cell at the original refinement + level. + + refinement_level: `int` + The new higher refinement level. + + :Returns: + + `None` + + """ + import dask.array as da + from cfdm import integer_dtype + from dask.array.core import normalize_chunks + + # Get any cached data values + cached = healpix_index.data._get_cached_elements().copy() + + # Get the Dask array (e.g. dx.shape is (48,)) + dx = healpix_index.data.to_dask_array(_force_mask_hardness=False) + + # Set the data type to allow for the largest possible HEALPix + # index at the new refinement level + dtype = integer_dtype(12 * (4**refinement_level) - 1) + if dx.dtype != dtype: + dx = dx.astype(dtype) + + # Change each original HEALPix index to the smallest new HEALPix + # index that the larger cell contains + dx = dx * ncells + + # Add a new size dimension just after the HEALPix dimension + # (e.g. dx.shape becomes (48, 1)) + new_axis = 1 + dx = da.expand_dims(dx, new_axis) + + # Modify the size of the new dimension to be the number of cells + # at the new refinement level which are contained in one cell at + # the original refinement level (e.g. size becomes 16) + shape = list(dx.shape) + shape[new_axis] = ncells + + # Work out what the chunks should be for the new dimension + chunks = list(dx.chunks) + chunks[new_axis] = "auto" + chunks = normalize_chunks(chunks, shape, dtype=dx.dtype) + + # Broadcast the data along the new dimension (e.g. dx.shape + # becomes (48, 16)) + dx = da.broadcast_to(dx, shape, chunks=chunks) + + # Increment the broadcast values along the new dimension, so that + # dx[i, :] contains all of the nested indices at the higher + # refinement level that correspond to HEALPix index value i at the + # original refinement level. + new_shape = [1] * dx.ndim + new_shape[new_axis] = shape[new_axis] + + dx += da.arange(ncells, chunks=chunks[new_axis]).reshape(new_shape) + + # Reshape the new array to combine the original HEALPix and + # broadcast dimensions into a single new HEALPix dimension + # (e.g. dx.shape becomes (768,)) + dx = dx.reshape(shape[0] * shape[new_axis]) + + healpix_index.set_data(dx, copy=False) + + # Set new cached data elements + data = healpix_index.data + if 0 in cached: + x = np.array(cached[0], dtype=dtype) * ncells + data._set_cached_elements({0: x, 1: x + 1}) + + if -1 in cached: + x = np.array(cached[-1], dtype=dtype) * ncells + (ncells - 1) + data._set_cached_elements({-1: x}) + + +def _healpix_indexing_scheme(f, hp, new_indexing_scheme): + """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 + Fast Analysis of Data Distributed on the Sphere. The Astrophysical + Journal, 2005, 622 (2), pp.759-771. + https://dx.doi.org/10.1086/427976 + + M. Reinecke and E. Hivon: Efficient data structures for masks on + 2D grids. A&A, 580 (2015) + A132. https://doi.org/10.1051/0004-6361/201526549 + + .. versionadded:: NEXTVERSION + + .. seealso:: `cf.Field.healpix_indexing_scheme` + + :Parameters: + + f: `Field` or `Domain` + The Field or Domain that contains the healpix_index + coordinates. *f* will be updated with the new HEALPix + indices in-place. + + hp: `dict` + The HEALPix info dictionary for *f*. + + new_indexing_scheme: `str` + The new indexing scheme. + + :Returns: + + `None` + + """ + from cfdm import integer_dtype + + from .data.dask_utils import cf_healpix_indexing_scheme + + healpix_index = hp["healpix_index"] + indexing_scheme = hp["indexing_scheme"] + refinement_level = hp.get("refinement_level") + + # Change the HEALPix indices + # + # `cf_healpix_indexing_scheme` has its own call to + # `cfdm_to_memory`, so we can set _force_to_memory=False. + dx = healpix_index.data.to_dask_array( + _force_mask_hardness=False, _force_to_memory=False + ) + + # 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) + else: + # The largest possible "nested" or "ring" index at refinement + # level N is 12*(4**N) - 1 + dtype = integer_dtype(12 * (4**refinement_level) - 1) + + dx = dx.map_blocks( + cf_healpix_indexing_scheme, + dtype=dtype, + meta=np.array((), dtype=dtype), + indexing_scheme=indexing_scheme, + new_indexing_scheme=new_indexing_scheme, + healpix_index_dtype=dtype, + refinement_level=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")): + f.dimension_to_auxiliary(hp["coordinate_key"], inplace=True) + + +def _healpix_locate(lat, lon, f): + """Locate HEALPix cells containing latitude-longitude locations. + + Returns the discrete axis indices of the cells containing the + latitude-longitude locations. + + If a single latitude is given then it is paired with each + longitude, and if a single longitude is given then it is paired + with each latitude. If multiple latitudes and multiple longitudes + are provided then they are paired element-wise. + + If a cell contains more than one of the given latitude-longitude + locations then that cell's index appears only once in the output. + + 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.locate` + + :Parameters: + + lat: (sequence of) number + The latitude(s), in degrees north, for which to find the + cell indices. Must be in the range [-90, 90], although + this is not checked. + + lon: (sequence of) number + The longitude(s), in degrees east, for which to find the + cell indices. + + f: `Field` or `Domain` + The Field or Domain containing the HEALPix grid. + + :Returns: + + `numpy.ndarray` + Indices for the HEALPix axis that contain the + latitude-longitude locations. Note that these indices + identify locations along the HEALPix axis, and are not the + HEALPix indices defined by the indexing scheme. + + """ + import dask.array as da + + try: + import healpix + except ImportError as e: + raise ImportError( + f"{e}. Must install healpix (https://pypi.org/project/healpix) " + "to allow the location of HEALPix cells" + ) + + hp = f.healpix_info() + + healpix_index = hp.get("healpix_index") + if healpix_index is None: + raise ValueError( + f"Can't locate HEALPix cells for {f!r}: There are no " + "healpix_index coordinates" + ) + + indexing_scheme = hp.get("indexing_scheme") + match indexing_scheme: + case "nuniq": + # nuniq indexing scheme + index = [] + healpix_index = healpix_index.array + orders = healpix.uniq2pix(healpix_index, nest=True)[0] + orders = np.unique(orders) + for order in orders: + # For this refinement level, find the HEALPix nested + # indices of the cells that contain the lat-lon + # points. + nside = healpix.order2nside(order) + pix = healpix.ang2pix(nside, lon, lat, nest=True, lonlat=True) + # Remove duplicate indices + pix = np.unique(pix) + # Convert back to HEALPix nuniq indices + pix = healpix._chp.nest2uniq(order, pix, pix) + # Find where these HEALPix indices are located in the + # healpix_index coordinates + index.append(da.where(da.isin(healpix_index, pix))[0]) + + index = da.unique(da.concatenate(index, axis=0)) + + case "nested" | "ring": + # nested or ring indexing scheme + refinement_level = hp.get("refinement_level") + if refinement_level is None: + raise ValueError( + f"Can't locate HEALPix cells for {f!r}: refinement_level " + "has not been set in the healpix grid mapping coordinate " + "reference" + ) + + # Find the HEALPix indices of the cells that contain the + # lat-lon points + nest = indexing_scheme == "nested" + nside = healpix.order2nside(refinement_level) + pix = healpix.ang2pix(nside, lon, lat, nest=nest, lonlat=True) + # Remove duplicate indices + pix = np.unique(pix) + # Find where these HEALPix indices are located in the + # healpix_index coordinates + index = da.where(da.isin(healpix_index, pix))[0] + + case None: + raise ValueError( + f"Can't locate HEALPix cells for {f!r}: indexing_scheme has " + "not been set in the healpix grid mapping coordinate " + "reference" + ) + + case _: + from .constants import healpix_indexing_schemes + + raise ValueError( + f"Can't locate HEALPix cells for {f!r}: indexing_scheme in " + "the healpix grid mapping coordinate reference must be one " + f"of {healpix_indexing_schemes!r}. Got {indexing_scheme!r}" + ) + + # Return the cell locations as a numpy array of element indices + return index.compute() + + +def del_healpix_coordinate_reference(f): + """Remove a healpix grid mapping coordinate reference construct. + + If required, a new latitude_longitude grid mapping coordinate + reference will be created in-place to store any generic coordinate + conversion or datum parameters found in the healpix grid mapping + coordinate reference. + + .. versionadded:: NEXTVERSION + + :Parameters: + + f: `Field` or `Domain` + The Field or Domain from which to delete the healpix grid + mapping coordinate reference. + + :Returns: + + `CoordinateReference` or `None` + The removed healpix grid mapping coordinate reference + construct, or `None` if there wasn't one. + + """ + cr_key, cr = f.coordinate_reference( + "grid_mapping_name:healpix", item=True, default=(None, None) + ) + if cr is not None: + f.del_construct(cr_key) + + latlon = f.coordinate_reference( + "grid_mapping_name:latitude_longitude", default=None + ) + if latlon is None: + latlon = cr.copy() + cc = latlon.coordinate_conversion + cc.del_parameter("grid_mapping_name", None) + cc.del_parameter("indexing_scheme", None) + cc.del_parameter("refinement_level", None) + if cc.parameters() or latlon.datum.parameters(): + # The healpix coordinate reference contains generic + # coordinate conversion or datum parameters + latlon.coordinate_conversion.set_parameter( + "grid_mapping_name", "latitude_longitude" + ) + + # Remove healpix_index coordinates from the coordinate + # reference + for key in f.coordinates( + "healpix_index", filter_by_naxes=(1,), todict=True + ): + latlon.del_coordinate(key, None) + + f.set_construct(latlon) + + return cr + + +def healpix_info(f): + """Get information about the HEALPix grid, if there is one. + + .. versionadded:: NEXTVERSION + + :Parameters: + + f: `Field` or `Domain` + The field or domain. + + :Returns: + + `dict` + The information about the HEALPix axis, with some or all + of the following dictionary keys: + + * ``'coordinate_reference_key'``: The construct key of the + healpix coordinate + reference construct. + + * ``'grid_mapping_name:healpix'``: The healpix coordinate + reference construct. + + * ``'indexing_scheme'``: The HEALPix indexing scheme. + + * ``'refinement_level'``: The refinement level of the + HEALPix grid. + + * ``'domain_axis_key'``: The construct key of the HEALPix + domain axis construct. + + * ``'coordinate_key'``: The construct key of the + healpix_index coordinate + construct. + + * ``'healpix_index'``: The healpix_index coordinate + construct. + + The dictionary will be empty if there is no HEALPix axis. + + **Examples** + + >>> f = cf.example_field(12) + >>> healpix_info(f) + {'coordinate_reference_key': 'coordinatereference0', + 'grid_mapping_name:healpix': , + 'indexing_scheme': 'nested', + 'refinement_level': 1, + 'domain_axis_key': 'domainaxis1', + 'coordinate_key': 'auxiliarycoordinate0', + 'healpix_index': } + + >>> f = cf.example_field(0) + >>> healpix_info(f) + {} + + """ + info = {} + + cr_key, cr = f.coordinate_reference( + "grid_mapping_name:healpix", item=True, default=(None, None) + ) + if cr is not None: + info["coordinate_reference_key"] = cr_key + info["grid_mapping_name:healpix"] = cr + parameters = cr.coordinate_conversion.parameters() + for param in ("indexing_scheme", "refinement_level"): + value = parameters.get(param) + if value is not None: + info[param] = value + + hp_key, healpix_index = f.coordinate( + "healpix_index", + filter_by_naxes=(1,), + item=True, + default=(None, None), + ) + if healpix_index is not None: + info["domain_axis_key"] = f.get_data_axes(hp_key)[0] + info["coordinate_key"] = hp_key + info["healpix_index"] = healpix_index + + return info + + +def healpix_max_refinement_level(): + """Return the maximum permitted HEALPix refinement level. + + The maximum refinement level is the highest refinement level for + which all of HEALPix indices from any indexing scheme are + representable as double precision integers. + + 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 + + :Returns: + + `int` + The maximum permitted HEALPix refinement level. + + """ + try: + 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 healpix.nside2order(healpix._chp.NSIDE_MAX) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 60098d76cd..e770a794c7 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -21,6 +21,7 @@ logger = logging.getLogger(__name__) +_earth_radius = 6371229.0 _units_degrees = Units("degrees") @@ -296,6 +297,16 @@ def _indices(self, config, data_axes, ancillary_mask, kwargs): f"non-negative integer. Got {halo!r}" ) + # Create any latitude and longitude coordinates (e.g. as + # implied by a non-latitude_longitude grid mapping coordinate + # reference). This is so that latitude and longitude kwarg + # selections can work, even if the actual latitude and + # longitude coordinates are not currently part of the + # metadata. + # + # Do not do this in-place. + self = self.create_latlon_coordinates(cache=False) + domain_axes = self.domain_axes(todict=True) # Initialise the index for each axis @@ -323,6 +334,18 @@ def _indices(self, config, data_axes, ancillary_mask, kwargs): f"defined by {identity!r}" ) + # If the condition is a callable function, then call it + # with 'self' as the only argument and replace the + # condition with the result. + if callable(value): + try: + value = value(self) + except Exception as e: + raise RuntimeError( + f"Error encountered when calling condition " + f"{identity}={value}: {e}" + ) + if axes in parsed: # The axes are the same as an existing key parsed[axes].append((axes, key, construct, value, identity)) @@ -1611,8 +1634,15 @@ def auxiliary_to_dimension( axis = f.get_data_axes(key) dim = f._DimensionCoordinate(source=aux) - f.set_construct(dim, axes=axis) + dim_key = f.set_construct(dim, axes=axis, copy=False) + + # Update Coordinates listed in Coordinate References + for cr in f.coordinate_references().values(): + if key in cr.coordinates(): + cr.set_coordinate(dim_key) + f.del_construct(key) + return f def del_coordinate_reference( @@ -1621,66 +1651,66 @@ def del_coordinate_reference( """Remove a coordinate reference construct and all of its domain ancillary constructs. - .. versionadded:: 3.0.0 + .. versionadded:: 3.0.0 - .. seealso:: `del_construct` + .. seealso:: `del_construct` - :Parameters: + :Parameters: - identity: optional - Select the coordinate reference construct by one of: + identity: optional + Select the coordinate reference construct by one of: - * The identity of a coordinate reference construct. + * The identity of a coordinate reference construct. - {{construct selection identity}} + {{construct selection identity}} - * The key of a coordinate reference construct + * The key of a coordinate reference construct - * `None`. This is the default, which selects the - coordinate reference construct when there is only - one of them. + * `None`. This is the default, which selects the + coordinate reference construct when there is only + one of them. - *Parameter example:* - ``identity='standard_name:atmosphere_hybrid_height_coordinate'`` + *Parameter example:* + ``identity='standard_name:atmosphere_hybrid_height_coordinate'`` - *Parameter example:* - ``identity='grid_mapping_name:rotated_latitude_longitude'`` + *Parameter example:* + ``identity='grid_mapping_name:rotated_latitude_longitude'`` - *Parameter example:* - ``identity='transverse_mercator'`` + *Parameter example:* + ``identity='transverse_mercator'`` - *Parameter example:* - ``identity='coordinatereference1'`` + *Parameter example:* + ``identity='coordinatereference1'`` - *Parameter example:* - ``identity='key%coordinatereference1'`` + *Parameter example:* + ``identity='key%coordinatereference1'`` - *Parameter example:* - ``identity='ncvar%lat_lon'`` + *Parameter example:* + ``identity='ncvar%lat_lon'`` - *Parameter example:* - ``identity=cf.eq('rotated_pole')'`` + *Parameter example:* + ``identity=cf.eq('rotated_pole')'`` - *Parameter example:* - ``identity=re.compile('^rotated')`` + *Parameter example:* + ``identity=re.compile('^rotated')`` - construct: optional - TODO + construct: optional + TODO - default: optional - Return the value of the *default* parameter if the - construct can not be removed, or does not exist. + default: optional + Return the value of the *default* parameter if the + construct can not be removed, or does not exist. - {{default Exception}} + {{default Exception}} - :Returns: + :Returns: - The removed coordinate reference construct. + The removed coordinate reference construct. - **Examples** + **Examples** - >>> f.del_coordinate_reference('rotated_latitude_longitude') - + >>> f.del_coordinate_reference('rotated_latitude_longitude') + """ if construct is None: @@ -1947,6 +1977,611 @@ def coordinate_reference_domain_axes(self, identity=None): return set(axes) + def healpix_indexing_scheme(self, new_indexing_scheme, sort=False): + r"""Change the indexing scheme of HEALPix indices. + + Note that the Field data values are not changed, nor is the + Field Data array reordered. Only the "healpix_index" + coordinate values are changed, along with the corresponding + "healpix" grid mapping Coordinate Reference. + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + .. seealso:: `healpix_info`, `healpix_to_ugrid` + + :Parameters: + + new_indexing_scheme: `str` or `None` + The new HEALPix indexing scheme. One of ``'nested'``, + ``'ring'``, ``'nuniq'``, or `None`. If `None` then the + indexing scheme is unchanged. + + {{HEALPix indexing schemes}} + + sort: `bool`, optional + If True then re-order the HEALPix axis of the output + so that its HEALPix indices are monotonically + increasing, including when the indexing scheme is + unchanged. If False (the default) then don't do this. + + :Returns: + + `{{class}}` + The {{class}} with the HEALPix indices redefined for + the new scheme. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(48) = [0, ..., 47] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> f.healpix_info()['indexing_scheme'] + 'nested' + >>> print(f.coordinate('healpix_index').array) + [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 + 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47] + + >>> g = f.healpix_indexing_scheme('nuniq') + >>> g.healpix_info()['indexing_scheme'] + 'nuniq' + >>> print(g.coordinate('healpix_index').array) + [16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 + 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63] + + >>> g = f.healpix_indexing_scheme('ring') + >>> g.healpix_info()['indexing_scheme'] + 'ring' + >>> print(g.coordinate('healpix_index').array) + [13 5 4 0 15 7 6 1 17 9 8 2 19 11 10 3 28 20 27 12 30 22 21 14 + 32 24 23 16 34 26 25 18 44 37 36 29 45 39 38 31 46 41 40 33 47 43 42 35] + + >>> g = f.healpix_indexing_scheme('ring', sort=True) + >>> print(g.coordinate('healpix_index').array) + [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 + 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47] + >>> h = g.healpix_indexing_scheme('nested') + >>> print(h.coordinate('healpix_index').array) + [ 3 7 11 15 2 1 6 5 10 9 14 13 19 0 23 4 27 8 31 12 17 22 21 26 + 25 30 29 18 16 35 20 39 24 43 28 47 34 33 38 37 42 41 46 45 32 36 40 44] + >>> h = g.healpix_indexing_scheme(None, sort=True) + >>> print(h.coordinate('healpix_index').array) + [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 + 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47] + + """ + from ..constants import healpix_indexing_schemes + + f = self.copy() + + hp = f.healpix_info() + + if new_indexing_scheme not in healpix_indexing_schemes + (None,): + raise ValueError( + f"Can't change HEALPix index scheme of {f!r}: " + "new_indexing_scheme keyword must be None or one of " + f"{healpix_indexing_schemes!r}. Got {new_indexing_scheme!r}" + ) + + # Get the healpix_index coordinates + healpix_index = hp.get("healpix_index") + if healpix_index is None: + raise ValueError( + f"Can't change HEALPix index scheme of {f!r}: There are no " + "healpix_index coordinates" + ) + + indexing_scheme = hp.get("indexing_scheme") + if indexing_scheme is None: + raise ValueError( + f"Can't change HEALPix indexing scheme of {f!r}: " + "indexing_scheme has not been set in the healpix grid " + "mapping coordinate reference" + ) + + if indexing_scheme not in healpix_indexing_schemes: + raise ValueError( + f"Can't change HEALPix indexing scheme of {f!r}: " + "indexing_scheme in the healpix grid mapping coordinate " + "reference must be one of " + f"{healpix_indexing_schemes!r}. Got {new_indexing_scheme!r}" + ) + + if ( + new_indexing_scheme is not None + and new_indexing_scheme != indexing_scheme + ): + refinement_level = hp.get("refinement_level") + if ( + indexing_scheme in ("nested", "ring") + and refinement_level is None + ): + raise ValueError( + f"Can't change HEALPix indexing scheme of {f!r} from " + f"{indexing_scheme!r} to {new_indexing_scheme!r} when " + "refinement_level has not been set in the healpix grid " + "mapping coordinate reference" + ) + + # Update the Coordinate Reference + cr = hp["grid_mapping_name:healpix"] + cr.coordinate_conversion.set_parameter( + "indexing_scheme", new_indexing_scheme + ) + if new_indexing_scheme == "nuniq": + cr.coordinate_conversion.del_parameter("refinement_level") + elif indexing_scheme == "nuniq": + raise ValueError( + f"Can't change HEALPix indexing scheme of {f!r} from " + f"{indexing_scheme!r} to {new_indexing_scheme!r}" + ) + + # Change the HEALPix indices + from ..healpix import _healpix_indexing_scheme + + _healpix_indexing_scheme(f, hp, new_indexing_scheme) + + if sort: + # Sort the HEALPix axis so that the HEALPix indices are + # monotonically increasing. + d = healpix_index.data + if (d.diff() < 0).any(): + index = d.compute() + f = f.subspace(**{hp["domain_axis_key"]: np.argsort(index)}) + + # Now that the HEALPix indices are ordered, store them in + # a Dimension Coordinate. + if healpix_index.construct_type == "auxiliary_coordinate": + f.auxiliary_to_dimension(hp["coordinate_key"], inplace=True) + + return f + + def healpix_info(self): + """Get information about the HEALPix grid, if there is one. + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + :Returns: + + `dict` + + The information about the HEALPix axis, with some or + all of the following dictionary keys: + + * ``'coordinate_reference_key'``: The construct key of + the healpix coordinate reference construct. + + * ``'grid_mapping_name:healpix'``: The healpix + coordinate reference construct. + + * ``'indexing_scheme'``: The HEALPix indexing scheme. + + * ``'refinement_level'``: The refinement level of the + HEALPix grid. + + * ``'domain_axis_key'``: The construct key of the + HEALPix domain axis construct. + + * ``'coordinate_key'``: The construct key of the + healpix_index coordinate construct. + + * ``'healpix_index'``: The healpix_index coordinate + construct. + + The dictionary will be empty if there is no HEALPix axis. + + **Examples** + + >>> f = cf.example_field(12) + >>> f.healpix_info() + {'coordinate_reference_key': 'coordinatereference0', + 'grid_mapping_name:healpix': , + 'indexing_scheme': 'nested', + 'refinement_level': 1, + 'domain_axis_key': 'domainaxis1', + 'coordinate_key': 'auxiliarycoordinate0', + 'healpix_index': } + + >>> f = cf.example_field(0) + >>> healpix_info(f) + {} + + """ + from ..healpix import healpix_info + + return healpix_info(self) + + @_inplace_enabled(default=False) + def healpix_to_ugrid(self, cache=True, inplace=False): + """Convert a HEALPix domain to a UGRID domain. + + **References** + + {{HEALPix references}} + + .. versionadded:: NEXTVERSION + + .. seealso:: `healpix_info`, `create_latlon_coordinates` + + :Parameters: + + cache: `bool`, optional + + If True (the default) then cache in memory the first + and last of any newly-created UGRID coordinates and + bounds. This may slightly slow down the coordinate + creation process, but may greatly speed up, and reduce + the memory requirement of, a future inspection of the + coordinates and bounds. Even when *cache* is True, new + cached coordinate values can only be created if the + existing healpix_index coordinates themselves have + cached first and last values. + + {{inplace: `bool`, optional}} + + :Returns: + + `{{class}}` or `None` + The `{{class}}` converted to UGRID, or `None` if the + operation was in-place. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(48) = [0, ..., 47] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> print(f.healpix_to_ugrid()) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), ncdim%cell(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : height(1) = [1.5] m + Auxiliary coords: latitude(ncdim%cell(48)) = [19.47122063449069, ..., -19.47122063449069] degrees_north + : longitude(ncdim%cell(48)) = [45.0, ..., 315.0] degrees_east + Coord references: grid_mapping_name:latitude_longitude + Topologies : cell:face(ncdim%cell(48), 4) = [[774, ..., 3267]] + + """ + from ..healpix import del_healpix_coordinate_reference + + hp = self.healpix_info() + + axis = hp.get("domain_axis_key") + if axis is None: + raise ValueError( + f"Can't convert {self!r} from HEALPix to UGRID: There is no " + "HEALPix domain axis" + ) + + f = _inplace_enabled_define_and_cleanup(self) + + # If 1-d lat/lon coordinates do not exist, then derive them + # from the HEALPix indices. Setting the pole_longitude to + # something other than None - it doesn't matter what - ensures + # that the north (south) polar vertex comes out as a single + # node in the domain topology. + f.create_latlon_coordinates( + two_d=False, pole_longitude=0, cache=cache, inplace=True + ) + + # Get the lat/lon coordinates + x_key, x = f.auxiliary_coordinate( + "Y", + filter_by_axis=(axis,), + axis_mode="exact", + item=True, + default=(None, None), + ) + y_key, y = f.auxiliary_coordinate( + "X", + filter_by_axis=(axis,), + axis_mode="exact", + item=True, + default=(None, None), + ) + if x is None: + raise ValueError( + f"Can't convert {f!r} from HEALPix to UGRID: Not able to " + "find (nor create) longitude coordinates" + ) + + if y is None: + raise ValueError( + f"Can't convert {f!r} from HEALPix to UGRID: Not able to " + "find (or create) latitude coordinates" + ) + + bounds_y = y.get_bounds(None) + bounds_x = x.get_bounds(None) + if bounds_y is None: + raise ValueError( + f"Can't convert {f!r} from HEALPix to UGRID: No latitude " + "coordinate bounds" + ) + + if bounds_x is None: + raise ValueError( + f"Can't convert {f!r} from HEALPix to UGRID: No longitude " + "coordinate bounds" + ) + + # Create the UGRID Domain Topology construct, by creating a + # unique integer identifer for each unique node location. + bounds_y = bounds_y.data.to_dask_array(_force_mask_hardness=False) + bounds_x = bounds_x.data.to_dask_array(_force_mask_hardness=False) + + _, y_indices = np.unique(bounds_y, return_inverse=True) + _, x_indices = np.unique(bounds_x, return_inverse=True) + + nodes = y_indices * y_indices.size + x_indices + + domain_topology = f._DomainTopology(data=f._Data(nodes)) + domain_topology.set_cell("face") + domain_topology.set_property( + "long_name", "UGRID domain topology derived from HEALPix" + ) + f.set_construct(domain_topology, axes=axis, copy=False) + + # Remove the HEALPix index coordinates + f.del_construct(hp["coordinate_key"]) + + # Remove the HEALPix coordinate reference + del_healpix_coordinate_reference(f) + + return f + + @_inplace_enabled(default=False) + @_manage_log_level_via_verbosity + def create_latlon_coordinates( + self, + one_d=True, + two_d=True, + pole_longitude=None, + overwrite=False, + cache=True, + inplace=False, + verbose=None, + ): + """Create latitude and longitude coordinates. + + Creates 1-d or 2-d latitude and longitude coordinate + constructs that are implied by the {{class}} coordinate + reference constructs. By default (or if *overwrite* is False), + new coordinates are only created if the {{class}} doesn't + already include any latitude or longitude coordinates. + + When it is not possible to create latitude and longitude + coordinates, the reason why will be reported if the log level + is at ``2``/``'INFO'`` or higher (as set by `cf.log_level` or + the *verbose* parameter). + + .. versionadded:: NEXTVERSION + + .. seealso:: `healpix_to_ugrid` + + :Parameters: + + one_d: `bool`, optional` + If True (the default) then consider creating 1-d + latitude and longitude coordinates. If False then 1-d + coordinates will not be created. + + two_d: `bool`, optional` + If True (the default) then consider creating 2-d + latitude and longitude coordinates. If False then 2-d + coordinates will not be created. + + pole_longitude: `None` or number + Define the longitudes of coordinates or coordinate + bounds that lie exactly on the north or south pole. If + `None` (the default) then the longitudes of such + points are determined by whichever algorithm was used + to create the coordinates, which will likely result in + different points on a pole having different + longitudes. If set to a number, then the longitudes of + all points on the north or south pole will be given + the value *pole_longitude*. + + overwrite: `bool`, optional + If True then remove any existing latitude and + longitude coordinates, prior to attempting to create + new ones. If False (the default) then if any latitude + or longitude coordinates already exist, new ones will + not be created. + + cache: `bool`, optional + If True (the default) then cache in memory the first + and last of any newly-created coordinates and + bounds. This may slightly slow down the coordinate + creation process, but may greatly speed up, and reduce + the memory requirement of, a future inspection of the + coordinates and bounds. Even when *cache* is True, new + cached values can only be created if the existing + source coordinates (from which the newly-created + coordinates are calculated) themselves have cached + first and last values. + + {{inplace: `bool`, optional}} + + {{verbose: `int` or `str` or `None`, optional}} + + :Returns: + + `{{class}}` or `None` + A new {{class}}, with new latitude and longitude + constructs if any could be created. If the operation + was in-place then `None` is returned. + + **Examples** + + >>> f = cf.example_field(12) + >>> print(f) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(healpix_index(48)) = [0, ..., 47] + : height(1) = [1.5] m + Coord references: grid_mapping_name:healpix + >>> g = f.create_latlon_coordinates() + >>> print(g) + Field: air_temperature (ncvar%tas) + ---------------------------------- + Data : air_temperature(time(2), healpix_index(48)) K + Cell methods : time(2): mean area: mean + Dimension coords: time(2) = [2025-06-16 00:00:00, 2025-07-16 12:00:00] proleptic_gregorian + : healpix_index(48) = [0, ..., 47] + : height(1) = [1.5] m + Auxiliary coords: latitude(healpix_index(48)) = [19.47122063449069, ..., -19.47122063449069] degrees_north + : longitude(healpix_index(48)) = [45.0, ..., 315.0] degrees_east + Coord references: grid_mapping_name:healpix + + """ + f = _inplace_enabled_define_and_cleanup(self) + + # ------------------------------------------------------------ + # See if any lat/lon coordinates could be created + # ------------------------------------------------------------ + + # See if there are any existing latitude/longutude coordinates + latlon_coordinates = { + key: c + for key, c in f.coordinates(todict=True).items() + if c.Units.islatitude or c.Units.islongitude + } + if latlon_coordinates: + if overwrite: + # Remove existing latitude/longutude coordinates + # before carrying on + for key in latlon_coordinates: + f.del_construct(key) + else: + if is_log_level_info(logger): + logger.info( + "Can't create latitude and longitude coordinates: " + "overwrite=False and latitude and/or longitude " + "coordinates already exist: " + f"{', '.join(map(repr, latlon_coordinates.values()))}" + ) # pragma: no cover + + return f + + # Store all of the grid mapping Coordinate References in a + # dictionary + coordinate_references = { + cr.identity(None): cr + for cr in f.coordinate_references(todict=True).values() + } + coordinate_references.pop(None, None) + if not coordinate_references: + if is_log_level_info(logger): + logger.info( + "Can't create latitude and longitude coordinates: " + "There are no grid mapping coordinate references" + ) # pragma: no cover + + return f + + coordinate_references = { + identity: cr + for identity, cr in coordinate_references.items() + if identity.startswith("grid_mapping_name:") + } + + # Remove a 'latitude_longitude' grid mapping (if there is one) + # from the dictionary, saving it for later. + latlon_cr = coordinate_references.pop( + "grid_mapping_name:latitude_longitude", None + ) + if not coordinate_references: + if is_log_level_info(logger): + logger.info( + "Can't create latitude and longitude coordinates: There " + "is no non-latitude_longitude grid mapping coordinate " + "reference" + ) # pragma: no cover + + return f + + if len(coordinate_references) > 1: + if is_log_level_info(logger): + logger.info( + "Can't create latitude and longitude coordinates: There " + "is more than one non-latitude_longitude grid mapping " + "coordinate reference: " + f"{', '.join(map(repr, coordinate_references.values()))}" + ) # pragma: no cover + + return f + + # ------------------------------------------------------------ + # Still here? Then we might be able to create some lat/lon + # coordinates. + # ------------------------------------------------------------ + + # Initialize the flag that tells us if any new coordinates + # have been created + coords_created = False + + # Get the unique non-latitude_longitude grid mapping + identity, cr = coordinate_references.popitem() + + if one_d and not coords_created: + # -------------------------------------------------------- + # 1-d lat/lon coordinates + # -------------------------------------------------------- + if identity == "grid_mapping_name:healpix": + # ---------------------------------------------------- + # HEALPix + # ---------------------------------------------------- + from ..healpix import _healpix_create_latlon_coordinates + + lat_key, lon_key = _healpix_create_latlon_coordinates( + f, pole_longitude, cache + ) + coords_created = lat_key is not None + + if two_d and not coords_created: + # -------------------------------------------------------- + # 2-d lat/lon coordinates + # -------------------------------------------------------- + pass # For now ... + + # ------------------------------------------------------------ + # Update the appropriate coordinate reference with any new + # coordinate keys + # ------------------------------------------------------------ + if coords_created: + if latlon_cr is not None: + latlon_cr.set_coordinates((lat_key, lon_key)) + else: + cr.set_coordinates((lat_key, lon_key)) + + return f + def cyclic( self, *identity, iscyclic=True, period=None, config={}, **filter_kwargs ): @@ -2147,8 +2782,15 @@ def dimension_to_auxiliary( axis = f.get_data_axes(key) aux = f._AuxiliaryCoordinate(source=dim) - f.set_construct(aux, axes=axis) + aux_key = f.set_construct(aux, axes=axis, copy=False) + + # Update Coordinates listed in Coordinate References + for cr in f.coordinate_references().values(): + if key in cr.coordinates(): + cr.set_coordinate(aux_key) + f.del_construct(key) + return f @_deprecated_kwarg_check("axes", version="3.0.0", removed_at="4.0.0") @@ -2413,6 +3055,8 @@ def is_discrete_axis(self, *identity, **filter_kwargs): * An axis spanned by the domain topology construct of an unstructured grid. + * A HEALPix axis. + * The axis with geometry cells. .. versionaddedd:: 3.16.3 @@ -2482,6 +3126,15 @@ def is_discrete_axis(self, *identity, **filter_kwargs): ): return True + # HEALPix + if self.coordinates( + "healpix_index", + filter_by_axis=(axis,), + axis_mode="exact", + todict=True, + ): + return True + # Geometries for aux in self.auxiliary_coordinates( filter_by_axis=(axis,), axis_mode="exact", todict=True @@ -2578,6 +3231,112 @@ def _parse_axes(self, axes): return [self.domain_axis(x, key=True) for x in axes] + def radius(self, default=None): + """Return the radius of a latitude-longitude plane defined in + spherical polar coordinates. + + The radius is taken from the datums of any coordinate + reference constructs, but if and only if this is not possible + then a default value may be used instead. + + .. versionadded:: 3.0.2 + + .. seealso:: `bin`, `cell_area`, `collapse` + + :Parameters: + + default: optional + The radius is taken from the datums of any coordinate + reference constructs, but if and only if this is not + possible then the value set by the *default* parameter + is used. May be set to any numeric scalar object, + including `numpy` and `Data` objects. The units of the + radius are assumed to be metres, unless specified by a + `Data` object. If the special value ``'earth'`` is + given then the default radius taken as 6371229 + metres. If *default* is `None` an exception will be + raised if no unique datum can be found in the + coordinate reference constructs. + + *Parameter example:* + Five equivalent ways to set a default radius of + 6371200 metres: ``6371200``, + ``numpy.array(6371200)``, ``cf.Data(6371200)``, + ``cf.Data(6371200, 'm')``, ``cf.Data(6371.2, + 'km')``. + + :Returns: + + `Data` + The radius of the sphere, in units of metres. + + **Examples** + + >>> f.radius() + + + >>> g.radius() + ValueError: No radius found in coordinate reference constructs and no default provided + >>> g.radius('earth') + + >>> g.radius(1234) + + + """ + radii = [] + for cr in self.coordinate_references(todict=True).values(): + r = cr.datum.get_parameter("earth_radius", None) + if r is not None: + r = self._Data.asdata(r) + if not r.Units: + r.override_units("m", inplace=True) + + if r.size != 1: + radii.append(r) + continue + + got = False + for _ in radii: + if r == _: + got = True + break + + if not got: + radii.append(r) + + if len(radii) > 1: + raise ValueError( + "Multiple radii found from coordinate reference " + f"constructs: {radii!r}" + ) + + if not radii: + if default is None: + raise ValueError( + "No radius found from coordinate reference constructs " + "and no default provided" + ) + + if isinstance(default, str): + if default != "earth": + raise ValueError( + "The default radius must be numeric, 'earth', " + "or None" + ) + + return self._Data(_earth_radius, "m") + + r = self._Data.asdata(default).squeeze() + else: + r = self._Data.asdata(radii[0]).squeeze() + + if r.size != 1: + raise ValueError(f"Multiple radii: {r!r}") + + r.Units = Units("m") + r.dtype = float + return r + def replace_construct( self, *identity, new=None, copy=True, **filter_kwargs ): @@ -2992,10 +3751,13 @@ def set_coordinate_reference( for value in coordinate_reference.coordinates(): if value in coordinates: identity = coordinates[value].identity(strict=strict) - ckeys.append(self.coordinate(identity, key=True, default=None)) + ckey = self.coordinate(identity, key=True, default=None) + if ckey is not None: + ckeys.append(ckey) ref.clear_coordinates() - ref.set_coordinates(ckeys) + if ckeys: + ref.set_coordinates(ckeys) coordinate_conversion = coordinate_reference.coordinate_conversion diff --git a/cf/mixin/propertiesdata.py b/cf/mixin/propertiesdata.py index 2ac8e1ab2c..1c98d038cb 100644 --- a/cf/mixin/propertiesdata.py +++ b/cf/mixin/propertiesdata.py @@ -530,6 +530,8 @@ def __pos__(self): def __query_isclose__(self, value, rtol, atol): """Query interface method for an "is close" condition. + .. versionadded:: 3.15.2 + :Parameters: value: @@ -541,8 +543,6 @@ def __query_isclose__(self, value, rtol, atol): atol: number The tolerance on absolute numerical differences. - .. versionadded:: 3.15.2 - """ data = self.get_data(None, _fill_value=None) if data is None: @@ -5173,6 +5173,8 @@ def rechunk( {{balance: `bool`, optional}} + {{inplace: `bool`, optional}} + :Returns: `{{class}}` or `None` diff --git a/cf/mixin/propertiesdatabounds.py b/cf/mixin/propertiesdatabounds.py index 1150449a56..1b5dd008e2 100644 --- a/cf/mixin/propertiesdatabounds.py +++ b/cf/mixin/propertiesdatabounds.py @@ -945,8 +945,8 @@ def lower_bounds(self): return data.copy() raise AttributeError( - "Can't get lower bounds when there are no bounds nor coordinate " - "data" + f"Can't get lower bounds from {self!r} when there are no bounds " + "nor coordinate data" ) @property @@ -3666,6 +3666,8 @@ def rechunk( If True (the default) then rechunk an interior ring array, if one exists. + {{inplace: `bool`, optional}} + :Returns: `{{class}}` or `None` diff --git a/cf/query.py b/cf/query.py index 320bcbc92c..5551c01988 100644 --- a/cf/query.py +++ b/cf/query.py @@ -1918,9 +1918,9 @@ def contains(value, units=None): .. versionadded:: 3.0.0 - .. seealso:: `cf.Query.iscontains`, `cf.cellsize`, `cf.cellge`, - `cf.cellgt`, `cf.cellne`, `cf.cellle`, `cf.celllt`, - `cf.cellwi`, `cf.cellwo` + .. seealso:: `cf.locate`, `cf.cellsize`, `cf.cellge`, `cf.cellgt`, + `cf.cellne`, `cf.cellle`, `cf.celllt`, `cf.cellwi`, + `cf.cellwo`, `cf.Query.iscontains` :Parameters: diff --git a/cf/regrid/regrid.py b/cf/regrid/regrid.py index b67db47f8e..6091f305db 100644 --- a/cf/regrid/regrid.py +++ b/cf/regrid/regrid.py @@ -88,7 +88,9 @@ class Grid: dummy_size_2_dimension: bool = False # Whether or not the grid is a structured grid. is_grid: bool = False - # Whether or not the grid is a UGRID mesh. + # Whether or not the grid is a UGRID mesh (after any + # transformations are applied, such as converting HEALPix to + # UGRID). is_mesh: bool = False # Whether or not the grid is a location stream. is_locstream: bool = False @@ -120,6 +122,9 @@ class Grid: # The integer position in *coords* of a vertical coordinate. If # `None` then there are no vertical coordinates. z_index: Any = None + # The original field/domain before any transformations are applied + # (such as creating lat/lon coordinates, or converting to UGRID). + domain: Any = None def regrid( @@ -515,8 +520,10 @@ def regrid( ) if debug: + from pprint import pformat + logger.debug( - f"Source Grid:\n{src_grid}\n\nDestination Grid:\n{dst_grid}\n" + f"\n{pformat(src_grid)}\n\n{pformat(dst_grid)}\n" ) # pragma: no cover conform_coordinates(src_grid, dst_grid) @@ -526,6 +533,7 @@ def regrid( raise ValueError( f"{method!r} regridding is not available for 1-d regridding" ) + elif method in ("nearest_dtos", "nearest_stod"): if not has_coordinate_arrays(src_grid) and not has_coordinate_arrays( dst_grid @@ -693,7 +701,7 @@ def regrid( start_index=start_index, src_axes=src_axes, dst_axes=dst_axes, - dst=dst, + dst=dst_grid.domain, weights_file=weights_file if from_file else None, src_mesh_location=src_grid.mesh_location, dst_mesh_location=dst_grid.mesh_location, @@ -770,14 +778,16 @@ def regrid( # ---------------------------------------------------------------- # Update the regridded metadata # ---------------------------------------------------------------- - update_non_coordinates(src, dst, src_grid, dst_grid, regrid_operator) + cr_map = update_non_coordinates( + src, dst, src_grid, dst_grid, regrid_operator + ) - update_coordinates(src, dst, src_grid, dst_grid) + update_coordinates(src, dst, src_grid, dst_grid, cr_map) # ---------------------------------------------------------------- # Insert regridded data into the new field # ---------------------------------------------------------------- - update_data(src, regridded_data, src_grid) + update_data(src, regridded_data, src_grid, dst_grid) if coord_sys == "spherical" and dst_grid.is_grid: # Set the cyclicity of the longitude axis of the new field @@ -903,7 +913,7 @@ def spherical_coords_to_domain( elif coords["lat"].ndim == 2: message = ( "When 'dst' is a sequence containing 2-d latitude and longitude " - "coordinate constructs, 'dst_axes' must be dictionary with at " + "coordinate constructs, 'dst_axes' must be a dictionary with at " "least the keys {'X': 0, 'Y': 1} or {'X': 1, 'Y': 0}. " f"Got: {dst_axes!r}" ) @@ -1077,6 +1087,8 @@ def get_grid( .. versionadded:: 3.14.0 + :Parameters: + coord_sys: `str` The coordinate system of the source and destination grids. @@ -1089,6 +1101,11 @@ def get_grid( .. versionadded:: 3.16.2 + :Returns: + + `Grid` + The grid definition. + """ if coord_sys == "spherical": return spherical_grid( @@ -1181,6 +1198,16 @@ def spherical_grid( The grid definition. """ + domain = f.copy() + + # Create any implied lat/lon coordinates in-place + try: + # Try to convert a HEALPix grid to a UGRID grid (which will + # create 1-d lat/lon coordinates) + f.healpix_to_ugrid(inplace=True) + except ValueError: + f.create_latlon_coordinates(cache=False, inplace=True) + data_axes = f.constructs.data_axes() dim_coords_1d = False @@ -1324,7 +1351,7 @@ def spherical_grid( and (x_size == 1 or y_size == 1) ): raise ValueError( - f"Neither the X nor Y dimensions of the {name} field" + f"Neither the X nor Y dimensions of the {name} field " f"{f!r} can be of size 1 for spherical {method!r} regridding." ) @@ -1484,6 +1511,7 @@ def spherical_grid( z=z, ln_z=ln_z, z_index=z_index, + domain=domain, ) set_grid_type(grid) @@ -1501,8 +1529,8 @@ def Cartesian_grid(f, name=None, method=None, axes=None, z=None, ln_z=None): :Parameters: f: `Field` or `Domain` - The field or domain construct from which to get the - coordinates. + The field or domain construct from which to get the + coordinates. name: `str` A name to identify the grid. @@ -1529,6 +1557,8 @@ def Cartesian_grid(f, name=None, method=None, axes=None, z=None, ln_z=None): The grid definition. """ + domain = f.copy() + if not axes: if name == "source": raise ValueError( @@ -1718,6 +1748,7 @@ def Cartesian_grid(f, name=None, method=None, axes=None, z=None, ln_z=None): z=z, ln_z=ln_z, z_index=z_index, + domain=domain, ) set_grid_type(grid) @@ -3092,7 +3123,7 @@ def get_mask(f, grid): return mask -def update_coordinates(src, dst, src_grid, dst_grid): +def update_coordinates(src, dst, src_grid, dst_grid, cr_map): """Update the regrid axis coordinates. Replace the existing coordinate constructs that span the regridding @@ -3113,17 +3144,24 @@ def update_coordinates(src, dst, src_grid, dst_grid): dst: `Field` or `Domain` The field or domain containing the destination grid. - src_grid: `Grid` or `Mesh` + src_grid: `Grid` The definition of the source grid. - dst_grid: `Grid` or `Mesh` + dst_grid: `Grid` The definition of the destination grid. + cr_map `dict` + The mapping of destination coordinate reference identities + to source coordinate reference identities, as output by + `update_non_coordinates`. + :Returns: `None` """ + dst = dst_grid.domain + src_axis_keys = src_grid.axis_keys dst_axis_keys = dst_grid.axis_keys @@ -3178,19 +3216,25 @@ def update_coordinates(src, dst, src_grid, dst_grid): filter_by_axis=dst_axis_keys, axis_mode="subset", todict=True ).items(): axes = [axis_map[axis] for axis in dst_data_axes[key]] - src.set_construct(coord, axes=axes) + coord_key = src.set_construct(coord, axes=axes) + + # Update the coordinates in new source coordinate references + for dst_cr_key, src_cr_key in cr_map.items(): + dst_cr = dst.coordinate_reference(dst_cr_key) + if key in dst_cr.coordinates(): + src_cr = src.coordinate_reference(src_cr_key) + src_cr.set_coordinate(coord_key) # Copy domain topology and cell connectivity constructs from the # destination grid - if dst_grid.is_mesh: - for key, topology in dst.constructs( - filter_by_type=("domain_topology", "cell_connectivity"), - filter_by_axis=dst_axis_keys, - axis_mode="exact", - todict=True, - ).items(): - axes = [axis_map[axis] for axis in dst_data_axes[key]] - src.set_construct(topology, axes=axes) + for key, topology in dst.constructs( + filter_by_type=("domain_topology", "cell_connectivity"), + filter_by_axis=dst_axis_keys, + axis_mode="exact", + todict=True, + ).items(): + axes = [axis_map[axis] for axis in dst_data_axes[key]] + src.set_construct(topology, axes=axes) def update_non_coordinates(src, dst, src_grid, dst_grid, regrid_operator): @@ -3217,9 +3261,13 @@ def update_non_coordinates(src, dst, src_grid, dst_grid, regrid_operator): :Returns: - `None` + `dict` + The mapping of destination coordinate reference identities + to source coordinate reference identities. """ + dst = dst_grid.domain + src_axis_keys = src_grid.axis_keys dst_axis_keys = dst_grid.axis_keys @@ -3312,19 +3360,28 @@ def update_non_coordinates(src, dst, src_grid, dst_grid, regrid_operator): # ---------------------------------------------------------------- # Copy selected coordinate references from the destination grid + # + # Define the mapping of destination coordinate references to + # source coordinate references # ---------------------------------------------------------------- - dst_data_axes = dst.constructs.data_axes() + cr_map = {} - for ref in dst.coordinate_references(todict=True).values(): + dst_data_axes = dst.constructs.data_axes() + for dst_key, ref in dst.coordinate_references(todict=True).items(): axes = set() for c_key in ref.coordinates(): axes.update(dst_data_axes[c_key]) if axes and set(axes).issubset(dst_axis_keys): - src.set_coordinate_reference(ref, parent=dst, strict=True) + src_cr = ref.copy() + src_cr.clear_coordinates() + src_key = src.set_construct(src_cr) + cr_map[dst_key] = src_key + + return cr_map -def update_data(src, regridded_data, src_grid): +def update_data(src, regridded_data, src_grid, dst_grid): """Insert the regridded field data. .. versionadded:: 3.16.0 @@ -3343,6 +3400,11 @@ def update_data(src, regridded_data, src_grid): src_grid: `Grid` The definition of the source grid. + dst_grid: `Grid` + The definition of the destination grid. + + .. versionadded:: NEXTVERSION + :Returns: `None` @@ -3358,6 +3420,20 @@ def update_data(src, regridded_data, src_grid): index = data_axes.index(src_grid.axis_keys[0]) for axis in src_grid.axis_keys: data_axes.remove(axis) + for key, cm in src.cell_methods(todict=True).items(): + if axis in cm.axes: + # When regridding from one axis to two axes, we + # can safely rename cell method's axis to 'area', + # otherwise we have to delete the cell method to + # allow the domain axis itself to be deleted. + if ( + len(src_grid.axis_keys) == 1 + and len(dst_grid.axis_keys) == 2 + ): + cm.change_axes({axis: "area"}, inplace=True) + else: + src.del_construct(key) + src.del_construct(axis) data_axes[index:index] = src_grid.new_axis_keys diff --git a/cf/test/test_Data.py b/cf/test/test_Data.py index 2300a62e47..7cc2a6df44 100644 --- a/cf/test/test_Data.py +++ b/cf/test/test_Data.py @@ -4703,6 +4703,29 @@ def test_Data_to_units(self): with self.assertRaises(ValueError): e.to_units("degC") + def test_Data_coarsen(self): + """Test Data.coarsen.""" + d = cf.Data(np.arange(24).reshape((4, 6)), chunks=(3, 3)) + a = d.array + + self.assertIsNone(d.coarsen(np.min, axes={}, inplace=True)) + self.assertEqual(d.shape, a.shape) + self.assertTrue((d.array == a).all()) + + e = d.coarsen(np.min, {0: 2, 1: 3}) + self.assertIsInstance(e, cf.Data) + self.assertEqual(e.shape, (2, 2)) + self.assertTrue((e.array == [[0, 3], [12, 15]]).all()) + + # Non-full caorsening neighbourhood with trim_excess=True + e = d.coarsen(np.max, {-1: 5}, trim_excess=True) + self.assertEqual(e.shape, (4, 1)) + self.assertTrue((e.array == [[4], [10], [16], [22]]).all()) + + # Non-full caorsening neighbourhood with trim_excess=False + with self.assertRaises(ValueError): + d.coarsen(np.max, {-1: 5}) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_DimensionCoordinate.py b/cf/test/test_DimensionCoordinate.py index c727e2d1e0..b913d554a1 100644 --- a/cf/test/test_DimensionCoordinate.py +++ b/cf/test/test_DimensionCoordinate.py @@ -712,6 +712,12 @@ def test_DimensionCoordinate_create_regular(self): ) self.assertEqual(longitude_decreasing_no_bounds.units, "degrees_east") + # Size 1 with bounds + d = cf.DimensionCoordinate.create_regular((-180, 180, 360)) + self.assertEqual(d.size, 1) + self.assertTrue(np.allclose(d, 0)) + self.assertTrue(np.allclose(d.bounds, [-180, 180])) + def test_DimensionCoordinate_cell_characteristics(self): """Test the `cell_characteristic` DimensionCoordinate methods.""" d = self.dim.copy() diff --git a/cf/test/test_Domain.py b/cf/test/test_Domain.py index c42d806574..5d8bd6f47d 100644 --- a/cf/test/test_Domain.py +++ b/cf/test/test_Domain.py @@ -391,6 +391,14 @@ def test_Domain_create_regular(self): np.allclose(latitude_specific.array - y_points_specific, 0) ) + # Size 1 X axis + d = cf.Domain.create_regular((-180, 180, 360), (-90, 90, 18)) + self.assertEqual(d.size, 10) + self.assertTrue(d.cyclic()) + x = d.coordinate("X") + self.assertTrue(np.allclose(x, 0)) + self.assertTrue(np.allclose(x.bounds, [-180, 180])) + def test_Domain_del_construct(self): """Test the `del_construct` Domain method.""" # Test a domain without cyclic axes. These are equivalent tests to @@ -486,6 +494,51 @@ def test_Domain_cyclic_iscyclic(self): d2.cyclic("X", iscyclic=False) self.assertTrue(d2.iscyclic("X")) + def test_Domain_create_healpix(self): + """Test Domain.create_healpix.""" + d = cf.Domain.create_healpix(0) + self.assertEqual(len(d.constructs), 3) + self.assertEqual(len(d.domain_axes()), 1) + self.assertEqual(len(d.dimension_coordinates()), 1) + self.assertEqual(len(d.auxiliary_coordinates()), 0) + self.assertEqual(len(d.coordinate_references()), 1) + + hp = d.dimension_coordinate() + self.assertEqual(hp.data._get_cached_elements(), {0: 0, 1: 1, -1: 11}) + self.assertTrue(np.array_equal(hp, np.arange(12))) + + self.assertEqual( + d.coordinate_reference().coordinate_conversion.get_parameter( + "refinement_level" + ), + 0, + ) + + d = cf.Domain.create_healpix(0, "nuniq") + self.assertTrue( + np.array_equal(d.dimension_coordinate(), np.arange(4, 16)) + ) + self.assertIsNone( + d.coordinate_reference().datum.get_parameter("earth_radius", None) + ) + self.assertIsNone( + d.coordinate_reference().coordinate_conversion.get_parameter( + "refinement_level", None + ) + ) + + for radius in (1000, cf.Data(1, "km")): + d = cf.Domain.create_healpix(0, "ring", radius=radius) + self.assertEqual( + d.coordinate_reference().datum.get_parameter("earth_radius"), + 1000, + ) + + # Bad refinement_level specification + for level in (-1, 3.14, 30, "string"): + with self.assertRaises(ValueError): + cf.Domain.create_healpix(level) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 52fc191a3b..4c6ebfa0c2 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -79,6 +79,8 @@ class FieldTest(unittest.TestCase): f0 = cf.example_field(0) f1 = cf.example_field(1) + f12 = cf.example_field(12) + f13 = cf.example_field(13) def test_Field_creation_commands(self): for f in cf.example_fields(): @@ -965,10 +967,11 @@ def test_Field_radius(self): with self.assertRaises(ValueError): f.radius() - for default in ("earth", cf.field._earth_radius): + earth_radius = cf.Data(6371229, "m") + for default in ("earth", earth_radius): r = f.radius(default=default) self.assertEqual(r.Units, cf.Units("m")) - self.assertEqual(r, cf.field._earth_radius) + self.assertEqual(r, earth_radius) a = cf.Data(1234, "m") for default in ( @@ -2969,6 +2972,12 @@ def test_Field_auxiliary_to_dimension_to_auxiliary(self): self.assertIsNone(f.dimension_to_auxiliary("Y", inplace=True)) self.assertIsNone(g.auxiliary_to_dimension("Y", inplace=True)) + f = self.f12.copy() + g = f.dimension_to_auxiliary("healpix_index") + h = g.auxiliary_to_dimension("healpix_index") + self.assertFalse(f.equals(g)) + self.assertTrue(f.equals(h)) + f = cf.read("geometry_1.nc")[0] with self.assertRaises(ValueError): @@ -3105,6 +3114,355 @@ def test_Field_to_units(self): with self.assertRaises(ValueError): g.to_units("degC") + def test_Field_healpix_indexing_scheme(self): + """Test Field.healpix_indexing_scheme.""" + # HEALPix field + f = self.f12 + + # Null change + g = f.healpix_indexing_scheme(None) + self.assertTrue(g.equals(f)) + + # Null change + g = f.healpix_indexing_scheme("nested") + self.assertTrue(g.equals(f)) + + g = f.healpix_indexing_scheme("ring") + self.assertTrue( + np.array_equal(g.coordinate("healpix_index")[:4], [13, 5, 4, 0]) + ) + h = g.healpix_indexing_scheme("nested") + self.assertTrue( + np.array_equal(h.coordinate("healpix_index")[:4], [0, 1, 2, 3]) + ) + h = g.healpix_indexing_scheme("nuniq") + self.assertTrue( + np.array_equal(h.coordinate("healpix_index")[:4], [16, 17, 18, 19]) + ) + + g = f.healpix_indexing_scheme("ring", sort=True) + self.assertTrue( + np.array_equal(g.coordinate("healpix_index")[:4], [0, 1, 2, 3]) + ) + h = g.healpix_indexing_scheme("nested", sort=False) + self.assertTrue( + np.array_equal(h.coordinate("healpix_index")[:4], [3, 7, 11, 15]) + ) + h = g.healpix_indexing_scheme("nested", sort=True) + self.assertTrue( + np.array_equal(h.coordinate("healpix_index")[:4], [0, 1, 2, 3]) + ) + + g = f.healpix_indexing_scheme("nuniq") + self.assertTrue( + np.array_equal(g.coordinate("healpix_index")[:4], [16, 17, 18, 19]) + ) + h = g.healpix_indexing_scheme("nuniq") + self.assertTrue(h.equals(g)) + + # Can't change from 'nuniq' to 'nested' + with self.assertRaises(ValueError): + self.f13.healpix_indexing_scheme("nested") + + # Non-HEALPix field + with self.assertRaises(ValueError): + self.f0.healpix_indexing_scheme("ring") + + def test_Field_healpix_to_ugrid(self): + """Test Field.healpix_to_ugrid.""" + # HEALPix field + f = self.f12.copy() + + u = f.healpix_to_ugrid() + self.assertEqual(len(u.domain_topologies()), 1) + self.assertEqual(len(u.auxiliary_coordinates()), 2) + + topology = u.domain_topology().normalise() + self.assertEqual(np.unique(topology).size, 53) + self.assertTrue( + np.array_equal( + topology[:4], + [ + [14, 11, 13, 16], + [21, 14, 16, 20], + [8, 7, 11, 14], + [9, 8, 14, 21], + ], + ) + ) + + # North pole + self.assertTrue(np.allclose(topology[3:16:4, 0], 9)) + + # South pole + self.assertTrue(np.allclose(topology[32:48:4, 2], 3)) + + self.assertIsNone(f.healpix_to_ugrid(inplace=True)) + self.assertEqual(len(f.domain_topologies()), 1) + + self.assertEqual(len(u.coordinate_references()), 1) + cr = u.coordinate_reference() + self.assertEqual(cr.identity(), "grid_mapping_name:latitude_longitude") + + # Non-HEALPix field + with self.assertRaises(ValueError): + self.f0.healpix_to_ugrid() + + def test_Field_create_latlon_coordinates(self): + """Test Field.create_latlon_coordinates.""" + # ------------------------------------------------------------ + # HEALPix field + # ------------------------------------------------------------ + f = self.f12.copy() + self.assertEqual(len(f.auxiliary_coordinates()), 0) + self.assertEqual(len(f.dimension_coordinates("healpix_index")), 1) + + g = f.create_latlon_coordinates() + self.assertEqual(len(g.auxiliary_coordinates("X", "Y")), 2) + self.assertIsNone(f.create_latlon_coordinates(inplace=True)) + self.assertTrue(f.equals(g)) + + g = self.f12.healpix_indexing_scheme("nuniq") + g.create_latlon_coordinates(inplace=True) + for c in ("latitude", "longitude"): + self.assertTrue( + g.auxiliary_coordinate(c).equals(f.auxiliary_coordinate(c)) + ) + + # pole_longitude. Note that bounds index 0 is the + # northern-most vertex, and bounds index 2 is the + # southern-most vertex. + f = self.f12 + g = f.create_latlon_coordinates() + longitude = g.auxiliary_coordinate("X").bounds.array + # North pole + self.assertTrue( + np.allclose(longitude[3:16:4, 0], longitude[3:16:4, 2]) + ) + # South pole + self.assertTrue( + np.allclose(longitude[32:48:4, 2], longitude[32:48:4, 0]) + ) + + g = f.create_latlon_coordinates(pole_longitude=3.14) + longitude = g.auxiliary_coordinate("X").bounds.array + # North pole + self.assertTrue(np.allclose(longitude[3:16:4, 0], 3.14)) + # South pole + self.assertTrue(np.allclose(longitude[32:48:4, 2], 3.14)) + + # Multi-Order Coverage (MOC) grid with refinement levels 1 and + # 2 + m = self.f13 + m = m.create_latlon_coordinates() + + l1 = cf.Domain.create_healpix(1) + l2 = cf.Domain.create_healpix(2) + l1.create_latlon_coordinates(inplace=True) + l2.create_latlon_coordinates(inplace=True) + + for c in ("latitude", "longitude"): + mc = m.auxiliary_coordinate(c) + self.assertTrue(mc[:16].equals(l2.auxiliary_coordinate(c)[:16])) + self.assertTrue(mc[16:].equals(l1.auxiliary_coordinate(c)[4:])) + + def test_Field_healpix_subspace(self): + """Test Field.subspace for HEALPix grids""" + f = self.f12 + + index = [47, 3, 2, 0] + g = f.subspace(healpix_index=index) + self.assertTrue(np.array_equal(g.coordinate("healpix_index"), index)) + + g = f.subspace(X=cf.wi(40, 70), Y=cf.wi(-20, 30)) + self.assertTrue( + np.array_equal(g.coordinate("healpix_index"), [0, 22, 35]) + ) + g.create_latlon_coordinates(inplace=True) + self.assertTrue(np.allclose(g.coordinate("X"), [45.0, 67.5, 45.0])) + self.assertTrue( + np.allclose( + g.coordinate("Y"), [19.47122063449069, 0.0, -19.47122063449069] + ) + ) + + g = f.subspace(healpix_index=cf.locate(20, 46)) + self.assertEqual(g.coordinate("healpix_index").array, 0) + + g = f.subspace(healpix_index=cf.locate(20, 1)) + self.assertEqual(g.coordinate("healpix_index").array, 19) + + g = f.subspace(healpix_index=cf.locate(20, [1, 46])) + self.assertTrue(np.array_equal(g.coordinate("healpix_index"), [0, 19])) + + f = f.healpix_indexing_scheme("ring") + g = f.subspace(healpix_index=cf.locate(20, [1, 46])) + self.assertTrue( + np.array_equal(g.coordinate("healpix_index"), [13, 12]) + ) + + def test_Field_healpix_decrease_refinement_level(self): + """Test Field.healpix_decrease_refinement_level.""" + f = self.f12 + + # No change + for level in (None, 1): + g = f.healpix_decrease_refinement_level(level, "mean") + self.assertTrue(g.equals(f)) + + # Can't change refinement level when a larger cell is only + # partially covered by original cells + with self.assertRaises(ValueError): + f[:, 1:].healpix_decrease_refinement_level(0, "mean") + + g = f.healpix_decrease_refinement_level(0, "mean") + self.assertTrue( + np.array_equal(g.coord("healpix_index"), np.arange(12)) + ) + + for method, first_value in zip( + ( + "maximum", + "mean", + "minimum", + "standard_deviation", + "variance", + "sum", + "median", + ), + (293.5, 289.15, 285.3, 3.44201976, 11.8475, 1156.6, 288.9), + ): + g = f.healpix_decrease_refinement_level(0, method) + self.assertTrue(np.isclose(g[0, 0], first_value)) + + # Bad methods + for method in ("point", "range", "bad method", 3.14): + with self.assertRaises(ValueError): + f.healpix_decrease_refinement_level(0, method) + + def range_func(a, axis=None): + return np.max(a, axis=axis) - np.min(a, axis=axis) + + g = f.healpix_decrease_refinement_level(0, "range", range_func) + self.assertTrue(np.isclose(g[0, 0], 8.2)) + + def my_mean(a, axis=None): + return np.mean(a, axis=axis) + + g = f.healpix_decrease_refinement_level(0, "mean", my_mean) + self.assertTrue(np.isclose(g[0, 0], 289.15)) + + f = f.healpix_indexing_scheme("ring") + g = f.healpix_decrease_refinement_level(0, "maximum") + self.assertTrue( + np.array_equal(g.coord("healpix_index"), np.arange(12)) + ) + with self.assertRaises(ValueError): + f.healpix_decrease_refinement_level(0, "maximum", conform=False) + + f = f.healpix_indexing_scheme("ring", sort=True) + f = f.healpix_indexing_scheme("nested") + + g = f.healpix_decrease_refinement_level(0, "mean", conform=True) + self.assertTrue( + np.array_equal(g.coord("healpix_index"), np.arange(12)) + ) + + # Check that lat/lon coords get created when they're present + # in the original field + f = f.create_latlon_coordinates() + g = f.healpix_decrease_refinement_level(0, "mean") + self.assertEqual(g.auxiliary_coordinate("latitude"), (12,)) + self.assertEqual(g.auxiliary_coordinate("longitude"), (12,)) + + # Can't change refinement level when the HEALPix indices are + # not strictly monotonically increasing + with self.assertRaises(ValueError): + f.healpix_decrease_refinement_level(0, "mean", conform=False) + + # Bad results when check_healpix_index=False + h = f.healpix_decrease_refinement_level( + 0, "mean", conform=False, check_healpix_index=False + ) + self.assertFalse( + np.array_equal(h.coord("healpix_index"), np.arange(12)) + ) + + # Bad 'level' parameter + for level in (-1, 0.785, np.array(1), 2, "string"): + with self.assertRaises(ValueError): + f.healpix_decrease_refinement_level(level, "mean") + + # Can't change refinement level for a 'nuniq' field + # TODOHEALPIX + with self.assertRaises(ValueError): + self.f13.healpix_decrease_refinement_level(0, "mean") + + # Non-HEALPix field + with self.assertRaises(ValueError): + self.f0.healpix_decrease_refinement_level(0, "mean") + + def test_Field_healpix_increase_refinement_level(self): + """Test Field.healpix_increase_refinement_level.""" + f = self.f12 + f = f.rechunk((None, 17)) + f.coordinate("healpix_index").rechunk(13, inplace=True) + + g = f.healpix_increase_refinement_level(2, "intensive") + self.assertTrue( + np.array_equal(g.coord("healpix_index"), np.arange(192)) + ) + + self.assertEqual(g.shape, (2, 192)) + self.assertEqual(g.array.shape, (2, 192)) + + # Check selected data values for intensive and extensive + # quantities + n = 4 ** (2 - 1) + for i in (0, 1, 24, 46, 47): + self.assertTrue( + np.allclose(g[:, i * n : (i + 1) * n], f[:, i : i + 1]) + ) + + g = f.healpix_increase_refinement_level(2, "extensive") + for i in (0, 1, 24, 46, 47): + self.assertTrue( + np.allclose(g[:, i * n : (i + 1) * n], f[:, i : i + 1] / n) + ) + + # Cached values + f.coordinate("healpix_index").data._del_cached_elements() + with cf.chunksize(256 * (2**30)): + g = f.healpix_increase_refinement_level(16, "intensive") + self.assertEqual(g.data.npartitions, 6) + self.assertEqual( + g.coordinate("healpix_index").data._get_cached_elements(), {} + ) + # Create cached elements + _ = str(f.coordinate("healpix_index").data) + g = f.healpix_increase_refinement_level(16, "intensive") + self.assertEqual( + g.coordinate("healpix_index").data._get_cached_elements(), + {0: 0, 1: 1, -1: 51539607551}, + ) + + # Bad 'quantity' parameter + with self.assertRaises(ValueError): + f.healpix_increase_refinement_level(2, "bad quantity") + + # Bad 'level' parameter + for level in (-1, 0, np.array(2), 3.14, 30, "string"): + with self.assertRaises(ValueError): + f.healpix_increase_refinement_level(level, "intensive") + + # Can't change refinement level for a 'nuniq' field + with self.assertRaises(ValueError): + self.f13.healpix_increase_refinement_level(2, "intensive") + + # Non-HEALPix field + with self.assertRaises(ValueError): + self.f0.healpix_increase_refinement_level(2, "intensive") + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_collapse.py b/cf/test/test_collapse.py index c7c866c8b6..22db0bceb8 100644 --- a/cf/test/test_collapse.py +++ b/cf/test/test_collapse.py @@ -3,7 +3,7 @@ import os import unittest -import numpy +import numpy as np faulthandler.enable() # to debug seg faults and timeouts @@ -350,7 +350,7 @@ def test_Field_collapse(self): self.assertEqual(g.shape, (12, 4, 5)) for m in range(1, 13): - a = numpy.empty((5, 4, 5)) + a = np.empty((5, 4, 5)) for i, year in enumerate( f.subspace(T=cf.month(m)).coord("T").year.unique() ): @@ -360,7 +360,7 @@ def test_Field_collapse(self): a[i] = x.array a = a.min(axis=0) - self.assertTrue(numpy.allclose(a, g.array[m % 12])) + self.assertTrue(np.allclose(a, g.array[m % 12])) g = f.collapse("T: mean", group=360) @@ -791,6 +791,30 @@ def test_Field_collapse_non_positive_weights(self): # compute time g.array + def test_Field_collapse_HEALPix(self): + """Test HEALPix collapses.""" + f0 = cf.example_field(12) + f1 = cf.example_field(13) + + g0 = f0.collapse("area: mean", weights=False) + g1 = f1.collapse("area: mean", weights=False) + self.assertFalse(np.allclose(g0, g1)) + + g0 = f0.collapse("area: mean", weights=True) + g1 = f1.collapse("area: mean", weights=True) + self.assertTrue(g1.equals(g0)) + self.assertTrue( + g0.coordinate_reference("grid_mapping_name:latitude_longitude") + ) + + g0 = f0[:, 0].collapse("area: mean", weights=True) + g1 = f1[:, :4].collapse("area: mean", weights=True) + self.assertTrue(np.allclose(g0, g1)) + self.assertTrue(g0.coordinate_reference("grid_mapping_name:healpix")) + self.assertTrue( + g1.coordinate_reference("grid_mapping_name:latitude_longitude") + ) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_functions.py b/cf/test/test_functions.py index 1a26a1afdc..ca25b4607d 100644 --- a/cf/test/test_functions.py +++ b/cf/test/test_functions.py @@ -430,6 +430,32 @@ def test_normalize_slice(self): with self.assertRaises(IndexError): cf.normalize_slice(index, 8, cyclic=True) + def test_locate(self): + """Test cf.locate""" + # HEALPix + f = cf.example_field(12) + self.assertEqual(cf.locate(20, 90, f), 23) + self.assertEqual(cf.locate(-70, 90, f), 36) + self.assertEqual(cf.locate(20, [280, 280.001], f), 31) + + self.assertTrue(np.array_equal(cf.locate([-70, 20], 90, f), [23, 36])) + self.assertTrue( + np.array_equal(cf.locate([-70, 20], [90, 280], f), [31, 36]) + ) + self.assertTrue( + np.array_equal(cf.locate([-70, 20], [90, 280])(f), [31, 36]) + ) + + # Bad latitudes + for lat in (-91, 91): + with self.assertRaises(ValueError): + cf.locate(lat, 30, f) + + # Invalid grid types (regulat lat/lon, geometry, UGRID) + for f in cf.example_fields(0, 6, 8): + with self.assertRaises(ValueError): + cf.locate(60, 30, f) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_regrid_healpix.py b/cf/test/test_regrid_healpix.py new file mode 100644 index 0000000000..475cc0ea4e --- /dev/null +++ b/cf/test/test_regrid_healpix.py @@ -0,0 +1,120 @@ +import datetime +import faulthandler +import os +import unittest + +faulthandler.enable() # to debug seg faults and timeouts + +import numpy as np + +import cf + +esmpy_imported = True +try: + import esmpy # noqa: F401 +except ImportError: + esmpy_imported = False + + +all_methods = ( + "linear", + "conservative", + "conservative_2nd", + "nearest_dtos", + "nearest_stod", + "patch", +) + + +# Set numerical comparison tolerances +atol = 0 +rtol = 0 + + +class RegridMeshTest(unittest.TestCase): + # Get the test source and destination fields + src_mesh_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "ugrid_global_1.nc" + ) + dst_mesh_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "ugrid_global_2.nc" + ) + src_mesh = cf.read(src_mesh_file)[0] + dst_mesh = cf.read(dst_mesh_file)[0] + + def setUp(self): + """Preparations called immediately before each test method.""" + # Disable log messages to silence expected warnings + cf.log_level("DISABLE") + # Note: to enable all messages for given methods, lines or + # calls (those without a 'verbose' option to do the same) + # e.g. to debug them, wrap them (for methods, start-to-end + # internally) as follows: + # cfdm.log_level('DEBUG') + # < ... test code ... > + # cfdm.log_level('DISABLE') + + @unittest.skipUnless(esmpy_imported, "Requires esmpy package.") + def test_Field_regrid_mesh_to_healpix(self): + # Check that UGRID -> healpix is the same as UGRID -> UGRUD + self.assertFalse(cf.regrid_logging()) + + dst = cf.Domain.create_healpix(3) # 768 cells + dst_ugrid = dst.healpix_to_ugrid() + src = self.src_mesh.copy() + + for src_masked in (False, True): + if src_masked: + src = src.copy() + src[100:200] = cf.masked + + # Loop over whether or not to use the destination grid + # masked points + for method in all_methods: + x = src.regrids(dst, method=method) + y = src.regrids(dst_ugrid, method=method) + a = x.array + b = y.array + self.assertTrue(np.allclose(b, a, atol=atol, rtol=rtol)) + + if np.ma.isMA(a): + self.assertTrue((b.mask == a.mask).all()) + + # Check that the result is a HEALPix grid + self.assertTrue(cf.healpix.healpix_info(x)) + + @unittest.skipUnless(esmpy_imported, "Requires esmpy package.") + def test_Field_regrid_healpix_to_mesh(self): + # Check that healpix -> UGRID is the same as UGRID -> UGRUD + self.assertFalse(cf.regrid_logging()) + + src = cf.Field(source=cf.Domain.create_healpix(3)) # 768 cells + src.set_data(np.arange(768)) + + src_ugrid = src.healpix_to_ugrid() + + dst = self.dst_mesh.copy() + + for src_masked in (False, True): + if src_masked: + src[100:200] = cf.masked + src_ugrid[100:200] = cf.masked + + # Loop over whether or not to use the destination grid + # masked points + for method in all_methods: + x = src.regrids(dst, method=method) + y = src_ugrid.regrids(dst, method=method) + a = x.array + b = y.array + self.assertTrue(np.allclose(b, a, atol=atol, rtol=rtol)) + + if np.ma.isMA(a): + self.assertTrue((b.mask == a.mask).all()) + + +if __name__ == "__main__": + print("Run date:", datetime.datetime.now()) + cf.environment() + print("") + unittest.main(verbosity=2) diff --git a/cf/test/test_weights.py b/cf/test/test_weights.py index bd202672a3..5043a6b796 100644 --- a/cf/test/test_weights.py +++ b/cf/test/test_weights.py @@ -217,6 +217,15 @@ def test_weights_polygon_area_ugrid(self): self.assertTrue((w.array == correct_weights).all()) self.assertEqual(w.Units, cf.Units("m2")) + # For a UGRID derived from a global HEALPix grid, check that + # the global sum of cell areas is correct + u = cf.example_field(12).healpix_to_ugrid() + u_weights = u.weights( + measure=True, components=True, great_circle=True + )[(1,)] + global_area = 4 * np.pi * (u.radius() ** 2) + self.assertTrue(np.allclose(u_weights.sum(), global_area)) + def test_weights_line_length_geometry(self): # Spherical line geometry gls = gps.copy() @@ -336,6 +345,23 @@ def test_weights_exceptions(self): ): f.weights("area") + def test_weights_healpix(self): + """Test HEALPix weights.""" + # HEALPix grid with Multi-Order Coverage (a combination of + # refinement level 1 and 2 cells) + f = cf.example_field(13) + + w = f.weights(components=True)[(1,)].array + self.assertTrue(np.allclose(w[:16], 1 / (4**2))) + self.assertTrue(np.allclose(w[16:], 1 / (4**1))) + + w = f.weights(measure=True, components=True)[(1,)].array + x = 4 * np.pi * (f.radius() ** 2) / 12 + self.assertTrue(np.allclose(w[:16], x / (4**2))) + self.assertTrue(np.allclose(w[16:], x / (4**1))) + # Total global area + self.assertTrue(np.allclose(w.sum(), x * 12)) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/weights.py b/cf/weights.py index c0d50f3fbc..0a1ef29e44 100644 --- a/cf/weights.py +++ b/cf/weights.py @@ -1,4 +1,5 @@ import cfdm +import numpy as np from .mixin_container import Container @@ -635,8 +636,6 @@ def _polygon_area_geometry(cls, f, x, y, aux_X, aux_Y, n_nodes, spherical): square metres. """ - import numpy as np - x_units = x.Units y_units = y.Units @@ -798,8 +797,6 @@ def _polygon_area_ugrid(cls, f, x, y, n_nodes, spherical): metres. """ - import numpy as np - y_units = y.Units n_nodes0 = n_nodes.item(0) @@ -1259,7 +1256,6 @@ def _interior_angles(cls, f, lon, lat, interior_rings=None): units of radians. """ - import numpy as np from .query import lt @@ -1916,3 +1912,191 @@ def _spherical_polygon_areas(cls, f, x, y, N, interior_rings=None): areas = interior_angles.sum(-1, squeeze=True) - (N - 2) * pi areas.override_units(Units("m2"), inplace=True) return areas + + @classmethod + def healpix_area( + cls, + f, + domain_axis, + weights, + weights_axes, + auto=False, + measure=False, + radius=None, + return_areas=False, + methods=False, + ): + """Creates area weights for HEALPix cells. + + .. versionadded:: NEXTVERSION + + :Parameters: + + f: `Field` + The field for which the weights are being created. + + domain_axis: `str` or `None` + If set to a domain axis identifier + (e.g. ``'domainaxis1'``) then only accept cells that + recognise the given axis. If `None` then the cells may + span any axis. + + {{weights weights: `dict`}} + + {{weights weights_axes: `set`}} + + {{weights auto: `bool`, optional}} + + {{weights measure: `bool`, optional}} + + {{radius: optional}} + + {{weights methods: `bool`, optional}} + + :Returns: + + `bool` or `Data` + `True` if weights were created, otherwise `False`. If + *return_areas* is True and weights were created, then + the weights are returned. + + """ + axis = f.domain_axis("healpix_index", key=True, default=None) + if axis is None: + if auto: + return False + + if domain_axis is None: + raise ValueError("No HEALPix axis") + + raise ValueError( + "No HEALPix cells for " + f"{f.constructs.domain_axis_identity(domain_axis)!r} axis" + ) + + if domain_axis is not None and domain_axis != axis: + if auto: + return False + + raise ValueError( + "No HEALPix cells for " + f"{f.constructs.domain_axis_identity(domain_axis)!r} axis" + ) + + if axis in weights_axes: + if auto: + return False + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + cr = f.coordinate_reference("grid_mapping_name:healpix", default=None) + if cr is None: + # No healpix grid mapping + if auto: + return False + + raise ValueError( + "Can't create weights: No HEALPix grid mapping for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + parameters = cr.coordinate_conversion.parameters() + indexing_scheme = parameters.get("indexing_scheme") + if indexing_scheme not in ("nested", "ring", "nuniq"): + if auto: + return False + + raise ValueError( + "Can't create weights: Invalid HEALPix indexing_scheme for " + f"{f.constructs.domain_axis_identity(axis)!r} axis: " + f"{indexing_scheme!r}" + ) + + refinement_level = parameters.get("refinement_level") + if refinement_level is None and indexing_scheme != "nuniq": + # No refinement_level + if auto: + return False + + raise ValueError( + "Can't create weights: No HEALPix refinement_level for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + if measure and not methods and radius is not None: + radius = f.radius(default=radius) + + # Create weights for 'nuniq' indexed cells + if indexing_scheme == "nuniq": + if methods: + weights[(axis,)] = "HEALPix Multi-Order Coverage" + return True + + healpix_index = f.coordinate( + "healpix_index", + filter_by_axis=(axis,), + axis_mode="exact", + default=None, + ) + if healpix_index is None: + if auto: + return False + + raise ValueError( + "Can't create weights: Missing healpix_index coordinates" + ) + + if measure: + units = radius.Units**2 + r = radius.array + else: + units = "1" + r = None + + from .data.dask_utils import cf_healpix_weights + + dx = healpix_index.to_dask_array() + dx = dx.map_blocks( + cf_healpix_weights, + meta=np.array((), dtype="float64"), + indexing_scheme="nuniq", + measure=measure, + radius=r, + ) + area = f._Data(dx, units=units, copy=False) + + if return_areas: + return area + + weights[(axis,)] = area + weights_axes.add(axis) + return True + + # Still here? Then create weights for 'nested' or 'ring' + # indexed cells. + if methods: + if measure: + weights[(axis,)] = "HEALPix equal area" + + return True + + if not measure: + # Weights are all equal, so no need to create any. + return True + + r2 = radius**2 + area = f._Data.full( + (f.domain_axis(axis).get_size(),), + np.pi * float(r2) / (3.0 * (4**refinement_level)), + units=r2.Units, + ) + + if return_areas: + return area + + weights[(axis,)] = area + weights_axes.add(axis) + return True diff --git a/docs/source/class/cf.Data.rst b/docs/source/class/cf.Data.rst index b89f262de9..bd4503806b 100644 --- a/docs/source/class/cf.Data.rst +++ b/docs/source/class/cf.Data.rst @@ -134,6 +134,7 @@ Changing data shape ~cf.Data.flatten ~cf.Data.reshape + ~cf.Data.coarsen Transpose-like operations ^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/source/class/cf.Domain.rst b/docs/source/class/cf.Domain.rst index 7f75c0b3ac..9b2a898185 100644 --- a/docs/source/class/cf.Domain.rst +++ b/docs/source/class/cf.Domain.rst @@ -180,7 +180,9 @@ Miscellaneous ~cf.Domain.apply_masking ~cf.Domain.climatological_time_axes ~cf.Domain.copy + ~cf.Domain.create_latlon_coordinates ~cf.Domain.create_regular + ~cf.Domain.create_healpix ~cf.Domain.creation_commands ~cf.Domain.equals ~cf.Domain.fromconstructs @@ -191,7 +193,14 @@ Miscellaneous ~cf.Domain.get_original_filenames ~cf.Domain.close ~cf.Domain.persist - ~cf.Domain.uncompress + ~cf.Domain.uncompresshas_bounds + ~cf.Domain.has_data + ~cf.Domain.has_geometry + ~cf.Domain.apply_masking + ~cf.Domain.get_original_filenames + ~cf.Domain.close + ~cf.Domain.persist + ~cf.Domain.radius Domain axes ----------- @@ -221,6 +230,19 @@ Subspacing ~cf.Domain.indices ~cf.Domain.subspace +HEALPix grids +------------- + +.. autosummary:: + :nosignatures: + :toctree: ../method/ + :template: method.rst + + ~cf.Domain.healpix_info + ~cf.Domain.healpix_indexing_scheme + ~cf.Domain.healpix_to_ugrid + ~cf.Domain.create_healpix + NetCDF ------ diff --git a/docs/source/class/cf.Field.rst b/docs/source/class/cf.Field.rst index 1c307f362e..ede916be8d 100644 --- a/docs/source/class/cf.Field.rst +++ b/docs/source/class/cf.Field.rst @@ -386,6 +386,7 @@ Miscellaneous :template: method.rst ~cf.Field.copy + ~cf.Field.create_latlon_coordinates ~cf.Field.compute_vertical_coordinates ~cf.Field.dataset_compliance ~cf.Field.equals @@ -713,6 +714,22 @@ Regridding operations ~cf.Field.regridc ~cf.Field.regrids + ~cf.Field.healpix_decrease_refinement_level + ~cf.Field.healpix_increase_refinement_level + +HEALPix grids +------------- + +.. autosummary:: + :nosignatures: + :toctree: ../method/ + :template: method.rst + + ~cf.Field.healpix_info + ~cf.Field.healpix_indexing_scheme + ~cf.Field.healpix_to_ugrid + ~cf.Field.healpix_decrease_refinement_level + ~cf.Field.healpix_increase_refinement_level Date-time operations -------------------- diff --git a/docs/source/field_analysis.rst b/docs/source/field_analysis.rst index 32018df19c..08566f7442 100644 --- a/docs/source/field_analysis.rst +++ b/docs/source/field_analysis.rst @@ -1395,7 +1395,7 @@ on the surface of the sphere, rather than in :ref:`Euclidean space Spherical regridding can occur between source and destination grids that comprise any pairing of `Latitude-longitude`_, `Rotated latitude-longitude`_, `Plane projection`_, `Tripolar`_, `UGRID mesh`_, -and `DSG feature type`_ coordinate systems. +`HEALPix_`, and `DSG feature type`_ coordinate systems. The most convenient usage is when the destination domain exists in another field construct. In this case, all you need to specify is the diff --git a/docs/source/function.rst b/docs/source/function.rst index 8892f325c1..0cdff526c6 100644 --- a/docs/source/function.rst +++ b/docs/source/function.rst @@ -101,6 +101,7 @@ Condition constructors :template: function.rst cf.contains + cf.contains_latlon cf.cellsize cf.cellgt cf.cellge diff --git a/docs/source/installation.rst b/docs/source/installation.rst index b029966c48..e1da50c126 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -290,6 +290,15 @@ environments for which these features are not required. must use a copy of the ``pyfive`` branch of the https://github.com/NCAS-CMS/PyActiveStorage repository. +.. rubric:: HEALPix manipulations + +* `healpix `_, version 2025.1 or + newer. This package is not required to read and write HEALPix + datasets, but may be needed for particular manipulations with + HEALPix grids, such as creating latitude and longitude coordinates, + regridding, some changes to the refinement level, and some + collapses. + ---- .. _Tests: diff --git a/docs/source/recipes/plot_08_recipe.py b/docs/source/recipes/plot_08_recipe.py index 63427f62a7..6045f51448 100644 --- a/docs/source/recipes/plot_08_recipe.py +++ b/docs/source/recipes/plot_08_recipe.py @@ -9,11 +9,10 @@ # 1. Import cf-python, cf-plot, numpy and scipy.stats: import cfplot as cfp -import cf - import numpy as np import scipy.stats as stats +import cf # %% # 2. Three functions are defined: diff --git a/docs/source/recipes/plot_12_recipe.py b/docs/source/recipes/plot_12_recipe.py index b09db0b29f..5304194b19 100644 --- a/docs/source/recipes/plot_12_recipe.py +++ b/docs/source/recipes/plot_12_recipe.py @@ -13,8 +13,8 @@ # %% # 1. Import cf-python, cf-plot and matplotlib.pyplot: -import matplotlib.pyplot as plt import cfplot as cfp +import matplotlib.pyplot as plt import cf diff --git a/docs/source/recipes/plot_13_recipe.py b/docs/source/recipes/plot_13_recipe.py index bf0398713e..9b658597d8 100644 --- a/docs/source/recipes/plot_13_recipe.py +++ b/docs/source/recipes/plot_13_recipe.py @@ -18,13 +18,11 @@ # in next steps. import cartopy.crs as ccrs -import matplotlib.patches as mpatches - import cfplot as cfp +import matplotlib.patches as mpatches import cf - # %% # 2. Read and select the SST by index and look at its contents: sst = cf.read("~/recipes/ERA5_monthly_averaged_SST.nc")[0] diff --git a/docs/source/recipes/plot_17_recipe.py b/docs/source/recipes/plot_17_recipe.py index c94769e2ba..a66c90b518 100644 --- a/docs/source/recipes/plot_17_recipe.py +++ b/docs/source/recipes/plot_17_recipe.py @@ -11,8 +11,8 @@ # %% # 1. Import cf-python and cf-plot: -import matplotlib.pyplot as plt import cfplot as cfp +import matplotlib.pyplot as plt import cf diff --git a/docs/source/recipes/plot_18_recipe.py b/docs/source/recipes/plot_18_recipe.py index f0eae36e35..3beb9d0db9 100644 --- a/docs/source/recipes/plot_18_recipe.py +++ b/docs/source/recipes/plot_18_recipe.py @@ -10,15 +10,15 @@ """ +import cfplot as cfp + # %% # 1. Import cf-python, cf-plot and other required packages: import matplotlib.pyplot as plt import scipy.stats.mstats as mstats -import cfplot as cfp import cf - # %% # 2. Read the data in and unpack the Fields from FieldLists using indexing. # In our example We are investigating the influence of the land height on diff --git a/docs/source/recipes/plot_19_recipe.py b/docs/source/recipes/plot_19_recipe.py index 02d493dc21..ceb9db1c5c 100644 --- a/docs/source/recipes/plot_19_recipe.py +++ b/docs/source/recipes/plot_19_recipe.py @@ -9,10 +9,11 @@ maxima. """ +import cfplot as cfp + # %% # 1. Import cf-python, cf-plot and other required packages: import matplotlib.pyplot as plt -import cfplot as cfp import cf diff --git a/docs/source/recipes/plot_22_recipe.py b/docs/source/recipes/plot_22_recipe.py index 377313c899..fe329cda9d 100644 --- a/docs/source/recipes/plot_22_recipe.py +++ b/docs/source/recipes/plot_22_recipe.py @@ -11,10 +11,11 @@ # %% # 1. Import cf-python, Dask.array, NumPy, and Matplotlib: -import cf import dask.array as da -import numpy as np import matplotlib.pyplot as plt +import numpy as np + +import cf # %% # 2. Read the field constructs and load the wind speed component fields: diff --git a/docs/source/recipes/plot_23_recipe.py b/docs/source/recipes/plot_23_recipe.py index 2499b0d875..29537803af 100644 --- a/docs/source/recipes/plot_23_recipe.py +++ b/docs/source/recipes/plot_23_recipe.py @@ -18,12 +18,12 @@ # sphinx_gallery_thumbnail_number = 2 # sphinx_gallery_end_ignore -import matplotlib.pyplot as plt import cfplot as cfp -import cf - -import numpy as np import dask.array as da +import matplotlib.pyplot as plt +import numpy as np + +import cf # %% # 2. Read example data field constructs, and set region for our plots: diff --git a/setup.py b/setup.py index f5ddbed1f7..d97787d01f 100755 --- a/setup.py +++ b/setup.py @@ -186,6 +186,10 @@ def compile(): * write and append field constructs to netCDF datasets on disk, +* read, create, and manipulate UGRID mesh topologies, + +* read, write, and manipulate HEALPix grids, + * read, write, and create coordinates defined by geometry cells, * read netCDF and CDL datasets containing hierarchical groups, @@ -217,8 +221,8 @@ def compile(): * perform histogram, percentile and binning operations on field constructs, -* regrid structured grid, mesh and DSG field constructs with - (multi-)linear, nearest neighbour, first- and second-order +* regrid structured grid, UGRID, HEALPix, and DSG field constructs + with (multi-)linear, nearest neighbour, first- and second-order conservative and higher order patch recovery methods, including 3-d regridding, and large-grid support,