[QGIS-Developer] Exporting contour geometry from mesh layer

[Brian Haynes]

I want to be able to display certain datasets within a QgsMeshLayer as contours.  My code is set up to create a vector layer with the memory provider.  Despite the mesh layer being valid, I am not able to return valid linestrings.  Is there another way of going about this that doesn't include running a QProcess (gdal_contour)?   The part of the code in question:

QgsMeshLayer *mesh = new QgsMeshLayer("/path/to/file.nc","MyData","mdal");
QString contourField = "MSLP";

if ( mesh->isValid() )
{
     QgsMeshRendererSettings rens = mesh->rendererSettings();

     for ( int &i: mesh->datasetGroupsIndexes() )
     {
           QgsMeshDatasetIndex did = i;
           QgsMeshDatasetGroupMetadata gmd = mesh->datasetGroupMetadata(did);

            if ( gmd.name() == contourField )
            {
                 rens.setActiveScalarDatasetGroup(i);
                 mesh->setRendererSettings(rens);
                 QgsMeshContours *contours = new QgsMeshContours(mesh);
                 double cint = 2.0;
                 double max = gmd.maximum();
                 double start = gmd.minimum()-cint;
                 QgsVectorLayer *vec = new QgsVectorLayer("Linestringcrs=EPSG:4326",gmd.name(),"memory");

                 if ( vec->isValid() )
                 {
                      vec->startEditing();
                      QgsVectorDataProvider *vdp = vec->dataProvider();

                       while ( start < max )
                       {
                             start = start+cint;
                             QgsFeature feat;
                             QgsGeometry line = contours->exportLines(i,start,QgsMeshRendererScalarSettings::NeighbourAverage);

                              if ( line.isEmpty() )
                              {
                                    cout << "No valid geometry for value = " << start << endl;
                              } else
                              {
                                   //.....'create attributes/fields for feature'...
                               }
                         }
                         //...create linesymbol rendering, add layer to canvas...'
                    }
              }
       }

The output confirms that the dataset values are not empty or NaN (at least most of it), as I am returned the stdout "No valid geometry for value = 998.0, 1000.0, 1002.0,"..etc.  I've also tried to debug this by only entering a value and leaving the dataset index and resampling method out and looped it from a large negative to large positive value, nothing valid is returned.  Also to add to this exchange, if there exists any more documentation on QgsMeshLayers I'd be interested in having a link to it; seems like only rasters and vector layers have much in the way of programming tips (e.g. PyQGIS cookbook).
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-developer/attachments/20220125/a04a8be5/attachment-0001.html>