Hfrhyu

I created a plugin few years ago in Python. The plugin uses QGIS interpolation and works for QGIS 2.18. But because the QGIS interpolation creates raster (grid) and I need only Triangulation (TIN as lines) I checked the QGIS code and call function int QgsTinInterpolator::interpolatePoint( double x, double y, double &result, QgsFeedback * ) via PyQGIS.

The Python code for 2.18 QGIS API:

ld_layer = QgsInterpolator.LayerData()
ld_layer.vectorLayer = layer
ld_layer.zCoordInterpolation = False
ld_layer.interpolationAttribute = index
ld_layer.mInputType = 1

ld2_layer = QgsInterpolator.LayerData()
ld2_layer.vectorLayer = layer2
ld2_layer.zCoordInterpolation = False
ld2_layer.interpolationAttribute = index
ld2_layer.mInputType = 1

itp = QgsTINInterpolator(lds)
itp.setExportTriangulationToFile(True)

itp.setTriangulationFilePath(self.path)

rect = lds[0].vectorLayer.extent()

current_y_value = rect.yMaximum() - self.resolution / 2.0
current_x_value = rect.xMinimum() + self.resolution / 2.0
itp.interpolatePoint(current_x_value, current_y_value)

I created LayerData objects and then I called setExportTriangulationToFile(True) and setTriangulationFilePath(self.path).
After calling itp.interpolatePoint(current_x_value, current_y_value) the Triangulation filed was created. It is because of the void QgsTinInterpolator::initialize() function. This function is inside the qgstininterpolator.cpp and because it is private function I can't call that function via PyQGIS API, so I decided call interpolatePoint() which calls initialize() if initialization is not configured.

But now mainly because of in the plugin there is a performance need and because I am really interesting in exploring the world of C/C++ I decided few months ago learn C++ and rewrite the plugin. There are some changes in API:

• there is no setExportTriangulationToFile(True)
  


• there is no setTriangulationFilePath(self.path)
  


• instead of there is setTriangulationSink(QgsFeatureSink)
  

So I started hacking with QGIS and I wrote some code in C++:

QgsVectorLayer *exp5Layer = new QgsVectorLayer("/home/username/Downloads/data/testing_data/exp_q5_1.shp", "exp_q5", "ogr");
QgsVectorLayer *exp20Layer = new QgsVectorLayer("/home/username/Downloads/data/testing_data/exp_q20_1.shp", "exp_q20", "ogr");

const int index = 0;
QgsInterpolator::LayerData ld5;
QgsInterpolator::LayerData ld20;

ld5.source = exp5Layer;
ld5.interpolationAttribute = index;
ld5.sourceType = QgsInterpolator::SourceType::SourceStructureLines;
ld5.valueSource = QgsInterpolator::ValueSource::ValueZ;

ld20.source = exp20Layer;
ld20.interpolationAttribute = index;
ld20.sourceType = QgsInterpolator::SourceType::SourceStructureLines;
ld20.valueSource = QgsInterpolator::ValueSource::ValueZ;

QList<QgsInterpolator::LayerData> lds;
lds.append(ld5);
lds.append(ld20);

QgsTinInterpolator *itp = new QgsTinInterpolator(lds);

QgsWkbTypes::Type wkbType = exp5Layer->wkbType();
QgsFields fields = exp5Layer->fields();
QgsVectorLayer *sink = new QgsVectorLayer("/home/username/Downloads/data/testing_data/triangulation.shp", "triangulation", "ogr");
itp->setTriangulationSink(sink);

double current_y_value = rect.yMaximum() - 0.1 / 2.0;
double current_x_value = rect.xMinimum() + 0.1 / 2.0;
double interpolatedValue;
QgsFeedback *feedback = new QgsFeedback();

itp->interpolatePoint(current_x_value, current_y_value, interpolatedValue, feedback);
sink->~QgsVectorFileWriter();

The code is pretty similiar to Python code but of course there are changes in the API of QGIS 3.7. But the problem is the triangulation file (triangulation.shp) is created but there are not features.

It is because in bool DualEdgeTriangulation::saveTriangulation( QgsFeatureSink *sink, QgsFeedback *feedback ) const function in DualEdgeInterpolation.cpp there is a cycle (at line 2997) and because of mHalfEdge.size() == 0 addFeature function sink->addFeature( edgeLineFeature, QgsFeatureSink::FastInsert ); (at line 3031) is not called).

There is no error and triangulation file is created but without features. If I run interpolation from the QGIS GUI with same parameters the triangulation file is created with features correctly.

Is there any problem with new API (like more changes that I didn't notice) or there is any other problem? I can't figure why mHalfEdge is empty.

I continued with the investigating the problem and I found there was some big changes in QgsTINInterpolator::initialize() function (commit). Maybe this can lead to solution or some clue.

I continue with going through the code and it is not a bug but change in the API (and mainly my fault) which leads to the problems with triangulation.

In python there is zCoordInterpolation = False property but in the new API there is an enumeration so instead of ValueZ (which I wrongly set up) you need set up the different value from enumeration (of course based on your data). So in the C++ it is ValueAttribute (in my specific case):

layerDataVariable.valueSource = QgsInterpolator::ValueSource::ValueAttribute;