Interpolation of non-rectangular data failed

Question

I want to interpolate a type of data which is on a triangular lattice in order to make DensityPlot faster(ListDensityPlot is so slow). However, the interoplation failed with error, no matter setting InterpolationOrder->1 or InterpolationOrder->All

Here is a less dense data

data = Import["https://pastebin.com/raw/6XmFDzmf"];
plotData = 
  Partition[
   StringCases[data, x : NumberString :> Internal`StringToDouble@x], 
   3];
Interpolation[plotData]

it will raise several errors

Interpolation::udeg: Interpolation on unstructured grids is currently
  only supported for InterpolationOrder->1 or InterpolationOrder->All.
  Order will be reduced to 1.
  
  Interpolation::femimq: The element mesh has insufficient quality of
  0.`. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements.
  
  Interpolation::fememtlq: The quality 0.of the underlying mesh is too
  low. The quality needs to be larger than 0..

What does it mean? How to interpolate such data or general non-rectangular data?

PS: the density plot is

yohbs · Accepted Answer

The problem is that for regular (but not rectangular) meshes the Delaunay mesh is unstable. It's a bug, also mentioned in this question, which was about rectangular grids. The same workaround works here -- just jiggle the points a tiny bit around their perfect lattice positions:

epsilon = 10^-7;
jiggledPlotData = {#1 + RandomReal[epsilon {-1, 1}], #2 + 
      RandomReal[epsilon {-1, 1}], #3} & @@@ plotData;
reg = ConvexHullMesh[jiggledPlotData[[All, 1 ;; 2]]];
f = Interpolation[jiggledPlotData, InterpolationOrder -> 1];
DensityPlot[f[x, y], {x, y} [Element] reg]

matheorem · Answer

A better method than jiggle original grid which I learned recently is dealing with mesh explicitly.
First, we need a package
Needs["NDSolve`FEM`"];

We can generate Delaunay mesh by
originalMesh = ToElementMesh[plotData[[;; , 1 ;; 2]]]

but this will throw a warning message

ToElementMesh::femimq: The element mesh has insufficient quality of
0.`. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements.

this is the same femimq problem, though the mesh looks fine at first sight as(enlarged)

The ElementMesh provides powerful tools to examing what is going on.
originalMesh["Quality"][[1]] // Histogram

The histogram shows most of the triangle has good quality close to 1. But there are plenty of triangle has quality near 0 which are bad enough.
using highlightBadMeshElement(code at the end), we can show the location of these bad triangles.
highlightBadMeshElement[originalMesh, 0.2]

They are all located at the edge shown by red marker.
We can delete these redundent bad triangles and form a refined mesh using refineMesh(code at the end)
mesh = refineMesh[originalMesh, 0.2];
f = ElementMeshInterpolation[{mesh}, plotData[[;; , -1]]];
DensityPlot[f[x, y], {x, y} [Element] mesh]

gives

functions:
highlightBadMeshElement[mesh_, qualityThreshold_] := Module[{},
  pos = Position[mesh["Quality"], _?(# <= qualityThreshold &)];
  Show[
   mesh["Wireframe"[
     "ElementMeshDirective" -> 
      Directive[EdgeForm[GrayLevel[.6]], FaceForm[]]]],
   mesh["Wireframe"[pos, "MeshElement" -> "MeshElements", 
     "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
   , Boxed -> False]]

refineMesh[mesh_, qualityThreshold_] := Module[{},
   qualityList = First@mesh["Quality"];
   pos = Flatten@Position[qualityList, _?(# > qualityThreshold &)];
   If[Length@mesh["MeshElements"] != 1,
    Message[refineMesh::MoreThanOneTypeOfElement,
     Head /@ mesh["MeshElements"]];Abort[],
    elementHead = Head@First@mesh["MeshElements"];
    ToElementMesh["Coordinates" -> mesh["Coordinates"], 
     "MeshElements" -> {elementHead[
        mesh["MeshElements"][[1, 1, pos]]]}]]];
refineMesh::MoreThanOneTypeOfElement = 
  "More than one type of element. `1`";

Interpolation of non-rectangular data failed

2 Answers

Add your own answers!

Ask a Question