  Interpolation of data represented in curvilinear coordinates and possibly
having some non-trivial, typically Riemannian or semi-Riemannian geometry is an
ubiquitous task in all of physics. In this work we present a covariant
generalization of the barycentric coordinates and the barycentric interpolation
method for Riemannian and semi-Riemannian spaces of arbitrary dimension. We
show that our new method preserves the linear accuracy property of barycentric
interpolation in a coordinate-invariant sense. In addition, we show how the
method can be used to interpolate constrained quantities so that the given
constraint is automatically respected. We showcase the method with two
astrophysics related examples situated in the curved Kerr spacetime. The first
problem is interpolating a locally constant vector field, in which case
curvature effects are expected to be maximally important. The second example is
a General Relativistic Magnetohydrodynamics simulation of a turbulent accretion
flow around a black hole, wherein high intrinsic variability is expected to be
at least as important as curvature effects.
