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

Answer 1

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]

Answer 2

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`";