In the context of interpolation, a distance metric is an underlying assumption of an interpolation method, as it is implicit in the calculation of the interpolating weights. Interpolating over a spherical surface (dimensionality 2) and for a generic space unit  the weights are calculated in a \( [L]^2 \)  unit, to which corresponds a distance \( [L] \) . Because the metric is implicit to the method implementation, the metric discussion is here in terms of the specific methods.

Classical bi-linear method

The classical bi-linear method (EMOSLIB, MIR) interpolates using spherical coordinates, with curvilinear quadrilaterals as the supporting elements; calculating the interpolating weights (and distance) is fairly equivalent to (weights are de-normalized)

  • \( w_{bi} \propto (\Delta lon)^2 + (\Delta lat)^2 \)

  • \( d_{bi} = \sqrt{ w_{bi} } \)

For input regular grids (regular_ll and regular_ll) this operates over a spherical surface tessellation — ie. weights are not overlapping for different grid elements. However, for reduced grids (reduced_gg, or quasi-regular) this as an approximation because the component \( \Delta^2 lon \) assumes the supporting elements are meridian-aligned. In practice, when the number of points per latitude changes (for a reduced_gg octahedral grid every parallel has a different number of points) calculating \( w_{bi} \)  will not exactly match the supporting quadrilateral element, effectively blocking a possible local conservativity property.

Additionally, because the real distance (on the spherical surface) in longitude depends on the latitude parallel –the real distance at specific latitudes decreases approaching the poles– the weights are always distorted pole-wards away from the equator. This effect always present, however minimized due to the smaller latitude increments towards the poles of the Gaussian grids.

Triangular linear method

The triangular linear interpolation method (MIR-only) uses the 3D Euclidean space, tessellating a spherical surface with planar triangles. Weights are calculated as (note similarity to Euclidean distance):

  • \( w_{3D} \propto (\Delta x)^2 + (\Delta y)^2 + (\Delta z)^2 \)

  • \( d_{3D} = \sqrt{ w_{3D} } \)

The supporting elements form a spherical surface tessellation irrespective of the input grid format; Because the elements are described in 3D there is no special concern for the poles (the poles singularity is effectively avoided by operating in the dimension above). Using planar supporting elements does however introduce an approximation compared to the (real) spherical surface distance, an error minimized with resolution increase.

2 Comments

  1. Tim Stockdale

    Am I correct that the weights used in the interpolation are linear functions of the distance d, and are entirely different from the w described above, which are used to calculate d but are not directly or inversely used as weights? A slightly fuller description or reference might be helpful: as I understand it, the triangular linear method has a very simple geometrical interpretation.

    1. Pedro Maciel

      The description of the distance metrics is only to explain the MIR/EMOSLIB geometry differences. The practical difference is, as you say, in the calculation of the interpolation weigths. I only intended to show that there is a fundamental difference between MIR/EMOSLIB, from which all other differences derive.

      MIR interpolation is Finite Element-based, using 'flat' 3D triangles supporting linear shape functions (Lagrange polynomials of order 1), which is where the linear interpolation name comes from; Because it is FE-based, the supporting elements are a surface tessellation (no element overlaps). Quadrilaterals are used when suitable (such as for regular Gaussian grids, but that is a detail). The interpolation weights (w1, w2, w3) can be interpreted as area ratios w1 = A1 / A (not distances), also called the Lagrangian weights. The supporting Euclidean geometry (3D) has a fundamental difference in the distance calculation (different geometries) from which follows the calculations of these area ratios, and that's what I was intending to explain.

      EMOSLIB interpolation, instead, is supported by a spherical surface geometry, where interpolation weights can also be interpreted as area ratios wi = Ai / A (not distances), but always for quadrilaterals (where the bilinear interpolation name comes from) and where wi = Ai / A, with Ai = Δ lati2 + lati2 . The first problem with this geometry is that there is distortion towards the poles (a latitude difference represents always the same physical distance range over a sphere surface, but a longitude difference represents less physical distance closer to the poles); The second problem is, for meshes composed of non-regular quadrilaterals (general case for reduced Gaussian grids), the calculation of Δ lati and Δ loni do create area overlaps between different quadrilateral elements, therefore not constituting a surface tessellation, and therefore cannot support locally conservative interpolation.

      There is plenty of FE literature for this, and I can also add this explanation to the page. I only intended to show a fundamental difference rather than explain the calculation of the interpolation weights, but indeed it is probably very relevant to include. What do you think?