diff --git a/Changelog.rst b/Changelog.rst index 1d8a4bcb34..c6e6b25d5c 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,6 +3,12 @@ version NEXT **2024-??-??** +* Added spherical regridding to discrete sampling geometry destination + grids (https://github.com/NCAS-CMS/cf-python/issues/716) +* Added 3-d spherical regridding to `cf.Field.regrids`, and the option + to regrid the vertical axis in logarithmic coordinates to + `cf.Field.regrids` and `cf.Field.regridc` + (https://github.com/NCAS-CMS/cf-python/issues/715) * Fix misleading error message when it is not possible to create area weights requested from `cf.Field.collapse` (https://github.com/NCAS-CMS/cf-python/issues/731) diff --git a/cf/data/dask_regrid.py b/cf/data/dask_regrid.py index 3b1326032b..510591c339 100644 --- a/cf/data/dask_regrid.py +++ b/cf/data/dask_regrid.py @@ -54,8 +54,7 @@ def regrid( source grid cells i. Note that it is typical that for a given j most w_ji will be zero, reflecting the fact only a few source grid cells intersect a particular destination - grid cell. I.e. *weights* is typically a very sparse - matrix. + grid cell. I.e. *weights* is usually a very sparse matrix. If the destination grid has masked cells, either because it spans areas outside of the source grid, or by selection @@ -314,6 +313,27 @@ def regrid( # [0,1,3,4,2] raxis0, raxis = axis_order[-2:] axis_order = [i if i <= raxis else i - 1 for i in axis_order[:-1]] + elif n_src_axes == 3 and n_dst_axes == 1: + # The regridding operation decreased the number of data axes + # by 2 => modify 'axis_order' to remove the removed axes. + # + # E.g. regular Z-lat-lon -> DSG could change 'axis_order' from + # [0,2,5,1,3,4] to [0,2,3,1], or [0,2,4,5,3,1] to + # [0,1,2,3] + raxis0, raxis1 = axis_order[-2:] + if raxis0 > raxis1: + raxis0, raxis1 = raxis1, raxis0 + + new = [] + for i in axis_order[:-2]: + if i <= raxis0: + new.append(i) + elif raxis0 < i <= raxis1: + new.append(i - 1) + else: + new.append(i - 2) + + axis_order = new elif n_src_axes != n_dst_axes: raise ValueError( f"Can't (yet) regrid from {n_src_axes} dimensions to " diff --git a/cf/data/data.py b/cf/data/data.py index 9c730b72d7..fd84d14bb4 100644 --- a/cf/data/data.py +++ b/cf/data/data.py @@ -3840,7 +3840,7 @@ def _regrid( The positions of the regrid axes in the data, given in the relative order expected by the regrid operator. For spherical regridding this order is [Y, - X]. + X] or [Z, Y, X]. *Parameter example:* ``[2, 3]`` @@ -3902,7 +3902,7 @@ def _regrid( new_axis.extend(range(n + 1, n + n_sizes)) n += n_sizes - 1 else: - regridded_chunks.extend(c) + regridded_chunks.append(c) n += 1 @@ -3948,6 +3948,16 @@ def _regrid( min_weight=min_weight, ) + # Performance note: + # + # The function 'regrid_func' is copied into every Dask + # task. If we included the large 'weights_dst_mask' in the + # 'partial' definition then it would also be copied to every + # task, which "will start to be a pain in a few parts of the + # pipeline" definition. Instead we can pass it in via a + # keyword argument to 'map_blocks'. + # github.com/pangeo-data/pangeo/issues/334#issuecomment-403787663 + dx = dx.map_blocks( regrid_func, weights_dst_mask=weights_dst_mask, diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 28a7360df4..f87e6be374 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -115,43 +115,45 @@ * ``'bilinear'``: Deprecated alias for ``'linear'``. * ``'conservative_1st'``: First order conservative - interpolation. Preserves the area integral of the - data across the interpolation from source to - destination. It uses the proportion of the area of - the overlapping source and destination cells to - determine appropriate weights. + interpolation. Preserves the integral of the source + field across the regridding. Weight calculation is + based on the ratio of source cell area overlapped + with the corresponding destination cell area. * ``'conservative'``: Alias for ``'conservative_1st'`` * ``'conservative_2nd'``: Second-order conservative - interpolation. As with first order conservative - interpolation, preserves the area integral of the - field between source and destination using a - weighted sum, with weights based on the - proportionate area of intersection. In addition the - second-order conservative method takes the source - gradient into account, so it yields a smoother - destination field that typically better matches the - source data. - - * ``'patch'`` Patch recovery interpolation. A second - degree 2-d polynomial regridding method, which uses - a least squares algorithm to calculate the - polynomials. This method typically results in - better approximations to values and derivatives when + interpolation. Preserves the integral of the source + field across the regridding. Weight calculation is + based on the ratio of source cell area overlapped + with the corresponding destination cell area. The + second-order conservative calculation also includes + the gradient across the source cell, so in general + it gives a smoother, more accurate representation of + the source field. This is particularly true when + going from a coarse to finer grid. + + * ``'patch'`` Patch recovery interpolation. Patch + rendezvous method of taking the least squares fit of + the surrounding surface patches. This is a higher + order method that may produce interpolation weights + that may be slightly less than 0 or slightly greater + than 1. This method typically results in better + approximations to values and derivatives when compared to bilinear interpolation. - * ``'nearest_stod'``: Nearest neighbour interpolation - for which each destination point is mapped to the - closest source point. Useful for extrapolation of - categorical data. Some destination cells may be - unmapped. + * ``'nearest_stod'``: Nearest neighbour source to + destination interpolation for which each destination + point is mapped to the closest source point. A + source point can be mapped to multiple destination + points. Useful for regridding categorical data. - * ``'nearest_dtos'``: Nearest neighbour interpolation - for which each source point is mapped to the - destination point. Useful for extrapolation of - categorical data. All destination cells will be - mapped. + * ``'nearest_dtos'``: Nearest neighbour destination to + source interpolation for which each source point is + mapped to the closest destination point. A + destination point can be mapped to multiple source + points. Some destination points may not be + mapped. Useful for regridding of categorical data. * `None`: This is the default and can only be used when *dst* is a `RegridOperator`.""", @@ -493,7 +495,9 @@ The computation of the weights can be much more costly than the regridding itself, in which case reading - pre-calculated weights can improve performance.""", + pre-calculated weights can improve performance. + + Ignored if *dst* is a `RegridOperator`.""", # aggregated_units "{{aggregated_units: `str` or `None`, optional}}": """aggregated_units: `str` or `None`, optional The units of the aggregated array. Set to `None` to @@ -587,6 +591,20 @@ "{{weights auto: `bool`, optional}}": """auto: `bool`, optional If True then return `False` if weights can't be found, rather than raising an exception.""", + # ln_z + "{{ln_z: `bool` or `None`, optional}}": """ln_z: `bool` or `None`, optional + If True when *z*, *src_z* or *dst_z* are also set, + calculate the vertical component of the regridding + weights using the natural logarithm of the vertical + coordinate values. This option should be used if the + quantity being regridded varies approximately linearly + with logarithm of the vertical coordinates. If False, + then the weights are calculated using unaltered + vertical values. If `None`, the default, then an + exception is raised if any of *z*, *src_z* or *dst_z* + have also been set. + + Ignored if *dst* is a `RegridOperator`.""", # pad_width "{{pad_width: sequence of `int`, optional}}": """pad_width: sequence of `int`, optional Number of values to pad before and after the edges of @@ -603,30 +621,30 @@ True, or a tuple of both if *item* is True.""", # regrid RegridOperator "{{regrid RegridOperator}}": """* `RegridOperator`: The grid is defined by a regrid - operator that has been returned by a previous call - with the *return_operator* parameter set to True. - - Unlike the other options, for which the regrid - weights need to be calculated, the regrid operator - already contains the weights. Therefore, for cases - where multiple fields with the same source grids - need to be regridded to the same destination grid, - using a regrid operator can give performance - improvements by avoiding having to calculate the - weights for each source field. Note that for the - other types of *dst* parameter, the calculation of - the regrid weights is not a lazy operation. - - .. note:: The source grid of the regrid operator is - immediately checked for compatibility with - the grid of the source field. By default - only the computationally cheap tests are - performed (checking that the coordinate - system, cyclicity and grid shape are the - same), with the grid coordinates not being - checked. The coordinates check will be - carried out, however, if the - *check_coordinates* parameter is True.""", + operator that has been returned by a previous call + with the *return_operator* parameter set to True. + + Unlike the other options, for which the regrid weights + need to be calculated, the regrid operator already + contains the weights. Therefore, for cases where + multiple fields with the same source grids need to be + regridded to the same destination grid, using a regrid + operator can give performance improvements by avoiding + having to calculate the weights for each source + field. Note that for the other types of *dst* + parameter, the calculation of the regrid weights is + not a lazy operation. + + .. note:: The source grid of the regrid operator is + immediately checked for compatibility with + the grid of the source field. By default + only the computationally cheap tests are + performed (checking that the coordinate + system, cyclicity and grid shape are the + same), with the grid coordinates not being + checked. The coordinates check will be + carried out, however, if the + *check_coordinates* parameter is True.""", # Returns cfa_file_substitutions "{{Returns cfa_file_substitutions}}": """The CFA-netCDF file name substitutions in a dictionary whose key/value pairs are the file name parts to be diff --git a/cf/field.py b/cf/field.py index c9456f27a2..571a6c1fb1 100644 --- a/cf/field.py +++ b/cf/field.py @@ -13477,6 +13477,10 @@ def regrids( check_coordinates=False, min_weight=None, weights_file=None, + src_z=None, + dst_z=None, + z=None, + ln_z=None, verbose=None, inplace=False, i=False, @@ -13486,9 +13490,11 @@ def regrids( {{regridding overview}} - The 2-d regridding takes place on a sphere, with the grid - being defined by latitude and longitude spherical polar - coordinates. + The 2-d or 3-d regridding takes place on a sphere, with the + grid being defined by latitude and longitude spherical polar + coordinates, and any available vertical coordinates. In the + 3-d case, the regridding may be done assuming linear or log + linear weights in the vertical. **Latitude and longitude coordinates** @@ -13522,7 +13528,12 @@ def regrids( Data defined on UGRID face or node cells may be regridded to any other latitude-longitude grid, including other UGRID - meshes. + meshes and DSG feature types. + + **DSG feature types* + + Data on any latitude-longitude grid (including tripolar and + UGRID meshes) may be regridded to any DSG feature type. **Cyclicity of the X axis** @@ -13654,8 +13665,57 @@ def regrids( {{weights_file: `str` or `None`, optional}} + Ignored if *dst* is a `RegridOperator`. + .. versionadded:: 3.15.2 + src_z: optional + If `None`, the default, then the regridding is 2-d in + the latitude-longitude plane. + + If not `None` then 3-d spherical regridding is enabled + by identifying the source grid vertical coordinates + from which to derive the vertical component of the + regridding weights. The vertical coordinate construct + may be 1-d or 3-d and is defined by the unique + construct returned by ``f.coordinate(src_z)`` + + Ignored if *dst* is a `RegridOperator`. + + .. versionadded:: 3.17.0 + + dst_z: optional + If `None`, the default, then the regridding is 2-d in + the latitude-longitude plane. + + If not `None` then 3-d spherical regridding is enabled + by identifying the destination grid vertical + coordinates from which to derive the vertical + component of the regridding weights. The vertical + coordinate construct may be 1-d or 3-d. + + Ignored if *dst* is a `RegridOperator`. + + .. versionadded:: 3.17.0 + + z: optional + The *z* parameter is a convenience that may be used to + replace both *src_z* and *dst_z* when they would + contain identical values. If not `None` then 3-d + spherical regridding is enabled. See *src_z* and + *dst_z* for details. + + Ignored if *dst* is a `RegridOperator`. + + *Example:* + ``z='Z'`` is equivalent to ``src_z='Z', dst_z='Z'``. + + .. versionadded:: 3.17.0 + + {{ln_z: `bool` or `None`, optional}} + + .. versionadded:: 3.17.0 + {{verbose: `int` or `str` or `None`, optional}} .. versionadded:: 3.16.0 @@ -13749,6 +13809,10 @@ def regrids( check_coordinates=check_coordinates, min_weight=min_weight, weights_file=weights_file, + src_z=src_z, + dst_z=dst_z, + z=z, + ln_z=ln_z, inplace=inplace, ) @@ -13768,6 +13832,10 @@ def regridc( check_coordinates=False, min_weight=None, weights_file=None, + src_z=None, + dst_z=None, + z=None, + ln_z=None, inplace=False, i=False, _compute_field_mass=None, @@ -13905,6 +13973,42 @@ def regridc( .. versionadded:: 3.15.2 + src_z: optional + If not `None` then *src_z* specifies the identity of a + vertical coordinate construct of the source grid. On + its own this make no difference to the result, but it + allows the setting of *ln_z* to True. + + Ignored if *dst* is a `RegridOperator`. + + .. versionadded:: 3.17.0 + + dst_z: optional + If not `None` then *dst_z* specifies the identity of a + vertical coordinate construct of the destination + grid. On its own this make no difference to the + result, but it allows the setting of *ln_z* to True. + + Ignored if *dst* is a `RegridOperator`. + + .. versionadded:: 3.17.0 + + z: optional + The *z* parameter is a convenience that may be used to + replace both *src_z* and *dst_z* when they would + contain identical values. + + Ignored if *dst* is a `RegridOperator`. + + *Example:* + ``z='Z'`` is equivalent to ``src_z='Z', dst_z='Z'``. + + .. versionadded:: 3.17.0 + + {{ln_z: `bool` or `None`, optional}} + + .. versionadded:: 3.17.0 + {{inplace: `bool`, optional}} axis_order: sequence, optional @@ -13993,6 +14097,10 @@ def regridc( check_coordinates=check_coordinates, min_weight=min_weight, weights_file=weights_file, + src_z=src_z, + dst_z=dst_z, + z=z, + ln_z=ln_z, inplace=inplace, ) diff --git a/cf/regrid/regrid.py b/cf/regrid/regrid.py index d24682a285..4523072e56 100644 --- a/cf/regrid/regrid.py +++ b/cf/regrid/regrid.py @@ -1,6 +1,7 @@ """Worker functions for regridding.""" import logging from dataclasses import dataclass, field +from datetime import datetime from typing import Any import dask.array as da @@ -52,6 +53,12 @@ class Grid: """ + # Identify the grid as 'source' or 'destination'. + name: str = "" + # The coordinate system of the grid. + coord_sys: str = "" + # The type of the grid. E.g. 'structured grid' + type: str = "" # The domain axis identifiers of the regrid axes, in the order # expected by `Data._regrid`. E.g. ['domainaxis3', 'domainaxis2'] axis_keys: list = field(default_factory=list) @@ -66,8 +73,8 @@ class Grid: n_regrid_axes: int = 0 # The dimensionality of the regridding on this grid. Generally the # same as *n_regrid_axes*, but for a regridding a UGRID mesh axis - # *n_regrid_axes* is 1 and *regridding_dimensionality* is 2. - regridding_dimensionality: int = 0 + # *n_regrid_axes* is 1 and *dimensionality* is 2. + dimensionality: int = 0 # The shape of the regridding axes, in the same order as the # 'axis_keys' attribute. E.g. (96, 73) or (1243,) shape: tuple = None @@ -82,12 +89,8 @@ class Grid: # Only used if `mesh` is False. For spherical regridding, whether # or not the longitude axis is cyclic. cyclic: Any = None - # The coordinate system of the grid. - coord_sys: str = "" # The regridding method. method: str = "" - # Identify the grid as 'source' or 'destination'. - name: str = "" # If True then, for 1-d regridding, the esmpy weights are generated # for a 2-d grid for which one of the dimensions is a size 2 dummy # dimension. @@ -96,6 +99,10 @@ class Grid: is_grid: bool = False # Whether or not the grid is a UGRID mesh. is_mesh: bool = False + # Whether or not the grid is a location stream. + is_locstream: bool = False + # The type of grid. + type: str = "unknown" # The location on a UGRID mesh topology of the grid cells. An # empty string means that the grid is not a UGRID mesh # topology. E.g. '' or 'face'. @@ -103,6 +110,10 @@ class Grid: # A domain topology construct that spans the regrid axis. A value # of None means that the grid is not a UGRID mesh topology. domain_topology: Any = None + # The featureType of a discrete sampling geometry grid. An empty + # string means that the grid is not a DSG. E.g. '' or + # 'trajectory'. + featureType: str = "" # The domain axis identifiers of new axes that result from the # regridding operation changing the number of data dimensions # (e.g. by regridding a source UGRID (1-d) grid to a destination @@ -110,6 +121,14 @@ class Grid: # the regridding did not change the number of data axes. E.g. [], # ['domainaxis3', 'domainaxis4'], ['domainaxis4'] new_axis_keys: list = field(default_factory=list) + # Specify vertical regridding coordinates. E.g. 'air_pressure', + # 'domainaxis0' + z: Any = None + # Whether or not to use ln(z) when calculating vertical weights + ln_z: bool = False + # The integer position in *coords* of a vertical coordinate. If + # `None` then there are no vertical coordinates. + z_index: Any = None def regrid( @@ -124,6 +143,10 @@ def regrid( src_axes=None, dst_axes=None, axes=None, + src_z=None, + dst_z=None, + z=None, + ln_z=None, ignore_degenerate=True, return_operator=False, check_coordinates=False, @@ -258,6 +281,47 @@ def regrid( .. versionadded:: 3.15.2 + src_z: optional + The identity of the source grid vertical coordinates used + to calculate the weights. If `None` then no vertical axis + is identified, and in the spherical case regridding will + be 2-d. + + See `cf.Field.regrids` (for spherical regridding) or + `cf.Field.regridc` (for Cartesian regridding) for details. + + .. versionadded:: 3.17.0 + + dst_z: optional + The identity of the destination grid vertical coordinates + used to calculate the weights. If `None` then no vertical + axis is identified, and in the spherical case regridding + will be 2-d. + + See `cf.Field.regrids` (for spherical regridding) or + `cf.Field.regridc` (for Cartesian regridding) for details. + + .. versionadded:: 3.17.0 + + z: optional + The *z* parameter is a convenience that may be used to + replace both *src_z* and *dst_z* when they would contain + identical values. + + See `cf.Field.regrids` (for spherical regridding) or + `cf.Field.regridc` (for Cartesian regridding) for details. + + .. versionadded:: 3.17.0 + + ln_z: `bool` or `None`, optional + Whether or not the weights are to be calculated with the + natural logarithm of vertical coordinates. + + See `cf.Field.regrids` (for spherical regridding) or + `cf.Field.regridc` (for Cartesian regridding) for details. + + .. versionadded:: 3.17.0 + :Returns: `Field` or `None` or `RegridOperator` or `esmpy.Regrid` @@ -281,10 +345,13 @@ def regrid( # ---------------------------------------------------------------- if isinstance(dst, RegridOperator): regrid_operator = dst - dst = regrid_operator.dst + dst = regrid_operator.dst.copy() method = regrid_operator.method dst_cyclic = regrid_operator.dst_cyclic dst_axes = regrid_operator.dst_axes + dst_z = regrid_operator.dst_z + src_z = regrid_operator.src_z + ln_z = regrid_operator.ln_z if regrid_operator.coord_sys == "Cartesian": src_axes = regrid_operator.src_axes @@ -292,6 +359,31 @@ def regrid( else: create_regrid_operator = True + # Parse the z, src_z, dst_z, and ln_z parameters + if z is not None: + if dst_z is None: + dst_z = z + else: + raise ValueError("Can't set both 'z' and 'dst_z'") + + if src_z is None: + src_z = z + else: + raise ValueError("Can't set both 'z' and 'src_z'") + + elif (src_z is None) != (dst_z is None): + raise ValueError( + "Must set both 'src_z' and 'dst_z', or neither of them" + ) + + if ln_z is None and src_z is not None: + raise ValueError( + "When 'z', 'src_z', or 'dst_z' have been set, " + "'ln_z' cannot be None." + ) + + ln_z = bool(ln_z) + if method not in esmpy_methods: raise ValueError( "Can't regrid: Must set a valid regridding method from " @@ -353,7 +445,7 @@ def regrid( except TypeError: raise TypeError( "The 'dst' parameter must be one of Field, Domain, " - f"RegridOperator, or a sequence of Coordinate. Got {dst!r}" + f"RegridOperator, or a sequence of Coordinate. Got: {dst!r}" ) except IndexError: # This error will get trapped in one of the following @@ -364,42 +456,54 @@ def regrid( # destination grid into a Domain object use_dst_mask = False if spherical: - dst, dst_axes = spherical_coords_to_domain( + dst, dst_axes, dst_z = spherical_coords_to_domain( dst, dst_axes=dst_axes, cyclic=dst_cyclic, + dst_z=dst_z, domain_class=src._Domain, ) else: # Cartesian - dst, dst_axes = Cartesian_coords_to_domain( - dst, domain_class=src._Domain + dst, dst_axes, dst_z = Cartesian_coords_to_domain( + dst, dst_z=dst_z, domain_class=src._Domain ) # ---------------------------------------------------------------- # Create descriptions of the source and destination grids # ---------------------------------------------------------------- dst_grid = get_grid( - coord_sys, dst, "destination", method, dst_cyclic, axes=dst_axes + coord_sys, + dst, + "destination", + method, + dst_cyclic, + axes=dst_axes, + z=dst_z, + ln_z=ln_z, ) src_grid = get_grid( - coord_sys, src, "source", method, src_cyclic, axes=src_axes + coord_sys, + src, + "source", + method, + src_cyclic, + axes=src_axes, + z=src_z, + ln_z=ln_z, ) if is_log_level_debug(logger): logger.debug( - f"Source Grid:\n{src_grid}\nDestination Grid:\n{dst_grid}\n" + f"Source Grid:\n{src_grid}\n\nDestination Grid:\n{dst_grid}\n" ) # pragma: no cover - conform_coordinate_units(src_grid, dst_grid) + conform_coordinates(src_grid, dst_grid) if method in ("conservative_2nd", "patch"): - if not ( - src_grid.regridding_dimensionality == 2 - and dst_grid.regridding_dimensionality == 2 - ): + if not (src_grid.dimensionality >= 2 and dst_grid.dimensionality >= 2): raise ValueError( - f"{method!r} regridding is only available for 2-d regridding" + 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( @@ -410,19 +514,31 @@ def regrid( "source and destination grids have coordinate arrays" ) - if method == "nearest_dtos" and ( - src_grid.is_mesh is not dst_grid.is_mesh - ): - raise ValueError( - f"{method!r} regridding is (at the moment) only available " - "when neither or both the source and destination grids are " - "a UGRID mesh" - ) + if method == "nearest_dtos": + if src_grid.is_mesh is not dst_grid.is_mesh: + raise ValueError( + f"{method!r} regridding is (at the moment) only available " + "when neither or both the source and destination grids " + "are a UGRID mesh" + ) + + if src_grid.is_locstream or dst_grid.is_locstream: + raise ValueError( + f"{method!r} regridding is (at the moment) only available " + "when neither the source and destination grids are " + "DSG featureTypes." + ) elif cartesian and (src_grid.is_mesh or dst_grid.is_mesh): raise ValueError( - "Cartesian regridding is (at the moment) only available when " - "neither the source nor destination grid is a UGRID mesh" + "Cartesian regridding is (at the moment) not available when " + "either the source or destination grid is a UGRID mesh" + ) + + elif cartesian and (src_grid.is_locstream or dst_grid.is_locstream): + raise ValueError( + "Cartesian regridding is (at the moment) not available when " + "either the source or destination grid is a DSG featureType" ) if create_regrid_operator: @@ -480,7 +596,7 @@ def regrid( if is_log_level_debug(logger): logger.debug( - f"Source ESMF Grid:\n{src_esmpy_grid}\nDestination ESMF Grid:\n{dst_esmpy_grid}\n" + f"Source ESMF Grid:\n{src_esmpy_grid}\n\nDestination ESMF Grid:\n{dst_esmpy_grid}\n" ) # pragma: no cover esmpy_regrid_operator = [] if return_esmpy_regrid_operator else None @@ -524,6 +640,7 @@ def regrid( col, coord_sys=coord_sys, method=method, + dimensionality=src_grid.dimensionality, src_shape=src_grid.shape, dst_shape=dst_grid.shape, src_cyclic=src_grid.cyclic, @@ -538,6 +655,10 @@ def regrid( dst=dst, weights_file=weights_file if from_file else None, src_mesh_location=src_grid.mesh_location, + dst_featureType=dst_grid.featureType, + src_z=src_grid.z, + dst_z=dst_grid.z, + ln_z=ln_z, ) else: if weights_file is not None: @@ -613,7 +734,7 @@ def regrid( def spherical_coords_to_domain( - dst, dst_axes=None, cyclic=None, domain_class=None + dst, dst_axes=None, cyclic=None, dst_z=None, domain_class=None ): """Convert a sequence of `Coordinate` to spherical grid definition. @@ -643,18 +764,43 @@ def spherical_coords_to_domain( inferred from the coordinates, defaulting to `False` if it can not be determined. + dst_z: optional + If `None`, the default, then it assumed that none of the + coordinate consructs in *dst* are vertical coordinates. + + Otherwise identify the destination grid vertical + coordinate construct as the unique construct returned by + ``d.coordinate(dst_z)``, where ``d`` is the `Domain` + returned by this function. + + .. versionadded:: 3.17.0 + domain_class: `Domain` class The domain class used to create the new `Domain` instance. :Returns: - `Domain`, `dict` - The new domain containing the grid; and a dictionary - identifying the domain axis identifiers of the X and Y - regrid axes (as defined by the *dst_axes* parameter of - `cf.Field.regrids`). + 3-`tuple` + * The new domain containing the grid + * A dictionary identifying the domain axis identifiers of + the regrid axes (as defined by the *dst_axes* parameter + of `cf.Field.regrids`) + * The value of *dst_z*. Either `None`, or replaced with + its construct identifier in the output `Domain`. """ + if dst_z is None: + if len(dst) != 2: + raise ValueError( + "Expected a sequence of latitude and longitude " + f"coordinate constructs. Got: {dst!r}" + ) + elif len(dst) != 3: + raise ValueError( + "Expected a sequence of latitude, longitude, and vertical " + f"coordinate constructs. Got: {dst!r}" + ) + coords = {} for c in dst: try: @@ -665,13 +811,21 @@ def spherical_coords_to_domain( elif c.Units.islongitude: c.standard_name = "longitude" coords["lon"] = c + elif dst_z is not None: + coords["Z"] = c except AttributeError: pass - if len(coords) != 2: + if "lat" not in coords or "lon" not in coords: + raise ValueError( + "Expected a sequence that includes latitude and longitude " + f"coordinate constructs. Got: {dst!r}" + ) + + if dst_z is not None and "Z" not in coords: raise ValueError( - "When 'dst' is a sequence it must be of latitude and " - f"longitude coordinate constructs. Got: {dst!r}" + "Expected a sequence that includes vertical " + f"coordinate constructs. Got: {dst!r}" ) coords_1d = False @@ -681,31 +835,33 @@ def spherical_coords_to_domain( axis_sizes = [coords["lat"].size, coords["lon"].size] if coords["lon"].ndim != 1: raise ValueError( - "When 'dst' is a sequence of latitude and longitude " - "coordinates, they must have the same number of dimensions." + "When 'dst' is a sequence containing latitude and " + "longitude coordinate constructs, they must have the " + f"same shape. Got: {dst!r}" ) elif coords["lat"].ndim == 2: message = ( - "When 'dst' is a sequence of 2-d latitude and longitude " - "coordinates, then 'dst_axes' must be either " - "{'X': 0, 'Y': 1} or {'X': 1, 'Y': 0}" + "When 'dst' is a sequence containing 2-d latitude and longitude " + "coordinate constructs, 'dst_axes' must be dictionary with at " + "least the keys {'X': 0, 'Y': 1} or {'X': 1, 'Y': 0}. " + f"Got: {dst_axes!r}" ) if dst_axes is None: raise ValueError(message) axis_sizes = coords["lat"].shape - if dst_axes["X"] == 0: + if dst_axes.get("X") == 0 and dst_axes.get("Y") == 1: axis_sizes = axis_sizes[::-1] - elif dst_axes["Y"] != 0: + elif not (dst_axes.get("X") == 1 and dst_axes.get("Y") == 0): raise ValueError(message) if coords["lat"].shape != coords["lon"].shape: raise ValueError( "When 'dst' is a sequence of latitude and longitude " - "coordinates, they must have the same shape." + f"coordinates, they must have the same shape. Got: {dst!r}" ) - else: + elif coords["lat"].ndim > 2: raise ValueError( "When 'dst' is a sequence of latitude and longitude " "coordinates, they must be either 1-d or 2-d." @@ -719,10 +875,16 @@ def spherical_coords_to_domain( key = d.set_construct(d._DomainAxis(size), copy=False) axis_keys.append(key) + coord_axes = [] if coords_1d: # Set 1-d coordinates + coord_axes = [] for key, axis in zip(("lat", "lon"), axis_keys): - d.set_construct(coords[key], axes=axis, copy=False) + da_key = d.set_construct(coords[key], axes=axis, copy=False) + coord_axes.append(da_key) + + if dst_axes and "X" in dst_axes and dst_axes["X"] == 0: + coord_axes = coord_axes[::-1] else: # Set 2-d coordinates coord_axes = axis_keys @@ -738,10 +900,46 @@ def spherical_coords_to_domain( dst_axes = {"Y": axis_keys[0], "X": axis_keys[1]} - return d, dst_axes + if dst_z is not None: + # ------------------------------------------------------------ + # Deal with Z coordinates + # ------------------------------------------------------------ + z_coord = coords["Z"] + if z_coord.ndim == 1: + z_axis = d.set_construct(d._DomainAxis(z_coord.size), copy=False) + d.set_construct(z_coord, axes=z_axis, copy=False) + elif z_coord.ndim == 3: + if dst_axes is None or "Z" not in dst_axes or dst_axes["Z"] != 2: + raise ValueError( + "When 'dst' is a sequence containing a 3-d vertical " + "coordinate construct, 'dst_axes' must be either " + "{'X': 0, 'Y': 1, 'Z': 2} or {'X': 1, 'Y': 0, 'Z': 2}. " + f"Got: {dst_axes!r}" + ) + + z_axis = d.set_construct( + d._DomainAxis(z_coord.shape[2]), copy=False + ) + z_key = d.set_construct( + z_coord, axes=coord_axes + (z_axis,), copy=False + ) + + # Check that z_coord is indeed a vertical coordinate + # construct, and replace 'dst_z' with its construct + # identifier. + key = d.coordinate(dst_z, key=True, default=None) + if key != z_key: + raise ValueError( + f"Could not find destination {dst_z!r} vertical coordinates" + ) + + dst_z = key + dst_axes["Z"] = z_axis + return d, dst_axes, dst_z -def Cartesian_coords_to_domain(dst, domain_class=None): + +def Cartesian_coords_to_domain(dst, dst_z=None, domain_class=None): """Convert a sequence of `Coordinate` to Cartesian grid definition. .. versionadded:: 3.14.0 @@ -755,20 +953,31 @@ def Cartesian_coords_to_domain(dst, domain_class=None): field regridding axes defined elsewhere by the *src_axes* or *axes* parameter. + dst_z: optional + If `None`, the default, then it assumed that none of the + coordinate consructs in *dst* are vertical coordinates. + + Otherwise identify the destination grid vertical + coordinate construct as the unique construct returned by + ``d.coordinate(dst_z)``, where ``d`` is the `Domain` + returned by this function. + + .. versionadded:: 3.17.0 + domain_class: `Domain` class The domain class used to create the new `Domain` instance. :Returns: - `Domain`, `list` - The new domain containing the grid; and a list identifying - the domain axis identifiers of the regrid axes (as defined - by the *dst_axes* parameter of `cf.Field.regridc`). + 3-`tuple` + * The new domain containing the grid + * A list identifying the domain axis identifiers of the + regrid axes (as defined by the *dst_axes* parameter of + `cf.Field.regridc`) + * The value of *dst_z*. Either `None`, or replaced with + its construct identifier in the output `Domain`. """ - if not dst: - pass - d = domain_class() axis_keys = [] @@ -777,10 +986,30 @@ def Cartesian_coords_to_domain(dst, domain_class=None): d.set_construct(coord, axes=axis, copy=True) axis_keys.append(axis) - return d, axis_keys + if dst_z is not None: + # Check that there are vertical coordinates, and replace + # 'dst_z' with its construct identifier. + z_key = d.coordinate(dst_z, key=True, default=None) + if z_key is None: + raise ValueError( + f"Could not find destination {dst_z!r} vertical coordinates" + ) + + dst_z = z_key + + return d, axis_keys, dst_z -def get_grid(coord_sys, f, name=None, method=None, cyclic=None, axes=None): +def get_grid( + coord_sys, + f, + name=None, + method=None, + cyclic=None, + axes=None, + z=None, + ln_z=None, +): """Get axis and coordinate information for regridding. See `spherical_grid` and `Cartesian_grid` for details. @@ -790,17 +1019,42 @@ def get_grid(coord_sys, f, name=None, method=None, cyclic=None, axes=None): coord_sys: `str` The coordinate system of the source and destination grids. + ln_z: `bool` or `None`, optional + Whether or not the weights are to be calculated with the + natural logarithm of vertical coordinates. + + See `cf.Field.regrids` (for spherical regridding) or + `cf.Field.regridc` (for Cartesian regridding) for details. + + .. versionadded:: 3.17.0 + """ if coord_sys == "spherical": return spherical_grid( - f, name=name, method=method, cyclic=cyclic, axes=axes + f, + name=name, + method=method, + cyclic=cyclic, + axes=axes, + z=z, + ln_z=ln_z, ) # Cartesian - return Cartesian_grid(f, name=name, method=method, axes=axes) + return Cartesian_grid( + f, name=name, method=method, axes=axes, z=z, ln_z=ln_z + ) -def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): +def spherical_grid( + f, + name=None, + method=None, + cyclic=None, + axes=None, + z=None, + ln_z=None, +): """Get latitude and longitude coordinate information. Retrieve the latitude and longitude coordinates of a field, as @@ -841,6 +1095,25 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): ``0`` and ``1`` are axis positions in the 2-d coordinates of *f*. + z: optional + If `None`, the default, then the regridding is 2-d in + the latitude-longitude plane. + + If not `None` then 3-d spherical regridding is enabled by + identifying the grid vertical coordinates from which to + derive the vertical component of the regridding + weights. The vertical coordinate construct may be 1-d or + 3-d and is defined by the unique construct returned by + ``f.coordinate(src_z)`` + + .. versionadded:: 3.17.0 + + ln_z: `bool` or `None`, optional + Whether or not the weights are to be calculated with the + natural logarithm of vertical coordinates. + + .. versionadded:: 3.17.0 + :Returns: `Grid` @@ -852,8 +1125,12 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): dim_coords_1d = False aux_coords_2d = False aux_coords_1d = False - domain_topology = None - mesh_location = "" + domain_topology, mesh_location, axis1 = get_mesh(f) + + if not mesh_location: + featureType, axis1 = get_dsg(f) + else: + featureType = None # Look for 1-d X and Y dimension coordinates lon_key_1d, lon_1d = f.dimension_coordinate( @@ -873,6 +1150,8 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): lat = lat_1d if not dim_coords_1d: + domain_topology, mesh_location, mesh_axis = get_mesh(f) + # Look for 2-d X and Y auxiliary coordinates lon_key, lon = f.auxiliary_coordinate( "X", filter_by_naxes=(2,), item=True, default=(None, None) @@ -920,45 +1199,52 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): f"spanned by the {name} grid latitude and longitude " "coordinates" ) - else: - domain_topology, mesh_location, mesh_axis = get_mesh(f) - if mesh_location in ("face", "point"): - lon = f.auxiliary_coordinate( - "X", - filter_by_axis=(mesh_axis,), - axis_mode="exact", - default=None, - ) - lat = f.auxiliary_coordinate( - "Y", - filter_by_axis=(mesh_axis,), - axis_mode="exact", - default=None, - ) - if ( - lon is not None - and lat is not None - and lon.Units.islongitude - and lat.Units.islatitude - ): - # Found 1-d latitude and longitude auxiliary - # coordinates for a UGRID mesh topology - aux_coords_1d = True - x_axis = mesh_axis - y_axis = mesh_axis - elif mesh_location: + + elif mesh_location is not None or featureType is not None: + if mesh_location and mesh_location not in ("face", "point"): raise ValueError( f"Can't regrid {'from' if name == 'source' else 'to'} " - f"a {name} grid comprising an unstructured mesh of " + f"a {name} unstructured mesh of " f"{mesh_location!r} cells" ) + if featureType and conservative_regridding(method): + raise ValueError( + f"Can't do {method} regridding " + f"{'from' if name == 'source' else 'to'} " + f"a {name} DSG featureType" + ) + + lon = f.auxiliary_coordinate( + "X", + filter_by_axis=(axis1,), + axis_mode="exact", + default=None, + ) + lat = f.auxiliary_coordinate( + "Y", + filter_by_axis=(axis1,), + axis_mode="exact", + default=None, + ) + if ( + lon is not None + and lat is not None + and lon.Units.islongitude + and lat.Units.islatitude + ): + # Found 1-d latitude and longitude auxiliary + # coordinates for a UGRID mesh topology + aux_coords_1d = True + x_axis = axis1 + y_axis = axis1 + if not (dim_coords_1d or aux_coords_2d or aux_coords_1d): raise ValueError( "Could not find 1-d nor 2-d latitude and longitude coordinates" ) - if not mesh_location and x_axis == y_axis: + if not (mesh_location or featureType) and x_axis == y_axis: raise ValueError( "The X and Y axes must be distinct, but they are " f"the same for {name} field {f!r}." @@ -991,7 +1277,7 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): coords[dim] = coords[dim].transpose(esmpy_order) # Set cyclicity of X axis - if mesh_location: + if mesh_location or featureType: cyclic = None elif cyclic is None: cyclic = f.iscyclic(x_axis) @@ -1001,15 +1287,100 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): axes = {"X": x_axis, "Y": y_axis} axis_keys = [y_axis, x_axis] shape = (y_size, x_size) - if mesh_location: + + if mesh_location or featureType: axis_keys = axis_keys[0:1] shape = shape[0:1] n_regrid_axes = len(axis_keys) regridding_dimensionality = n_regrid_axes - if mesh_location: + if mesh_location or featureType: regridding_dimensionality += 1 + if z is not None: + # ------------------------------------------------------------ + # 3-d spherical regridding + # ------------------------------------------------------------ + if conservative_regridding(method): + raise ValueError(f"Can't do {method} 3-d spherical regridding") + + if mesh_location: + raise ValueError( + "Can't do 3-d spherical regridding " + f"{'from' if name == 'source' else 'to'} a {name} " + f"unstructured mesh of {mesh_location!r} cells" + ) + elif featureType: + # Look for 1-d Z auxiliary coordinates + z_key, z_1d = f.auxiliary_coordinate( + z, + filter_by_axis=(axis1,), + axis_mode="exact", + item=True, + default=(None, None), + ) + if z_1d is None: + raise ValueError( + f"Could not find {name} DSG {featureType} 1-d " + f"{z!r} coordinates" + ) + + z_coord = z_1d + z_axis = axis1 + else: + # Look for 1-d Z dimension coordinates + z_key, z_1d = f.dimension_coordinate( + z, item=True, default=(None, None) + ) + if z_1d is not None: + z_axis = data_axes[z_key][0] + z_coord = z_1d + else: + # Look for 3-d Z auxiliary coordinates + z_key, z_3d = f.auxiliary_coordinate( + z, + filter_by_naxes=(3,), + axis_mode="exact", + item=True, + default=(None, None), + ) + if z_3d is None: + raise ValueError( + f"Could not find {name} structured grid 1-d or 3-d " + f"{z!r} coordinates" + ) + + coord_axes = data_axes[z_key] + if x_axis not in coord_axes or y_axis not in coord_axes: + raise ValueError( + f"The {name} structured grid 3-d {z!r} " + "coordinates do not span the latitude and longitude " + "domain axes" + ) + + z_axis = [ + axis for axis in coord_axes if axis not in (x_axis, y_axis) + ][0] + + # Re-order 3-d Z coordinates to ESMF order + esmpy_order = [ + coord_axes.index(axis) for axis in (x_axis, y_axis, z_axis) + ] + z_coord = z_3d.transpose(esmpy_order) + + coords.append(z_coord) # esmpy order + + if not (mesh_location or featureType): + axes["Z"] = z_axis + axis_keys.insert(0, z_axis) + z_size = domain_axes[z_axis].size + shape = (z_size,) + shape + + regridding_dimensionality += 1 + z_index = 2 + else: + z_index = None + if f.construct_type == "domain": axis_indices = list(range(n_regrid_axes)) else: @@ -1026,28 +1397,38 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): axis_indices = [data_axes.index(key) for key in axis_keys] is_mesh = bool(mesh_location) + is_locstream = bool(featureType) + is_grid = not is_mesh and not is_locstream - return Grid( + grid = Grid( + name=name, + coord_sys="spherical", + method=method, axis_keys=axis_keys, axis_indices=axis_indices, axes=axes, n_regrid_axes=n_regrid_axes, - regridding_dimensionality=regridding_dimensionality, + dimensionality=regridding_dimensionality, shape=shape, coords=coords, bounds=get_bounds(method, coords, mesh_location), cyclic=cyclic, - coord_sys="spherical", - method=method, - name=name, + is_grid=is_grid, is_mesh=is_mesh, - is_grid=not is_mesh, + is_locstream=is_locstream, mesh_location=mesh_location, domain_topology=domain_topology, + featureType=featureType, + z=z, + ln_z=ln_z, + z_index=z_index, ) + set_grid_type(grid) + return grid -def Cartesian_grid(f, name=None, method=None, axes=None): + +def Cartesian_grid(f, name=None, method=None, axes=None, z=None, ln_z=None): """Retrieve the specified Cartesian dimension coordinates of the field and their corresponding keys. @@ -1068,6 +1449,18 @@ def Cartesian_grid(f, name=None, method=None, axes=None): Specifiers for the dimension coordinates to be retrieved. See `cf.Field.domain_axes` for details. + z: optional + If not `None` then *src_z* specifies the identity of a + vertical coordinate construct of the source grid. + + .. versionadded:: 3.17.0 + + ln_z: `bool` or `None`, optional + Whether or not the weights are to be calculated with the + natural logarithm of vertical coordinates. + + .. versionadded:: 3.17.0 + :Returns: `Grid` @@ -1092,11 +1485,17 @@ def Cartesian_grid(f, name=None, method=None, axes=None): f"for Cartesian regridding. Got: {axes!r}" ) + if z is not None and z not in axes: + raise ValueError( + f"Z coordinate {z!r} must match exactly an " + f"element of 'axes' ({axes!r})" + ) + # Find the axis keys, sizes and indices of the regrid axes, in the # order that they have been given by the 'axes' parameters. axis_keys = [] axis_sizes = [] - for axis in axes: + for i, axis in enumerate(axes): key, domain_axis = f.domain_axis(axis, item=True, default=(None, None)) if key is None: raise ValueError( @@ -1106,31 +1505,46 @@ def Cartesian_grid(f, name=None, method=None, axes=None): axis_keys.append(key) axis_sizes.append(domain_axis.size) - domain_topology, mesh_location, mesh_axis = get_mesh(f) - if mesh_location: + domain_topology, mesh_location, axis1 = get_mesh(f) + if not mesh_location: + featureType, axis1 = get_dsg(f) + else: + featureType, axis1 = None + + if mesh_location or featureType: # There is a domain topology axis - if list(set(axis_keys)) == [mesh_axis]: - # There is a unique regridding axis, and its the domain - # topology axis. - if mesh_location not in ("face", "point"): + if tuple(set(axis_keys)) == (axis1,): + # There is a unique regridding axis, and it's the discrete + # axis. + if mesh_location and mesh_location not in ("face", "point"): raise ValueError( - f"Can't regrid {'from' if name == 'source' else 'to'} " - f"a {name} grid comprising an unstructured mesh of " + f"Can't do Cartesian regridding " + f"{'from' if name == 'source' else 'to'} " + f"a {name} unstructured mesh of " f"{mesh_location!r} cells" ) + if featureType and conservative_regridding(method): + raise ValueError( + f"Can't do {method} Cartesian regridding " + f"{'from' if name == 'source' else 'to'} " + f"a {name} DSG featureType" + ) + axis_keys = axis_keys[0:1] axis_sizes = axis_sizes[0:1] - elif mesh_axis in axis_keys: + elif axis1 in axis_keys: raise ValueError( "Can't do Cartesian regridding for a combination of " - f"mesh and non-mesh axes: {axis_keys}" + f"discrete and non-discrete axes: {axis_keys}" ) else: - # None of the regridding axes have a domain topology + # None of the regridding axes have a domain topology or + # featureType domain_topology = None + featureType = None mesh_location = None - mesh_axis = None + axis1 = None if f.construct_type == "domain": axis_indices = list(range(len(axis_keys))) @@ -1148,28 +1562,26 @@ def Cartesian_grid(f, name=None, method=None, axes=None): axis_indices = [data_axes.index(key) for key in axis_keys] cyclic = False + coord_ids = axes[::-1] coords = [] - if mesh_location: - if n_axes == 1: - coord_ids = ["X", "Y"] - elif n_axes == 2: - coord_ids = axes[::-1] - else: - raise ValueError( - "Can't provide 3 or more axes for mesh axis regridding" - ) + if mesh_location or featureType: + raise ValueError( + "Cartesian regridding is (at the moment) not available when " + "either the source or destination grid is a UGRID mesh or " + "a DSG featureType" + ) + # This code may get used if we remove the above exception: for coord_id in coord_ids: aux = f.auxiliary_coordinate( coord_id, - filter_by_axis=(mesh_axis,), + filter_by_axis=(axis1,), axis_mode="exact", default=None, ) if aux is None: raise ValueError( - f"Could not find {coord_id!r} 1-d mesh auxiliary " - "coordinates" + f"Could not find {coord_id!r} 1-d auxiliary coordinates" ) coords.append(aux) @@ -1184,10 +1596,15 @@ def Cartesian_grid(f, name=None, method=None, axes=None): coords.append(dim) + if z is not None: + z_index = coord_ids.index(z) + else: + z_index = None + bounds = get_bounds(method, coords, mesh_location) dummy_size_2_dimension = False - if len(coords) == 1: + if not (mesh_location or featureType) and len(coords) == 1: # Create a dummy axis because esmpy doesn't like creating # weights for 1-d regridding data = np.array([-1.0, 1.0]) @@ -1203,33 +1620,43 @@ def Cartesian_grid(f, name=None, method=None, axes=None): n_regrid_axes = len(axis_keys) regridding_dimensionality = n_regrid_axes - if mesh_location: + if mesh_location or featureType: regridding_dimensionality += 1 is_mesh = bool(mesh_location) + is_locstream = bool(featureType) + is_grid = not is_mesh and not is_locstream - return Grid( + grid = Grid( + name=name, + coord_sys="Cartesian", + method=method, axis_keys=axis_keys, axis_indices=axis_indices, axes=axis_keys, n_regrid_axes=n_regrid_axes, - regridding_dimensionality=regridding_dimensionality, + dimensionality=regridding_dimensionality, shape=tuple(axis_sizes), coords=coords, bounds=bounds, cyclic=cyclic, - coord_sys="Cartesian", - method=method, - name=name, dummy_size_2_dimension=dummy_size_2_dimension, is_mesh=is_mesh, - is_grid=not is_mesh, + is_locstream=is_locstream, + is_grid=is_grid, mesh_location=mesh_location, domain_topology=domain_topology, + featureType=featureType, + z=z, + ln_z=ln_z, + z_index=z_index, ) + set_grid_type(grid) + return grid -def conform_coordinate_units(src_grid, dst_grid): + +def conform_coordinates(src_grid, dst_grid): """Make the source and destination coordinates have the same units. Modifies *src_grid* in-place so that its coordinates and bounds @@ -1237,7 +1664,7 @@ def conform_coordinate_units(src_grid, dst_grid): .. versionadded:: 3.14.0 - .. seealso:: `get_spherical_grid`, `regrid` + .. seealso:: `regrid` :Parameters: @@ -1252,11 +1679,6 @@ def conform_coordinate_units(src_grid, dst_grid): `None` """ - if src_grid.coord_sys == "spherical": - # For spherical coordinate systems, the units will have - # already been checked in `get_spherical_grid`. - return - for src, dst in zip( (src_grid.coords, src_grid.bounds), (dst_grid.coords, dst_grid.bounds) ): @@ -1279,6 +1701,14 @@ def conform_coordinate_units(src_grid, dst_grid): s.Units = d_units src[dim] = s + # Take the natural logarithm of spherical vertical Z coordinates + for grid in (src_grid, dst_grid): + index = grid.z_index + if grid.ln_z and index is not None: + grid.coords[index] = grid.coords[index].log() + if grid.bounds: + grid.bounds[index] = grid.bounds[index].log() + def check_operator(src, src_grid, regrid_operator, check_coordinates=False): """Check compatibility of the source field and regrid operator. @@ -1347,6 +1777,20 @@ def check_operator(src, src_grid, regrid_operator, check_coordinates=False): f"{src_grid.mesh_location} != {regrid_operator.src_mesh_location}" ) + if regrid_operator.src_featureType != src_grid.featureType: + raise ValueError( + f"Can't regrid {src!r} with {regrid_operator!r}: " + "Source grid DSG featureType mismatch: " + f"{src_grid.featureType} != {regrid_operator.src_featureType}" + ) + + if regrid_operator.dimensionality != src_grid.dimensionality: + raise ValueError( + f"Can't regrid {src!r} with {regrid_operator!r}: " + "Source grid regridding dimensionality: " + f"{src_grid.dimensionality} != {regrid_operator.dimensionality}" + ) + if not check_coordinates: return True @@ -1409,6 +1853,9 @@ def esmpy_initialise(): "patch": esmpy.RegridMethod.PATCH, } ) + reverse = {value: key for key, value in esmpy_methods.items()} + esmpy_methods.update(reverse) + # ... diverge from esmpy with respect to name for bilinear # method by using 'linear' because 'bi' implies 2D linear # interpolation, which could mislead or confuse for Cartesian @@ -1443,6 +1890,10 @@ def create_esmpy_grid(grid, mask=None): # Create an `esmpy.Mesh` return create_esmpy_mesh(grid, mask) + if grid.is_locstream: + # Create an `esmpy.LocStream` + return create_esmpy_locstream(grid, mask) + # Create an `esmpy.Grid` coords = grid.coords bounds = grid.bounds @@ -1450,37 +1901,43 @@ def create_esmpy_grid(grid, mask=None): num_peri_dims = 0 periodic_dim = 0 - spherical = False if grid.coord_sys == "spherical": spherical = True - lon, lat = 0, 1 + lon, lat, z = 0, 1, 2 coord_sys = esmpy.CoordSys.SPH_DEG if cyclic: num_peri_dims = 1 periodic_dim = lon else: # Cartesian + spherical = False coord_sys = esmpy.CoordSys.CART - # Parse coordinates for the esmpy.Grid and get its shape + # Parse coordinates for the esmpy.Grid, and get its shape. n_axes = len(coords) - coords_1d = coords[0].ndim == 1 - coords = [np.asanyarray(c) for c in coords] - if coords_1d: - # 1-d coordinates for N-d regridding - shape = [c.size for c in coords] - coords = [ - c.reshape([c.size if i == dim else 1 for i in range(n_axes)]) - for dim, c in enumerate(coords) - ] - elif n_axes == 2: - # 2-d coordinates for 2-d regridding - shape = coords[0].shape - else: - raise ValueError( - "Coordinates must be 1-d, or possibly 2-d for 2-d regridding" - ) + shape = [None] * n_axes + for dim, c in enumerate(coords[:]): + ndim = c.ndim + if ndim == 1: + # 1-d + shape[dim] = c.size + c = c.reshape([c.size if i == dim else 1 for i in range(n_axes)]) + elif ndim == 2: + # 2-d lat or lon + shape[:ndim] = c.shape + if n_axes == 3: + c = c.reshape(c.shape + (1,)) + elif ndim == 3: + # 3-d Z + shape[:ndim] = c.shape + else: + raise ValueError( + f"Can't create an esmpy.Grid from coordinates with {ndim} " + f"dimensions: {c!r}" + ) + + coords[dim] = c # Parse bounds for the esmpy.Grid if bounds: @@ -1511,22 +1968,24 @@ def create_esmpy_grid(grid, mask=None): f"{grid.method} regridding." ) - # Convert each bounds to a grid with no repeated values. - if coords_1d: - # Bounds for 1-d coordinates. - # - # E.g. if the esmpy.Grid is (X, Y) then for non-cyclic - # bounds - # we create a new bounds array with shape (97, 1); - # and for non-cyclic bounds we create a new bounds array with - # shape (1, 74). When multiplied, these arrays would - # create the 2-d (97, 74) bounds grid expected by - # esmpy.Grid. - # - # Note that if the X axis were cyclic, then its new - # bounds array would have shape (96, 1). - for dim, b in enumerate(bounds): + # Convert each bounds to a grid with no repeated values + for dim, b in enumerate(bounds[:]): + ndim = b.ndim + if ndim == 2: + # Bounds for 1-d coordinates. + # + # E.g. if the esmpy.Grid is (X, Y) then for non-cyclic + # bounds we create a new bounds array with + # shape (97, 1); and for non-cyclic bounds we + # create a new bounds array with shape (1, + # 74). When multiplied, these arrays would create + # the 2-d (97, 74) bounds grid expected by + # esmpy.Grid. + # + # Note that if the X axis were cyclic, then its + # new bounds array would have shape (96, 1). if spherical and cyclic and dim == lon: tmp = b[:, 0] else: @@ -1538,39 +1997,51 @@ def create_esmpy_grid(grid, mask=None): tmp = tmp.reshape( [tmp.size if i == dim else 1 for i in range(n_axes)] ) - bounds[dim] = tmp - else: - # Bounds for 2-d coordinates - # - # E.g. if the esmpy.Grid is (X, Y) then for bounds with a - # non-cyclic X axis, we create a new bounds array - # with shape (97, 74). - # - # Note that if the X axis were cyclic, then the new - # bounds array would have shape (96, 74). - if spherical and cyclic: - for dim, b in enumerate(bounds): - n, m = b.shape[0:2] + elif ndim == 3: + # Bounds for 2-d coordinates + # + # E.g. if the esmpy.Grid is (X, Y) then for bounds with + # a non-cyclic X axis, we create a new bounds + # array with shape (97, 74). + # + # Note that if the X axis were cyclic, then the + # new bounds array would have shape (96, 74). + n, m = b.shape[:2] + if spherical and cyclic: tmp = np.empty((n, m + 1), dtype=b.dtype) tmp[:, :m] = b[:, :, 0] if dim == lon: tmp[:, m] = b[:, -1, 0] else: tmp[:, m] = b[:, -1, 1] - - bounds[dim] = tmp - else: - for dim, b in enumerate(bounds): - n, m = b.shape[0:2] + else: tmp = np.empty((n + 1, m + 1), dtype=b.dtype) tmp[:n, :m] = b[:, :, 0] tmp[:n, m] = b[:, -1, 1] tmp[n, :m] = b[-1, :, 3] tmp[n, m] = b[-1, -1, 2] - bounds[dim] = tmp - # Define the esmpy.Grid stagger locations + if n_axes == 3: + tmp = tmp.reshape(tmp.shape + (1,)) + + elif ndim == 4: + # Bounds for 3-d coordinates + raise ValueError( + f"Can't do {grid.method} 3-d {grid.coord_sys} regridding " + f"with {grid.coord_sys} 3-d coordinates " + f"{coords[z].identity!r}." + ) + + bounds[dim] = tmp + + # Define the esmpy.Grid stagger locations. For details see + # + # 2-d: + # https://earthsystemmodeling.org/docs/release/latest/ESMF_refdoc/node5.html#fig:gridstaggerloc2d + # + # 3-d: + # https://earthsystemmodeling.org/docs/release/latest/ESMF_refdoc/node5.html#fig:gridstaggerloc3d if bounds: if n_axes == 3: staggerlocs = [ @@ -1641,6 +2112,7 @@ def create_esmpy_grid(grid, mask=None): # masked/unmasked elements. grid_mask[...] = np.invert(mask).astype("int32") + # print(esmpy_grid) return esmpy_grid @@ -1708,7 +2180,7 @@ def create_esmpy_mesh(grid, mask=None): node_owners = np.zeros(node_count) # Make sure that node IDs are >= 1, as needed by newer versions of - # esmpy + # esmpy. min_id = node_ids.min() if min_id < 1: node_ids += min_id + 1 @@ -1754,6 +2226,70 @@ def create_esmpy_mesh(grid, mask=None): return esmpy_mesh +def create_esmpy_locstream(grid, mask=None): + """Create an `esmpy.LocStream`. + + .. versionadded:: 3.17.0 + + .. seealso:: `create_esmpy_grid`, `create_esmpy_mesh` + + :Parameters: + + grid: `Grid` + The definition of the source or destination grid. + + mask: array_like, optional + The mesh mask. If `None` (the default) then there are no + masked cells, otherwise must be a 1-d Boolean array, with + True for masked elements. + + :Returns: + + `esmpy.LocStream` + The `esmpy.LocStream` derived from *grid*. + + """ + if grid.coord_sys == "spherical": + coord_sys = esmpy.CoordSys.SPH_DEG + keys = ("ESMF:Lon", "ESMF:Lat", "ESMF:Radius") + else: + # Cartesian + coord_sys = esmpy.CoordSys.CART + keys = ("ESMF:X", "ESMF:Y", "ESMF:Z") + + # Create an empty esmpy.LocStream + location_count = grid.shape[0] + esmpy_locstream = esmpy.LocStream( + location_count=location_count, + coord_sys=coord_sys, + name=grid.featureType, + ) + + # Add coordinates (must be of type float64) + for coord, key in zip(grid.coords, keys): + esmpy_locstream[key] = coord.array.astype(float) + + # Add mask (always required, and must be of type int32) + if mask is not None: + if mask.dtype != bool: + raise ValueError( + "'mask' must be None or a Boolean array. " + f"Got: dtype={mask.dtype}" + ) + + # Note: 'mask' has True/False for masked/unmasked elements, + # but the esmpy mask requires 0/1 for masked/unmasked + # elements. + mask = np.invert(mask).astype("int32") + else: + # No masked points + mask = np.full((location_count,), 1, dtype="int32") + + esmpy_locstream["ESMF:Mask"] = mask + + return esmpy_locstream + + def create_esmpy_weights( method, src_esmpy_grid, @@ -1862,20 +2398,20 @@ def create_esmpy_weights( from_file = False src_mesh_location = src_grid.mesh_location - if not src_mesh_location: - src_meshloc = None - elif src_mesh_location == "face": + if src_mesh_location == "face": src_meshloc = esmpy.api.constants.MeshLoc.ELEMENT elif src_mesh_location == "point": src_meshloc = esmpy.api.constants.MeshLoc.NODE + elif not src_mesh_location: + src_meshloc = None dst_mesh_location = dst_grid.mesh_location - if not dst_mesh_location: - dst_meshloc = None - elif dst_mesh_location == "face": + if dst_mesh_location == "face": dst_meshloc = esmpy.api.constants.MeshLoc.ELEMENT elif dst_mesh_location == "point": dst_meshloc = esmpy.api.constants.MeshLoc.NODE + elif not dst_mesh_location: + dst_meshloc = None src_esmpy_field = esmpy.Field( src_esmpy_grid, name="src", meshloc=src_meshloc @@ -1911,8 +2447,8 @@ def create_esmpy_weights( # # To do this, only retain the indices that correspond to # the top left quarter of the weights matrix in dense - # form. I.e. if w is the NxM (N, M both even) dense form - # of the weights, then this is equivalent to w[:N//2, + # form. I.e. if w is the NxM dense form of the weights (N, + # M both even), then this is equivalent to w[:N//2, # :M//2]. index = np.where( (row <= dst_esmpy_field.data.size // 2) @@ -1939,16 +2475,47 @@ def create_esmpy_weights( else: i_dtype = "i8" + upper_bounds = src_esmpy_grid.upper_bounds + if len(upper_bounds) > 1: + upper_bounds = upper_bounds[0] + + src_shape = tuple(upper_bounds) + + upper_bounds = dst_esmpy_grid.upper_bounds + if len(upper_bounds) > 1: + upper_bounds = upper_bounds[0] + + dst_shape = tuple(upper_bounds) + + regrid_method = f"{src_grid.coord_sys} {src_grid.method}" + if src_grid.ln_z: + regrid_method += f", ln {src_grid.method} in vertical" + _lock.acquire() nc = Dataset(weights_file, "w", format="NETCDF4") + nc.title = ( - f"{src_grid.coord_sys.capitalize()} {src_grid.method} " - f"regridding weights from source grid shape {src_grid.shape} " - f"to destination grid shape {dst_grid.shape}." + f"Regridding weights from source {src_grid.type} " + f"with shape {src_shape} to destination " + f"{dst_grid.type} with shape {dst_shape}" ) - nc.source = f"cf-python v{__version__}" + nc.source = f"cf v{__version__}, esmpy v{esmpy.__version__}" + nc.history = f"Created at {datetime.now()}" + nc.regrid_method = regrid_method + nc.ESMF_unmapped_action = r.unmapped_action + nc.ESMF_ignore_degenerate = int(r.ignore_degenerate) nc.createDimension("n_s", weights.size) + nc.createDimension("src_grid_rank", src_esmpy_grid.rank) + nc.createDimension("dst_grid_rank", dst_esmpy_grid.rank) + + v = nc.createVariable("src_grid_dims", i_dtype, ("src_grid_rank",)) + v.long_name = "Source grid shape" + v[...] = src_shape + + v = nc.createVariable("dst_grid_dims", i_dtype, ("dst_grid_rank",)) + v.long_name = "Destination grid shape" + v[...] = dst_shape v = nc.createVariable("S", weights.dtype, ("n_s",)) v.long_name = "Weights values" @@ -2467,6 +3034,37 @@ def get_mesh(f): ) +def get_dsg(f): + """Get domain discrete sampling geometry information. + + .. versionadded:: 3.17.0 + + :Parameters: + + f: `Field` or `Domain` + The construct from which to get the DSG information. + + :Returns: + + 2-`tuple` + If the field or domain is not a DSG then ``(None, None)`` + is returned. Otherwise the tuple contains: + + * The featureType (e.g. ``'trajectory'``) + * The identifier of domain axis construct that is spanned + by the DSG. + + """ + featureType = f.get_property("featureType", None) + if featureType is None or f.ndim != 1: + return (None, None) + + return ( + featureType, + f.get_data_axes()[0], + ) + + def has_coordinate_arrays(grid): """Whether all grid coordinates have representative arrays. @@ -2500,3 +3098,26 @@ def has_coordinate_arrays(grid): return False return True + + +def set_grid_type(grid): + """Set the ``type`` attribute of a `Grid` instance in-place. + + .. versionadded:: 3.17.0 + + :Parameters: + + grid: `Grid` + The definition of the grid. + + :Returns: + + `None` + + """ + if grid.is_grid: + grid.type = "structured grid" + elif grid.is_mesh: + grid.type = f"UGRID {grid.mesh_location} mesh" + elif grid.is_locstream: + grid.type = f"DSG {grid.featureType}" diff --git a/cf/regrid/regridoperator.py b/cf/regrid/regridoperator.py index a8f783f696..5e46b86e35 100644 --- a/cf/regrid/regridoperator.py +++ b/cf/regrid/regridoperator.py @@ -18,6 +18,8 @@ class RegridOperator(mixin_Container, Container): information, such as the grid shapes; the CF metadata for the destination grid; and the source grid coordinates. + .. versionadded:: 3.10.0 + """ def __init__( @@ -41,6 +43,13 @@ def __init__( dst=None, weights_file=None, src_mesh_location=None, + dst_mesh_location=None, + src_featureType=None, + dst_featureType=None, + dimensionality=None, + src_z=None, + dst_z=None, + ln_z=False, ): """**Initialisation** @@ -106,6 +115,11 @@ def __init__( or 1-based indexing. By default 0-based indexing is used. + If *row* and *col* are to be read from a weights file + and their netCDF variables have ``start_index`` + attributes, then these will be used in preference to + *start_index*. + parameters: Deprecated at version 3.14.0 Use keyword parameters instead. @@ -127,11 +141,56 @@ def __init__( src_mesh_location: `str`, optional The UGRID mesh element of the source grid - (e.g. ``'face'``). An empty string should be used for - a non-UGRID source grid. + (e.g. ``'face'``). .. versionadded:: 3.16.0 + dst_mesh_location: `str`, optional + The UGRID mesh element of the destination grid + (e.g. ``'face'``). + + .. versionadded:: 3.17.0 + + src_featureType: `str`, optional + The discrete sampling geometry (DSG) featureType of + the source grid (e.g. ``'trajectory'``). + + .. versionadded:: 3.17.0 + + dst_featureType: `str`, optional + The DSG featureType of the destination grid + (e.g. ``'trajectory'``). + + .. versionadded:: 3.17.0 + + src_z: optional + The identity of the source grid vertical coordinates + used to calculate the weights. If `None` then no + source grid vertical axis is identified. + + .. versionadded:: 3.17.0 + + dst_z: optional + The identity of the destination grid vertical + coordinates used to calculate the weights. If `None` + then no destination grid vertical axis is identified. + + .. versionadded:: 3.17.0 + + ln_z: `bool`, optional + Whether or not the weights were calculated with the + natural logarithm of vertical coordinates. + + .. versionadded:: 3.17.0 + + dimensionality: `int`, optional + The number of physical regridding dimensions. This may + differ from the corresponding number of storage + dimensions in the source or destination grids, if + either has an unstructured mesh or a DSG featureType. + + .. versionadded:: 3.17.0 + """ super().__init__() @@ -158,6 +217,13 @@ def __init__( self._set_component("dst", dst, copy=False) self._set_component("weights_file", weights_file, copy=False) self._set_component("src_mesh_location", src_mesh_location, copy=False) + self._set_component("dst_mesh_location", dst_mesh_location, copy=False) + self._set_component("src_featureType", src_featureType, copy=False) + self._set_component("dst_featureType", dst_featureType, copy=False) + self._set_component("dimensionality", dimensionality, copy=False) + self._set_component("src_z", src_z, copy=False) + self._set_component("dst_z", dst_z, copy=False) + self._set_component("ln_z", bool(ln_z), copy=False) def __repr__(self): """x.__repr__() <==> repr(x)""" @@ -184,6 +250,15 @@ def coord_sys(self): """ return self._get_component("coord_sys") + @property + def dimensionality(self): + """The number of physical regridding dimensions. + + .. versionadded:: 3.17.0 + + """ + return self._get_component("dimensionality") + @property def dst(self): """The definition of the destination grid. @@ -212,6 +287,15 @@ def dst_cyclic(self): cyclic.""" return self._get_component("dst_cyclic") + @property + def dst_featureType(self): + """The DSG featureType of the destination grid. + + .. versionadded:: 3.17.0 + + """ + return self._get_component("dst_featureType") + @property def dst_mask(self): """A destination grid mask to be applied to the weights matrix. @@ -229,6 +313,15 @@ def dst_mask(self): """ return self._get_component("dst_mask") + @property + def dst_mesh_location(self): + """The UGRID mesh element of the destination grid. + + .. versionadded:: 3.16.0 + + """ + return self._get_component("dst_mesh_location") + @property def dst_shape(self): """The shape of the destination grid. @@ -238,6 +331,24 @@ def dst_shape(self): """ return self._get_component("dst_shape") + @property + def dst_z(self): + """The identity of the destination grid vertical coordinates. + + .. versionadded:: 3.17.0 + + """ + return self._get_component("dst_z") + + @property + def ln_z(self): + """Whether or not vertical weights are based on ln(z). + + .. versionadded:: 3.17.0 + + """ + return self._get_component("ln_z") + @property def method(self): """The name of the regridding method. @@ -305,6 +416,15 @@ def src_cyclic(self): """ return self._get_component("src_cyclic") + @property + def src_featureType(self): + """The DSG featureType of the source grid. + + .. versionadded:: 3.17.0 + + """ + return self._get_component("src_featureType") + @property def src_mask(self): """The source grid mask that was applied during the weights @@ -340,10 +460,23 @@ def src_shape(self): """ return self._get_component("src_shape") + @property + def src_z(self): + """The identity of the source grid vertical coordinates. + + .. versionadded:: 3.17.0 + + """ + return self._get_component("src_z") + @property def start_index(self): """The start index of the row and column indices. + If `row` and `col` are to be read from a weights file and + their netCDF variables have ``start_index`` attributes, then + these will be used in preference to `start_index`. + .. versionadded:: 3.14.0 """ @@ -410,6 +543,12 @@ def copy(self): dst=self.dst.copy(), weights_file=self.weights_file, src_mesh_location=self.src_mesh_location, + dst_mesh_location=self.dst_mesh_location, + src_featureType=self.src_featureType, + dst_featureType=self.dst_featureType, + src_z=self.src_z, + dst_z=self.dst_z, + ln_z=self.ln_z, ) @_display_or_return @@ -439,6 +578,7 @@ def dump(self, display=True): for attr in ( "coord_sys", "method", + "dimensionality", "src_shape", "dst_shape", "src_cyclic", @@ -451,6 +591,12 @@ def dump(self, display=True): "src_axes", "dst_axes", "src_mesh_location", + "dst_mesh_location", + "src_featureType", + "dst_featureType", + "src_z", + "dst_z", + "ln_z", "dst", "weights", "row", @@ -550,6 +696,11 @@ def tosparse(self): any further modification of the weights to account for missing values in the source grid will always involve row-slicing. + If the weights are already in a sparse array format then no + action is taken. + + .. versionadded:: 3.13.0 + :Returns: `None` @@ -566,6 +717,10 @@ def tosparse(self): from scipy.sparse import csr_array + start_index = self.start_index + col_start_index = None + row_start_index = None + if weights is None: weights_file = self.weights_file if weights_file is not None: @@ -579,6 +734,17 @@ def tosparse(self): weights = nc.variables["S"][...] row = nc.variables["row"][...] col = nc.variables["col"][...] + + try: + col_start_index = nc.variables["col"].start_index + except AttributeError: + col_start_index = 1 + + try: + row_start_index = nc.variables["row"].start_index + except AttributeError: + row_start_index = 1 + nc.close() _lock.release() else: @@ -588,12 +754,17 @@ def tosparse(self): "be set" ) - # Convert to sprase array format - start_index = self.start_index - if start_index: - row = row - start_index + # Convert to sparse array format + if col_start_index: + col = col - col_start_index + elif start_index: col = col - start_index + if row_start_index: + row = row - row_start_index + elif start_index: + row = row - start_index + src_size = prod(self.src_shape) dst_size = prod(self.dst_shape) @@ -620,9 +791,12 @@ def tosparse(self): else: dst_mask = np.zeros((dst_size,), dtype=bool) - # Note: It is much more efficient to access 'weights.indptr' - # and 'weights.data' directly, rather than iterating - # over rows of 'weights' and using 'weights.getrow'. + # Performance note: + # + # It is much more efficient to access 'weights.indptr' and + # 'weights.data' directly, rather than iterating over rows of + # 'weights' and using 'weights.getrow'. + indptr = weights.indptr.tolist() data = weights.data for j, (i0, i1) in enumerate(zip(indptr[:-1], indptr[1:])): diff --git a/cf/test/create_test_files_2.npz b/cf/test/create_test_files_2.npz index 86adb54b23..dd3c2a1806 100644 Binary files a/cf/test/create_test_files_2.npz and b/cf/test/create_test_files_2.npz differ diff --git a/cf/test/create_test_files_2.py b/cf/test/create_test_files_2.py index 1433e6c6ab..336d0b6490 100644 --- a/cf/test/create_test_files_2.py +++ b/cf/test/create_test_files_2.py @@ -20,6 +20,15 @@ arrays = np.load(filename) +# -------------------------------------------------------------------- +# Add a new array +# -------------------------------------------------------------------- +# new_key = +# new_arrays = dict(arrays) +# new_arrays[new_key] = +# np.savez("create_test_files_2", **new_arrays) + + def _make_broken_bounds_cdl(filename): with open(filename, mode="w") as f: f.write( @@ -88,7 +97,7 @@ def _make_broken_bounds_cdl(filename): def _make_regrid_file(filename): - n = netCDF4.Dataset(filename, "w", format="NETCDF3_CLASSIC") + n = netCDF4.Dataset(filename, "w") n.Conventions = "CF-" + VN @@ -202,7 +211,7 @@ def _make_regrid_file(filename): def _make_cfa_file(filename): - n = netCDF4.Dataset(filename, "w", format="NETCDF4") + n = netCDF4.Dataset(filename, "w") n.Conventions = f"CF-{VN} CFA-0.6.2" n.comment = ( @@ -272,9 +281,132 @@ def _make_cfa_file(filename): return filename -broken_bounds_file = _make_broken_bounds_cdl("broken_bounds.cdl") +def _make_regrid_xyz_file(filename): + n = netCDF4.Dataset(filename, "w") + + n.Conventions = "CF-" + VN + + n.createDimension("time", 1) + n.createDimension("air_pressure", 6) + n.createDimension("bounds2", 2) + n.createDimension("latitude", 4) + n.createDimension("longitude", 4) + + latitude = n.createVariable("latitude", "f8", ("latitude",)) + latitude.standard_name = "latitude" + latitude.units = "degrees_north" + latitude.bounds = "latitude_bounds" + latitude[...] = [49.375, 50.625, 51.875, 53.125] + + longitude = n.createVariable("longitude", "f8", ("longitude",)) + longitude.standard_name = "longitude" + longitude.units = "degrees_east" + longitude.bounds = "longitude_bounds" + longitude[...] = [0.9375, 2.8125, 4.6875, 6.5625] + + longitude_bounds = n.createVariable( + "longitude_bounds", "f8", ("longitude", "bounds2") + ) + longitude_bounds[...] = [ + [0, 1.875], + [1.875, 3.75], + [3.75, 5.625], + [5.625, 7.5], + ] + + latitude_bounds = n.createVariable( + "latitude_bounds", "f8", ("latitude", "bounds2") + ) + latitude_bounds[...] = [ + [48.75, 50], + [50, 51.25], + [51.25, 52.5], + [52.5, 53.75], + ] + + time = n.createVariable("time", "f4", ("time",)) + time.standard_name = "time" + time.units = "days since 1860-1-1" + time.axis = "T" + time[...] = 183.041666666667 + + air_pressure = n.createVariable("air_pressure", "f4", ("air_pressure",)) + air_pressure.units = "hPa" + air_pressure.standard_name = "air_pressure" + air_pressure.axis = "Z" + air_pressure[...] = [1000, 955, 900, 845, 795, 745] + + ta = n.createVariable( + "ta", "f4", ("time", "air_pressure", "latitude", "longitude") + ) + ta.standard_name = "air_temperature" + ta.units = "K" + ta.cell_methods = "time: point" + ta[...] = arrays["ta"] + + n.close() + + return filename + + +def _make_dsg_trajectory_file(filename): + n = netCDF4.Dataset(filename, "w") + + n.Conventions = "CF-" + VN + n.featureType = "trajectory" + + n.createDimension("obs", 258) + + latitude = n.createVariable("latitude", "f4", ("obs",)) + latitude.standard_name = "latitude" + latitude.units = "degrees_north" + latitude[...] = arrays["dsg_latitude"] + + longitude = n.createVariable("longitude", "f4", ("obs",)) + longitude.standard_name = "longitude" + longitude.units = "degrees_east" + longitude[...] = arrays["dsg_longitude"] + time = n.createVariable("time", "f4", ("obs",)) + time.standard_name = "time" + time.units = "days since 1900-01-01 00:00:00" + time.axis = "T" + time[...] = arrays["dsg_time"] + + air_pressure = n.createVariable("air_pressure", "f4", ("obs",)) + air_pressure.units = "hPa" + air_pressure.standard_name = "air_pressure" + air_pressure.axis = "Z" + air_pressure[...] = arrays["dsg_air_pressure"] + + altitude = n.createVariable("altitude", "f4", ("obs",)) + altitude.units = "m" + altitude.standard_name = "altitude" + altitude[...] = arrays["dsg_altitude"] + + campaign = n.createVariable("campaign", str, ()) + campaign.cf_role = "trajectory_id" + campaign.long_name = "campaign" + campaign[...] = "FLIGHT" + + O3_TECO = n.createVariable("O3_TECO", "f8", ("obs",)) + O3_TECO.standard_name = "mole_fraction_of_ozone_in_air" + O3_TECO.units = "ppb" + O3_TECO.cell_methods = "time: point" + O3_TECO.coordinates = ( + "time altitude air_pressure latitude longitude campaign" + ) + O3_TECO[...] = arrays["dsg_O3_TECO"] + + n.close() + + return filename + + +broken_bounds_file = _make_broken_bounds_cdl("broken_bounds.cdl") regrid_file = _make_regrid_file("regrid.nc") +regrid_xyz_file = _make_regrid_xyz_file("regrid_xyz.nc") +dsg_trajectory_file = _make_dsg_trajectory_file("dsg_trajectory.nc") cfa_file = _make_cfa_file("cfa.nc") diff --git a/cf/test/test_RegridOperator.py b/cf/test/test_RegridOperator.py index b20cb0cfcf..dae6be8ce8 100644 --- a/cf/test/test_RegridOperator.py +++ b/cf/test/test_RegridOperator.py @@ -15,6 +15,7 @@ class RegridOperatorTest(unittest.TestCase): def test_RegridOperator_attributes(self): self.assertEqual(self.r.coord_sys, "spherical") self.assertEqual(self.r.method, "linear") + self.assertEqual(self.r.dimensionality, 2) self.assertEqual(self.r.start_index, 1) self.assertTrue(self.r.src_cyclic) self.assertFalse(self.r.dst_cyclic) @@ -30,7 +31,13 @@ def test_RegridOperator_attributes(self): self.assertIsNone(self.r.row) self.assertIsNone(self.r.col) self.assertIsNone(self.r.weights_file) - self.assertEqual(self.r.src_mesh_location, "") + self.assertFalse(self.r.src_mesh_location) + self.assertFalse(self.r.dst_mesh_location) + self.assertFalse(self.r.src_featureType) + self.assertFalse(self.r.dst_featureType) + self.assertIsNone(self.r.src_z) + self.assertIsNone(self.r.dst_z) + self.assertFalse(self.r.ln_z) def test_RegridOperator_copy(self): self.assertIsInstance(self.r.copy(), self.r.__class__) diff --git a/cf/test/test_formula_terms.py b/cf/test/test_formula_terms.py index 17f0d04e6d..035f85fbcb 100644 --- a/cf/test/test_formula_terms.py +++ b/cf/test/test_formula_terms.py @@ -888,6 +888,7 @@ def test_compute_vertical_coordinates(self): for standard_name in cf.formula_terms.FormulaTerms.standard_names: if standard_name == "atmosphere_hybrid_height_coordinate": continue + f, a, csn = _formula_terms(standard_name) g = f.compute_vertical_coordinates(verbose=None) diff --git a/cf/test/test_regrid.py b/cf/test/test_regrid.py index 52970a4982..fc429f2483 100644 --- a/cf/test/test_regrid.py +++ b/cf/test/test_regrid.py @@ -551,7 +551,6 @@ def test_Field_regridc_3d_field(self): """3-d Cartesian regridding with Field destination grid.""" methods = list(all_methods) methods.remove("conservative_2nd") - methods.remove("patch") dst = self.dst.copy() src0 = self.src.copy() @@ -643,7 +642,7 @@ def test_Field_regridc_3d_field(self): self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol)) # These methods aren't meant to work for 3-d regridding - for method in ("conservative_2nd", "patch"): + for method in ("conservative_2nd",): with self.assertRaises(ValueError): src.regridc(dst, method=method, axes=axes) diff --git a/cf/test/test_regrid_featureType.py b/cf/test/test_regrid_featureType.py new file mode 100644 index 0000000000..6efaa0f984 --- /dev/null +++ b/cf/test/test_regrid_featureType.py @@ -0,0 +1,271 @@ +import datetime +import faulthandler +import os +import unittest + +faulthandler.enable() # to debug seg faults and timeouts + +import numpy as np + +import cf + +# ESMF renamed its Python module to `esmpy` at ESMF version 8.4.0. Allow +# either for now for backwards compatibility. +esmpy_imported = False +try: + import esmpy + + esmpy_imported = True +except ImportError: + try: + # Take the new name to use in preference to the old one. + import ESMF as esmpy + + esmpy_imported = True + except ImportError: + pass + +disallowed_methods = ( + "conservative", + "conservative_2nd", + "nearest_dtos", +) + +methods = ( + "linear", + "nearest_stod", + "patch", +) + + +# Set numerical comparison tolerances +atol = 2e-12 +rtol = 0 + +meshloc = { + "face": esmpy.MeshLoc.ELEMENT, + "node": esmpy.MeshLoc.NODE, +} + + +def esmpy_regrid(coord_sys, method, src, dst, **kwargs): + """Helper function that regrids one dimension of Field data using + pure esmpy. + + Used to verify `cf.Field.regridc` + + :Returns: + + Regridded numpy masked array. + + """ + esmpy_regrid = cf.regrid.regrid( + coord_sys, + src, + dst, + method, + return_esmpy_regrid_operator=True, + **kwargs + ) + + src_meshloc = None + dst_meshloc = None + + domain_topology = src.domain_topology(default=None) + if domain_topology is not None: + src_meshloc = meshloc[domain_topology.get_cell()] + + domain_topology = dst.domain_topology(default=None) + if domain_topology is not None: + dst_meshloc = meshloc[domain_topology.get_cell()] + + src_field = esmpy.Field( + esmpy_regrid.srcfield.grid, meshloc=src_meshloc, name="src" + ) + dst_field = esmpy.Field( + esmpy_regrid.dstfield.grid, meshloc=dst_meshloc, name="dst" + ) + + fill_value = 1e20 + + array = np.squeeze(src.array) + if coord_sys == "spherical": + array = array.transpose() + + src_field.data[...] = np.ma.MaskedArray(array, copy=False).filled( + fill_value + ) + dst_field.data[...] = fill_value + + esmpy_regrid(src_field, dst_field, zero_region=esmpy.Region.SELECT) + + out = dst_field.data + + return np.ma.MaskedArray(out.copy(), mask=(out == fill_value)) + + +class RegridFeatureTypeTest(unittest.TestCase): + # Get the test source and destination fields + src_grid_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "regrid_xyz.nc" + ) + src_mesh_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "ugrid_global_1.nc" + ) + dst_featureType_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "dsg_trajectory.nc" + ) + src_grid = cf.read(src_grid_file)[0] + src_mesh = cf.read(src_mesh_file)[0] + dst_featureType = cf.read(dst_featureType_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/ESMF package.") + def test_Field_regrid_grid_to_featureType_3d(self): + self.assertFalse(cf.regrid_logging()) + + dst = self.dst_featureType.copy() + src = self.src_grid.copy() + + # Mask some destination grid points + dst[20:25] = cf.masked + + coord_sys = "spherical" + + for src_masked in (False, True): + if src_masked: + src = src.copy() + src[0, 2, 2, 3] = cf.masked + + # Loop over whether or not to use the destination grid + # masked points + for use_dst_mask in (False, True): + kwargs = { + "use_dst_mask": use_dst_mask, + "z": "air_pressure", + "ln_z": True, + } + for method in methods: + x = src.regrids(dst, method=method, **kwargs) + a = x.array + + y = esmpy_regrid(coord_sys, method, src, dst, **kwargs) + + self.assertEqual(y.size, a.size) + self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol)) + + if isinstance(a, np.ma.MaskedArray): + self.assertTrue((y.mask == a.mask).all()) + else: + self.assertFalse(y.mask.any()) + + @unittest.skipUnless(esmpy_imported, "Requires esmpy/ESMF package.") + def test_Field_regrid_grid_to_featureType_2d(self): + self.assertFalse(cf.regrid_logging()) + + dst = self.dst_featureType.copy() + src = self.src_grid.copy() + src = src[0, 0] + + # Mask some destination grid points + dst[20:25] = cf.masked + + coord_sys = "spherical" + + for src_masked in (False, True): + if src_masked: + src = src.copy() + src[0, 0, 2, 3] = cf.masked + + # Loop over whether or not to use the destination grid + # masked points + for use_dst_mask in (False, True): + kwargs = {"use_dst_mask": use_dst_mask} + for method in methods: + x = src.regrids(dst, method=method, **kwargs) + a = x.array + + y = esmpy_regrid(coord_sys, method, src, dst, **kwargs) + + self.assertEqual(y.size, a.size) + self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol)) + + if isinstance(a, np.ma.MaskedArray): + self.assertTrue((y.mask == a.mask).all()) + else: + self.assertFalse(y.mask.any()) + + @unittest.skipUnless(esmpy_imported, "Requires esmpy/ESMF package.") + def test_Field_regrid_mesh_to_featureType_2d(self): + self.assertFalse(cf.regrid_logging()) + + dst = self.dst_featureType.copy() + src = self.src_mesh.copy() + + # Mask some destination grid points + dst[20:25] = cf.masked + + coord_sys = "spherical" + + for src_masked in (False, True): + if src_masked: + src = src.copy() + src[20:30] = cf.masked + + # Loop over whether or not to use the destination grid + # masked points + for use_dst_mask in (False, True): + kwargs = {"use_dst_mask": use_dst_mask} + for method in methods: + x = src.regrids(dst, method=method, **kwargs) + a = x.array + + y = esmpy_regrid(coord_sys, method, src, dst, **kwargs) + + self.assertEqual(y.size, a.size) + self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol)) + + if isinstance(a, np.ma.MaskedArray): + self.assertTrue((y.mask == a.mask).all()) + else: + self.assertFalse(y.mask.any()) + + @unittest.skipUnless(esmpy_imported, "Requires esmpy/ESMF package.") + def test_Field_regrid_featureType_cartesian(self): + self.assertFalse(cf.regrid_logging()) + + # Cartesian regridding involving DSG featureTypes is not + # currently supported + src = self.src_grid + dst = self.dst_featureType + with self.assertRaises(ValueError): + src.regridc(dst, method="linear") + + src, dst = dst, src + with self.assertRaises(ValueError): + src.regridc(dst, method="linear") + + def test_Field_regrid_featureType_bad_methods(self): + dst = self.dst_featureType.copy() + src = self.src_grid.copy() + + for method in disallowed_methods: + with self.assertRaises(ValueError): + src.regrids(dst, method=method) + + +if __name__ == "__main__": + print("Run date:", datetime.datetime.now()) + cf.environment() + print("") + unittest.main(verbosity=2) diff --git a/docs/source/field_analysis.rst b/docs/source/field_analysis.rst index c3ea848470..ea588ccfae 100644 --- a/docs/source/field_analysis.rst +++ b/docs/source/field_analysis.rst @@ -1342,43 +1342,16 @@ Spherical regridding ^^^^^^^^^^^^^^^^^^^^ Regridding from and to spherical coordinate systems using the -`~cf.Field.regrids` method is only available for the 'X' and 'Y' axes -simultaneously. All other axes are unchanged. The calculation of the -regridding weights is based on areas and distances on the surface of -the sphere, rather in :ref:`Euclidean space `. +`~cf.Field.regrids` method is available for the 'X', 'Y' and (if +requested) 'Z' axes simultaneously. All other axes are unchanged. The +calculation of the regridding weights is based on areas and distances +on the surface of the sphere, rather than in :ref:`Euclidean space +`. -The following combinations of spherical source and destination domain -coordinate systems are available to the `~Field.regrids` method: - -============================== ============================== -Spherical source domain Spherical destination domain -============================== ============================== -`Latitude-longitude`_ `Latitude-longitude`_ -`Latitude-longitude`_ `Rotated latitude-longitude`_ -`Latitude-longitude`_ `Plane projection`_ -`Latitude-longitude`_ `Tripolar`_ -`Latitude-longitude`_ `UGRID mesh`_ -`Rotated latitude-longitude`_ `Latitude-longitude`_ -`Rotated latitude-longitude`_ `Rotated latitude-longitude`_ -`Rotated latitude-longitude`_ `Plane projection`_ -`Rotated latitude-longitude`_ `Tripolar`_ -`Rotated latitude-longitude`_ `UGRID mesh`_ -`Plane projection`_ `Latitude-longitude`_ -`Plane projection`_ `Rotated latitude-longitude`_ -`Plane projection`_ `Plane projection`_ -`Plane projection`_ `Tripolar`_ -`Plane projection`_ `UGRID mesh`_ -`Tripolar`_ `Latitude-longitude`_ -`Tripolar`_ `Rotated latitude-longitude`_ -`Tripolar`_ `Plane projection`_ -`Tripolar`_ `Tripolar`_ -`Tripolar`_ `UGRID mesh`_ -`UGRID mesh`_ `Latitude-longitude`_ -`UGRID mesh`_ `Rotated latitude-longitude`_ -`UGRID mesh`_ `Plane projection`_ -`UGRID mesh`_ `Tripolar`_ -`UGRID mesh`_ `UGRID mesh`_ -============================== ============================== +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. The most convenient usage is when the destination domain exists in another field construct. In this case, all you need to specify is the @@ -1541,13 +1514,11 @@ point will not participate in the regridding. Vertical regridding ^^^^^^^^^^^^^^^^^^^ -The only option for regridding along a vertical axis is to use -Cartesian regridding. However, care must be taken to ensure that the -vertical axis is transformed so that it's coordinate values vary -linearly. For example, to regrid data on one set of vertical pressure -coordinates to another set, the pressure coordinates may first be -transformed into the logarithm of pressure, and then changed back to -pressure coordinates after the regridding operation. +Vertical regridding may be incorporated into either spherical or +Cartesian regridding, but in both cases the vertical coordinates need +to be explicitly identified, and whether or not to calculate the +regridding weights according to the natural logarithm of the vertical +coordinate values stated. .. code-block:: python :caption: *Regrid a field construct from one set of pressure levels @@ -1566,19 +1537,8 @@ pressure coordinates after the regridding operation. Auxiliary coords: latitude(grid_latitude(11), grid_longitude(10)) = [[67.12, ..., 66.07]] degrees_north : longitude(grid_latitude(11), grid_longitude(10)) = [[-45.98, ..., -31.73]] degrees_east Coord references: grid_mapping_name:rotated_latitude_longitude - >>> z_p = v.construct('Z') - >>> print(z_p.array) - [850. 700. 500. 250. 50.] - >>> z_ln_p = z_p.log() - >>> z_ln_p.axis = 'Z' - >>> print(z_ln_p.array) - [6.74523635 6.55108034 6.2146081 5.52146092 3.91202301] - >>> _ = v.replace_construct('Z', new=z_ln_p) - >>> new_z_p = cf.DimensionCoordinate(data=cf.Data([800, 705, 632, 510, 320.], 'hPa')) - >>> new_z_ln_p = new_z_p.log() - >>> new_z_ln_p.axis = 'Z' - >>> new_v = v.regridc([new_z_ln_p], axes='Z', method='linear') - >>> new_v.replace_construct('Z', new=new_z_p) + >>> new_z = cf.DimensionCoordinate(data=cf.Data([800, 705, 632, 510, 320.], 'hPa')) + >>> new_v = v.regridc([new_z], axes='Z', method='linear', z='Z', ln_z=True) >>> print(new_v) Field: eastward_wind (ncvar%ua) ------------------------------- @@ -1592,11 +1552,6 @@ pressure coordinates after the regridding operation. : longitude(grid_latitude(11), grid_longitude(10)) = [[-45.98, ..., -31.73]] degrees_east Coord references: grid_mapping_name:rotated_latitude_longitude -Note that the `~Field.replace_construct` method of the field construct -is used to directly replace the vertical dimension coordinate construct, -without having to manually match up the corresponding domain axis -construct and construct key. - ---- .. _Mathematical-operations: @@ -2549,3 +2504,5 @@ See `cf.curl_xy` for details and examples. .. _Latitude-longitude: http://cfconventions.org/cf-conventions/cf-conventions.html#_latitude_longitude .. _Rotated latitude-longitude: http://cfconventions.org/cf-conventions/cf-conventions.html#_rotated_pole .. _Plane projection: http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings +.. _UGRID mesh: https://cfconventions.org/cf-conventions/cf-conventions.html#mesh-topology-variables +.. _DSG feature type: https://cfconventions.org/cf-conventions/cf-conventions.html#discrete-sampling-geometries