Feature Proposal: Postprocessor Class for Mesh-to-Mesh Interpolation

[vieirasaulo]

🚀 Feature Proposal: Postprocessor Class for Mesh-to-Mesh Interpolation

Context

It's common to work with multiple finite element meshes that represent the same or similar spatial domains but differ in resolution, refinement, or topology. A typical use case includes:

• Transferring simulation results from one mesh (e.g., coarse mesh from a solver) to another (e.g., refined mesh for visualization or post-processing).
• Performing mesh adaptivity or remeshing and needing to preserve solution fields.
• Comparing results across simulation runs with differing discretizations.

Currently, interpolating scalar/vector fields (both nodal and elemental) between two unstructured 2D triangular meshes requires custom logic involving:

1. Delaunay triangulation and barycentric coordinate evaluation for node interpolation (scipy.spatial.Delaunay)
2. Nearest-neighbor search (scipy.spatial.cKDTree).

These are non-trivial operations that can benefit from a standardized, reusable interface that in my opinion fits very well pandamesh.

---

Feature

I propose the introduction of a Postprocessor utility class that provides mesh-to-mesh interpolation capabilities. This would serve as a general-purpose post-processing tool to:

• Interpolate nodal attributes from one mesh onto another using barycentric interpolation.
• Interpolate elemental attributes using centroid-based nearest-neighbor mapping.
• Support multi-dimensional attributes (e.g., scalar fields, vector fields).
• Handle invalid or out-of-bounds queries gracefully.

---

Proposed API (Sketch)

class Postprocessor:
    def __init__(
        self,
        source_nodes: np.ndarray,
        source_elements: np.ndarray,
    ):
        ...

    def interpolate_to(
        self,
        nodal_attributes: np.ndarray,
        elemental_attributes: np.ndarray,
        target_nodes: np.ndarray,
        target_elements: np.ndarray,
    ) -> Tuple[np.ndarray, np.ndarray]:
        """
        Interpolate nodal and elemental attributes to a new mesh.
        Returns:
            interpolated_nodal_attributes: shape (N2, Dn)
            interpolated_elemental_attributes: shape (M2, De)
        """
        ...

I’ve already prototyped something and would be happy to contribute a pull request including:

• The Postprocessor class with interpolate_to method. This is what I can see for now, but likely other use cases might emerge.
• Unit tests and examples using basic geometries
• Clear docstrings

Please let me know if this aligns with your roadmap or if you'd like to discuss further before I submit the PR!

Thanks for the great library 🙌

[Huite]

Hey @vieirasaulo!

I fully agree that this is useful functionality. In my view, this has some relation to mesh generation, but more so to general handling and processing of unstructured grid data. I've placed most of the functionality of that in xugrid, under the name of "regridding". The essential idea of xugrid is to join data without spatial axes (like structured grids) in xarray form to the unstructured grid topology.

I've written two sections on regridding in the user guide:

• Overview
• Overlap regridding

For the overview here, there are three basic regridders:

• The centroid locator just does look up via the centroid of the destination face (can be triangle, can be any convex polygon)
• The barycentric interpolator does barycentric interpolation based on the centroids of the destination face through its (Voronoi) dual; xugrid can derive this relatively efficiently from the connectivity data https://deltares.github.io/xugrid/api/xugrid.Ugrid2d.voronoi_topology.html
• The overlap regridder computes overlap weights (so areas in 2D); this can be used for area weighted regridding. It also supports custom reduction (e.g. mean, sum, etc.).

They all support regridding from and to both structured and unstructured topologies (where structured topologies are plain xarray DataArrays with an x and a y axis). Xugrid also contains some tooling for projecting vector data like geodataframes onto an unstructured topology.

Unless you have specific objections, it probably makes most sense to continue the discussion in a xugrid issue. But I'll ramble on here a bit anyway...

For pandamesh, I should perhaps refer to xugrid more explicitly. Potentially, we should add a user guide which shows one of the workflows you mention, e.g.:

1. Generate a grid with pandamesh
2. Read in some example data (e.g. elevation raster)
3. Regrid to the generated grid

Faces and nodes

I should note that so far, all these regridders currently assume face associated data (so the triangle instead of the nodes). This is because my colleagues and I mostly work with finite volume / control volume finite difference models. Supporting node associated data should be straightforward, although it might require a separate class here and there?

I notice now that my API isn't fully future proof, take the barycentric interpolator -- you can provide the target topology as either a UgridDataArray or a Ugrid2d topology. (The reason I support a UgridDataArray is because you often have a template with data available.) In first case, we could check the dimension on the data and see whether it's a node or face dimension, but not so with a bare topology. But adding an argument seems reasonable. In general, working with the vertices instead of the faces is simpler.

The xugrid base regridder class also uses xarray apply_ufunc and such to generalize the regridding to orthogonal dimensions in a lazy evaluation manner.

Thinking out loud

• The overlap regridder only makes sense for face data (a vertex will always have an area of zero).
• The centroid locator only makes sense for face data as well -- a vertex is directly the point to search.
• The BarycentricInterpolator could accept an additional argument regarding the target facet of the topology (the source array determines whether it's interpolating between the voronoi tesselation of the centroids of face associated data or straight from the vertices). This logic is already provided by numba_celltree via compute_barycentric_weights; this generalizes to any convex polygon through Wachspress barycentric coordinates (barycentric coordinates are uniquely determined for triangles only, we need additional constraints for N-gons). CGAL has a nice overview
• I've thought about a NearestRegridder as well. This would work for both node and face associated data (again by taking the centroid as representative). The source array dimension again determines which source topology facet we need to consider, and we'd need a similar argument for the destination. I think a scipy kDTree is indeed the best spatial index. A minor consideration is how to handle vertices that fall outside of the source topology. We cannot rely on a single maximum distance; I think it's easiest to run the Celltree2d.locate_points method and discard anything that's outside (maybe configurable).

numba_celltree is primarily a spatial index, and a bit computational geometry kitchen sink.

Concluding

I think it makes most sense to have this in xugrid. It's something that should be added there anyway, and it'll be able to re-use the lazy evaluation / broadcasting machinery here: https://github.com/Deltares/xugrid/blob/38da2d7e71bd87550105cfa70805e2d51dee8649/xugrid/regrid/regridder.py#L95

All the building blocks are already available in efficient form: scipy kDTree and numba_celltree for barycentric weights and the inside/outside predicate.

The main challenge is that if you wanted to add a PR, there's a sizeable amount of code to familiarize yourself with. The changes are mostly straightforward (just use the vertices instead of face centroids), but it'll require changes in many places. My suggestion is:

• Take a look at xugrid. If you're not familiar with xarray, it might have a somewhat steep learning curve (but I'd say it's very worthwhile).
• If you still have a strong preference for adding to this pandamesh instead, I'd be curious to hear it. You could argue dependencies, but we'd want numba_celltree anyway for efficiency, and then xugrid will also work with quite a minimal set of dependencies.
• I'll make two draft PRs for xugrid anyway: first for adding support for vertex associated data in the BarycentricInterpolator, and a second one for a NearestRegridder. This'll allow me to think about some aspect more concretely (e.g. what are the interactions with structured data -- my guess is few, since a structured format (think GDAL) generally doesn't allow associating data with the (implicit) vertices). I'm not sure yet about the time frame -- I won't be working much the coming weeks due to a medical procedure for my 6-month old daughter.

[vieirasaulo]

Hi @Huite,

Thank you for the very thorough response and the detailed overview of the (definitely not basic!) regridders you've implemented. It’s clear a lot of thoughtful design has gone into those tools. I agree—it would make a lot of sense to incorporate this kind of interpolation functionality into xugrid.

Before moving further in that direction, though, I wanted to get a sense of how stable you consider the xugrid API at this stage, given that the project is still relatively young.

But before shifting the conversation entirely to xugrid, I’d love to touch on a few more technical aspects—especially as I’m still relatively new to xugrid, so apologies if I’ve misunderstood anything.

---

"For pandamesh, I should perhaps refer to xugrid more explicitly. Potentially, we should add a user guide which shows one of the workflows you mention..."

That sounds like a great idea. Regarding point 1—would it make sense to construct an xugrid.Ugrid2d directly from pandamesh output? If so, and if you’re open to it, I wonder if pandamesh could be listed as an optional dependency within xugrid, enabling tighter integration. Otherwise, a minimal bridge using numpy.ndarray would keep the dependency footprint small on both ends.

---

As for point 2: does xugrid currently support nodal data directly, or does it primarily operate on face-centered fields? From what I gather, the xarray integration should allow selecting which dimension to operate on (e.g., faces vs. nodes), but I’m not entirely clear on whether nodal attributes are fully first-class citizens in the current API.

In my own (far more limited) implementation, I’ve been using barycentric interpolation to transfer nodal values between meshes. The approach works well when transferring from a finer mesh to a coarser one (and vice versa, though that’s inherently more lossy). For elemental values, I’m relying on nearest-neighbor mapping based on triangle centroids.

Here’s a distilled example of the method I’ve been using:

def interpolate_to(
    mesh1_nodes: np.ndarray,
    mesh1_elements: np.ndarray,
    mesh2_nodes: np.ndarray,
    mesh2_elements: np.ndarray,
    mesh1_nodal_attributes: np.ndarray,
    mesh1_elemental_attributes: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Interpolates nodal and elemental attributes from one triangle mesh to another.

    Parameters
    ----------
    mesh1_nodes : (N1, 2) array
        Coordinates of nodes in the source mesh.
    mesh1_elements : (T1, 3) array
        Triangle connectivity for the source mesh.
    mesh2_nodes : (N2, 2) array
        Coordinates of nodes in the target mesh.
    mesh2_elements : (T2, 3) array
        Triangle connectivity for the target mesh.
    mesh1_nodal_attributes : (N1, D) array
        Nodal data on the source mesh (D = number of fields).
    mesh1_elemental_attributes : (T1, E) array
        Elemental data on the source mesh (E = number of fields).

    Returns
    -------
    interpolated_nodal_attributes : (N2, D) array
        Interpolated nodal values on the target mesh.
    interpolated_elemental_attributes : (T2, E) array
        Interpolated elemental values on the target mesh.

    Notes
    -----
    - Nodal interpolation uses barycentric coordinates based on Delaunay triangulation.
    - Elemental interpolation uses nearest-neighbor search in centroid space.
    """
    tri = Delaunay(mesh1_nodes)
    simplices = tri.find_simplex(mesh2_nodes)
    transform = tri.transform[simplices]
    delta = mesh2_nodes - transform[:, 2]
    bary = np.einsum('ijk,ik->ij', transform[:, :2], delta)
    bary_coords = np.c_[bary, 1 - bary.sum(axis=1)]

    vertices = tri.simplices[simplices]
    gathered = mesh1_nodal_attributes[vertices]  # (N2, 3, D)
    interpolated_nodal_attributes = np.einsum('nij,ni->nj', gathered, bary_coords)

    mesh2_centroids = np.mean(mesh2_nodes[mesh2_elements], axis=1)
    mesh1_centroids = np.mean(mesh1_nodes[mesh1_elements], axis=1)
    _, idx = cKDTree(mesh1_centroids).query(mesh2_centroids)
    interpolated_elemental_attributes = mesh1_elemental_attributes[idx]

    return interpolated_nodal_attributes, interpolated_elemental_attributes

---

I could imagine a high-level xugrid.regrid module that wraps functionality like this, with options for regridding either nodal or face data and supporting method choices (e.g., 'barycentric', 'nearest', 'idw', etc.).

Let me know what you think!

Thanks again for your time :)

[Huite]

API

Before moving further in that direction, though, I wanted to get a sense of how stable you consider the xugrid API at this stage, given that the project is still relatively young.

I'd say it's fairly stable -- certainly not less stable than pandamesh! But it's obviously bigger, since pandamesh primarily just wraps existing API's; while xugrid covers IO via the UGRID conventions, plotting, spatial selection, conversion from structured data / vector data, regridding, partioning & merging, reprojection, as well as some geometric and connectivity utilities.

Pandamesh --> Ugrid2d

That sounds like a great idea. Regarding point 1—would it make sense to construct an xugrid.Ugrid2d directly from pandamesh output? If so, and if you’re open to it, I wonder if pandamesh could be listed as an optional dependency within xugrid, enabling tighter integration. Otherwise, a minimal bridge using numpy.ndarray would keep the dependency footprint small on both ends.

Yes, this already exists! See: https://deltares.github.io/pandamesh/api/api/pandamesh.TriangleMesher.generate_ugrid.html

The relevant dependency import can be found here:

pandamesh/pandamesh/common.py

Line 531 in 76e1c89

def to_ugrid(vertices: FloatArray, faces: IntArray) -> "xugrid.Ugrid2d": # type: ignore # noqa pragma: no cover

It imports xugrid inside the function so it's an optional dependency for pandamesh.

Nodal data

(Apologies for the info dump here)

As for point 2: does xugrid currently support nodal data directly, or does it primarily operate on face-centered fields? From what I gather, the xarray integration should allow selecting which dimension to operate on (e.g., faces vs. nodes), but I’m not entirely clear on whether nodal attributes are fully first-class citizens in the current API.

Nodal data is directly supported yes, e.g. plotting: https://deltares.github.io/xugrid/examples/plotting.html#other-types-of-plot
One of the sample dataset is node-based: https://deltares.github.io/xugrid/sample_data/adh_san_diego.html#sphx-glr-sample-data-adh-san-diego-py

The way it works is that xugrid will match the dimension name of the xarray DataArray with the node, edge, or face dimension of the topology. So to initialize a (dummy) UgridDataArray with nodal data, you'd do:

node_uda = xugrid.UgridDataArray(
     xr.DataArray(data=np.ones(grid.n_node), dims=(grid.node_dimension,),
     grid,
)

Vs face data:

face_uda = xugrid.UgridDataArray(
     xr.DataArray(data=np.ones(grid.n_face), dims=(grid.face_dimension,),
     grid,
)

Spatial selection/subsetting works as well, although not all operations work/are obvious when dealing with nodal data. E.g. a cross-section:

da = uda.ugrid.sel(y=section_y)

is straightforward for face centered data: we get line segments, one for each intersected face; the value of the line segment is associated with the intersected face. For nodal data, we could project perpendicularly to the line (this might gives conflicting values though), try to interpolate at discrete sampling points, etc. It's less straightforward. Some methods only work on face associated data whereas they have natural extensions/generalizations to nodal data as well. E.g. interpolation will work on nodes and faces since it generalizes naturally (and we needed it at some point): https://deltares.github.io/xugrid/examples/connectivity.html#laplace-interpolation

At any rate, there are also ways to go from node to face or vice versa:

• The voronoi tesselation will make faces nodes (using the face centroids), or will make nodes faces: https://deltares.github.io/xugrid/examples/connectivity.html#centroidal-voronoi-tesselation
  (I.e. the dual, see: https://en.wikipedia.org/wiki/Mesh_generation#Duality)
• I've recently added to_face, to_edge, and to_node methods: https://deltares.github.io/xugrid/api/xugrid.UgridDataArrayAccessor.to_face.html; these are "lossy" though, since they aggregate (for example you might use it to give a triangle face the mean value of its nodes)

Regrid API

You might wonder why xugrid has Regridding classes rather than regrid functions. The reason is that computing regridding indices and weights require quite a lot of work. We have to construct a spatial index (generally some tree) and traverse. In case of overlaps, we need to run the intersection algorithm. And in many cases, we have datasets with similar data. E.g. model input is generally defined on the same topology, so if we'd like to move to another topology, it's useful to be able to apply the regridding weights across all variables without having to recompute the costly weights. (We can even store them in files to have them persist across Python sessions.)

We could wrap them in a single regrid function, but you can already write it as a one-liner currently:

new = xugrid.BarycentricInterpolator(source=old, target=new_mesh).regrid(old)

Admittedly, this is cleaner:

new = xugrid.regrid(source=old, target=new_mesh, method="barycentric")

But it doesn't allow re-use, and is barely shorter.

Example

The current BarycentricInterpolator constructor looks like:

 def __init__(
        self,
        source: UgridDataArray,
        target: UgridDataArray | UgridTopology,
        tolerance: Optional[float] = None,
    ):
     ...

And it would become something like:

 def __init__(
        self,
        source: UgridDataArray,
        target: UgridDataArray | UgridTopology,
        target_dim: Optional[str]=None,
        tolerance: Optional[float] = None,
    ):
    if isinstance(target_dim, UgridDataArray ) and facet is not None:
         raise ValueError("target_dim should only be provided for a UgridTopology target, not UgridDataArray")
    ...

Then, to use it in your pandamesh workflow, you'd do:

# Generate mesh
mesher = TriangleMesher(geoms)
grid = mesher.generate_ugrid()

# Get some data
da = xr.open_dataarray("my-elevation-raster.nc")

# Regrid
regridder = BarycentricInterpolator(source=da, target=grid, target_dim=grid.node_dimension)
result = regridder.regrid(da)

I think most xugrid users have face associated data, but obviously xugrid should fully support nodal data as well! Feel free to open issues or a discussion in case you run into missing functionality. (I'd call the current lack of regridding support for nodal data missing functionality.)

[vieirasaulo]

Thanks again for your detailed reply and the walkthrough of the new methods added to xugrid — very helpful!

Would it be appropriate for me to open an issue proposing nodal regridding using Barycentric interpolation as a built-in feature in xugrid?

Before doing so, I wanted to better understand how regridding of both nodal and face data could fit into the current framework. Right now, it seems that a practical short-term path might involve defining two separate regridder classes. However, it feels natural to abstract this further.

A functional approach for a user-friendly API could look like:

new = xugrid.regrid(
    source=old,
    target=new_mesh,
    nodal_method="barycentric",
    face_method="nearest",
)

Under the hood, nodal_method and face_method could accept Literal-based options (e.g., "nearest", "barycentric", "linear"), with method-specific arguments handled via **kwargs. Each method would lazily regrid the respective components and return a UgridDataset with all nodal distributions and face distributions regridded.

---

Motivation: Large Lazy Workflows

As a user, one of the most exciting applications of this would be to lazily transfer data from a large dataset—say, 200 distributions over 2M nodes and 200 over 3M faces—from an old mesh to a new one (e.g. from pandamesh). Previously, I’d imagined needing to:

1. Load each nodal/face attribute separately
2. Regrid manually
3. Save back to disk one-by-one

Now, I wonder if this could be streamlined using a single regridding interface that understands the mesh topology and handles nodal and face variables transparently. Example of very useful use case below:

import xarray as xr
import xugrid as xu

src = xu.open_dataset("input.nc", chunks={})  # Lazy load ???

target_grid = TriangleMesher(geoms).generate_ugrid()

regridder = xu.Regridder(
    src,
    target_grid,
    nodal_method="barycentric",
    face_method="nearest"
)

(
    regridder
    .regrid()  # could it return a lazily regridded dataset?
    .to_netcdf("regridded_output.nc", compute=True)
)

This would let users regrid any dataset—regardless of how many nodal/face variables it contains—without needing to manually loop through or distinguish them ahead of time. A unified Regridder object could lazily map and dispatch operations based on variable topology.

---

Next Steps

I'd be happy to open a feature request/issue describing this idea if you're open to it. While I might not be ready to submit a PR myself right away, I'd be glad to help refine and document the design if it aligns with your roadmap.

Thanks again for all your work on xugrid and being open for discussion! I wish your daughter all the best and good luck with the procedure.

Best,
Saulo

[Huite]

Would it be appropriate for me to open an issue proposing nodal regridding using Barycentric interpolation as a built-in feature in xugrid?

Yes!

And as mentioned this is an interesting but clearly scoped problem, so I'm curious to experiment a little myself in a draft PR!

Right now, it seems that a practical short-term path might involve defining two separate regridder classes. However, it feels natural to abstract this further.

A NearestRegridder is a new class yes, not sure about BarycentricInterpolator. As a general tenet, "types are your friend" and there's something to be said for having a unique type for each behavior (so node is different than face), but one issue that I foresee is that there are already four combinations when considering node and face both as sources and as targets. This may grow if we add full 3D support for xugrid. On the whole, we can extend the current BarycentricInterpolator and provide something intuitive to a user. Especially when they have a source and target arrays, the regridder will automatically infer the correct behavior based on the presence of the node or face dimensions.

As a user, one of the most exciting applications of this would be to lazily transfer data from a large dataset—say, 200 distributions over 2M nodes and 200 over 3M faces—from an old mesh to a new one (e.g. from pandamesh)

Yes, this type of operation is exactly what the regridder should enable! The combination of precomputing weights, storing them, and streaming the data via xarray (dask) is exactly what should make this convenient and efficient.

I've made an issue about dataset regridding some time ago, but haven't had much reason to prioritize it yet:
Deltares/xugrid#88

A regridder class contains a single set of weights so it won't be able to regrid any dataset, like those containing multiple topologies, or both face and node associated data. We can either solve this with a xugrid.regrid function that creates the necessary regridders internally or some form of composite regridder.

The API needs some thought, as it needs flexibility:

• Both node and face data
• As source and target
• For multiple topologies (note that xugrid also contains "1D" network topologies which may be present in a dataset)
• Extended in due time for new types (fully unstructured 3D), as well as layered 2D

The most general form (which allows full control) is to allow specifying the method per variable.

regridded = xugrid.regrid(
     source=old,
     target=template,
     var_methods={"a": BarycentricInterpolator, "b": NearestRegridder},
)

In case of one regridder, we allow var_methods=BarycentricInterpolator next to the dict mapping.

This gets a bit tedious if you just want to use the same method across variables.
Hence, per dimension would allow us to make a distinction between both node and face data, but also between different topologies (since multiple topologies in a single dataset require unique dimension names):

regridded = xugrid.regrid(
     source=old,
     target=template,
     dim_methods={"mesh2d_node": BarycentricInterpolator, "grid2d_face": NearestRegridder},
)

This is easy to extend, since variables and dimensions generalize to any type of topology.

Instead of literals, we could use the classes directly. Only the OverlapRegridder takes an additional method argument (for the reduction), but that's easily addressed with a factory function:

class OverlapRegridder:
    def __init__(self, source_grid, target_grid, method="mean"):
        self.source_grid = source_grid
        self.target_grid = target_grid
        self.method = method
        # Rest of initialization
    
    @classmethod
    def from_method(cls, method: Union[Callable, str]):
        """
        Returns a factory function that creates an OverlapRegridder with the specified method.
        
        Parameters
        ----------
        method : callable or str
            The method to use for regridding.
            
        Returns
        -------
        callable
            A function that takes source_grid and target_grid and returns an OverlapRegridder.
        """
        def factory(source_grid, target_grid):
            return cls(source_grid, target_grid, method=method)
        
        return factory

In use:

regridded = xugrid.regrid(
    source=old, 
    target=template, 
    var_methods={
        "a": BarycentricInterpolator,
        "b": NearestRegridder,
        "c": OverlapRegridder.from_method("mean")
    }
)

Internally, we are free to make as many regridders as needed. I think I prefer this form over some sort of a DatasetRegridder or CompositeRegridder because it'll get a messy (internal) state. Furthermore, such a class could always be added afterwards, while maintaining the functional form with minimal repetition or overhead.

[Huite]

Given the detail of the dicussion, we should probably move this issue wholesale to the xugrid issues 😅