QGIS API Documentation  2.99.0-Master (dcec6bb)
qgstransectsample.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstransectsample.cpp
3  ---------------------
4  begin : July 2013
5  copyright : (C) 2013 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 #include "qgstransectsample.h"
16 #include "qgsdistancearea.h"
17 #include "qgsfeatureiterator.h"
18 #include "qgsgeometry.h"
19 #include "qgsspatialindex.h"
20 #include "qgsvectorfilewriter.h"
21 #include "qgsvectorlayer.h"
22 #include "qgsproject.h"
23 #include <QProgressDialog>
24 #include <QFileInfo>
25 #ifndef _MSC_VER
26 #include <cstdint>
27 #endif
28 #include "mersenne-twister.h"
29 #include <limits>
30 
31 QgsTransectSample::QgsTransectSample( QgsVectorLayer *strataLayer, const QString &strataIdAttribute, const QString &minDistanceAttribute, const QString &nPointsAttribute, DistanceUnits minDistUnits,
32  QgsVectorLayer *baselineLayer, bool shareBaseline, const QString &baselineStrataId, const QString &outputPointLayer,
33  const QString &outputLineLayer, const QString &usedBaselineLayer, double minTransectLength,
34  double baselineBufferDistance, double baselineSimplificationTolerance )
35  : mStrataLayer( strataLayer )
36  , mStrataIdAttribute( strataIdAttribute )
37  , mMinDistanceAttribute( minDistanceAttribute )
38  , mNPointsAttribute( nPointsAttribute )
39  , mBaselineLayer( baselineLayer )
40  , mShareBaseline( shareBaseline )
41  , mBaselineStrataId( baselineStrataId )
42  , mOutputPointLayer( outputPointLayer )
43  , mOutputLineLayer( outputLineLayer )
44  , mUsedBaselineLayer( usedBaselineLayer )
45  , mMinDistanceUnits( minDistUnits )
46  , mMinTransectLength( minTransectLength )
47  , mBaselineBufferDistance( baselineBufferDistance )
48  , mBaselineSimplificationTolerance( baselineSimplificationTolerance )
49 {
50 }
51 
53  : mStrataLayer( nullptr )
54  , mBaselineLayer( nullptr )
55  , mShareBaseline( false )
56  , mMinDistanceUnits( Meters )
57  , mMinTransectLength( 0.0 )
58  , mBaselineBufferDistance( -1.0 )
59  , mBaselineSimplificationTolerance( -1.0 )
60 {
61 }
62 
63 int QgsTransectSample::createSample( QProgressDialog *pd )
64 {
65  Q_UNUSED( pd );
66 
67  if ( !mStrataLayer || !mStrataLayer->isValid() )
68  {
69  return 1;
70  }
71 
72  if ( !mBaselineLayer || !mBaselineLayer->isValid() )
73  {
74  return 2;
75  }
76 
77  //stratum id is not necessarily an integer
78  QVariant::Type stratumIdType = QVariant::Int;
79  if ( !mStrataIdAttribute.isEmpty() )
80  {
81  stratumIdType = mStrataLayer->fields().field( mStrataIdAttribute ).type();
82  }
83 
84  //create vector file writers for output
85  QgsFields outputPointFields;
86  outputPointFields.append( QgsField( QStringLiteral( "id" ), stratumIdType ) );
87  outputPointFields.append( QgsField( QStringLiteral( "station_id" ), QVariant::Int ) );
88  outputPointFields.append( QgsField( QStringLiteral( "stratum_id" ), stratumIdType ) );
89  outputPointFields.append( QgsField( QStringLiteral( "station_code" ), QVariant::String ) );
90  outputPointFields.append( QgsField( QStringLiteral( "start_lat" ), QVariant::Double ) );
91  outputPointFields.append( QgsField( QStringLiteral( "start_long" ), QVariant::Double ) );
92 
93  QgsVectorFileWriter outputPointWriter( mOutputPointLayer, QStringLiteral( "utf-8" ), outputPointFields, QgsWkbTypes::Point,
94  mStrataLayer->crs() );
95  if ( outputPointWriter.hasError() != QgsVectorFileWriter::NoError )
96  {
97  return 3;
98  }
99 
100  outputPointFields.append( QgsField( QStringLiteral( "bearing" ), QVariant::Double ) ); //add bearing attribute for lines
101  QgsVectorFileWriter outputLineWriter( mOutputLineLayer, QStringLiteral( "utf-8" ), outputPointFields, QgsWkbTypes::LineString,
102  mStrataLayer->crs() );
103  if ( outputLineWriter.hasError() != QgsVectorFileWriter::NoError )
104  {
105  return 4;
106  }
107 
108  QgsFields usedBaselineFields;
109  usedBaselineFields.append( QgsField( QStringLiteral( "stratum_id" ), stratumIdType ) );
110  usedBaselineFields.append( QgsField( QStringLiteral( "ok" ), QVariant::String ) );
111  QgsVectorFileWriter usedBaselineWriter( mUsedBaselineLayer, QStringLiteral( "utf-8" ), usedBaselineFields, QgsWkbTypes::LineString,
112  mStrataLayer->crs() );
113  if ( usedBaselineWriter.hasError() != QgsVectorFileWriter::NoError )
114  {
115  return 5;
116  }
117 
118  //debug: write clipped buffer bounds with stratum id to same directory as out_point
119  QFileInfo outputPointInfo( mOutputPointLayer );
120  QString bufferClipLineOutput = outputPointInfo.absolutePath() + "/out_buffer_clip_line.shp";
121  QgsFields bufferClipLineFields;
122  bufferClipLineFields.append( QgsField( QStringLiteral( "id" ), stratumIdType ) );
123  QgsVectorFileWriter bufferClipLineWriter( bufferClipLineOutput, QStringLiteral( "utf-8" ), bufferClipLineFields, QgsWkbTypes::LineString, mStrataLayer->crs() );
124 
125  //configure distanceArea depending on minDistance units and output CRS
126  QgsDistanceArea distanceArea;
127  distanceArea.setSourceCrs( mStrataLayer->crs() );
128  if ( mMinDistanceUnits == Meters )
129  {
130  distanceArea.setEllipsoid( QgsProject::instance()->ellipsoid() );
131  }
132  else
133  {
134  distanceArea.setEllipsoid( GEO_NONE );
135  }
136 
137  //possibility to transform output points to lat/long
139 
140  //init random number generator
141  mt_srand( QTime::currentTime().msec() );
142 
144  fr.setSubsetOfAttributes( QStringList() << mStrataIdAttribute << mMinDistanceAttribute << mNPointsAttribute, mStrataLayer->fields() );
145  QgsFeatureIterator strataIt = mStrataLayer->getFeatures( fr );
146 
147  QgsFeature fet;
148  int nTotalTransects = 0;
149  int nFeatures = 0;
150 
151  if ( pd )
152  {
153  pd->setMaximum( mStrataLayer->featureCount() );
154  }
155 
156  while ( strataIt.nextFeature( fet ) )
157  {
158  if ( pd )
159  {
160  pd->setValue( nFeatures );
161  }
162  if ( pd && pd->wasCanceled() )
163  {
164  break;
165  }
166 
167  if ( !fet.hasGeometry() )
168  {
169  continue;
170  }
171  QgsGeometry strataGeom = fet.geometry();
172 
173  //find baseline for strata
174  QVariant strataId = fet.attribute( mStrataIdAttribute );
175  QgsGeometry baselineGeom = findBaselineGeometry( strataId.isValid() ? strataId : -1 );
176  if ( baselineGeom.isNull() )
177  {
178  continue;
179  }
180 
181  double minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
182  double minDistanceLayerUnits = minDistance;
183  //if minDistance is in meters and the data in degrees, we need to apply a rough conversion for the buffer distance
184  double bufferDist = bufferDistance( minDistance );
185  if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QgsUnitTypes::DistanceDegrees )
186  {
187  minDistanceLayerUnits = minDistance / 111319.9;
188  }
189 
190  QgsGeometry clippedBaseline = strataGeom.intersection( baselineGeom );
191  if ( !clippedBaseline || clippedBaseline.wkbType() == QgsWkbTypes::Unknown )
192  {
193  continue;
194  }
195  QgsGeometry *bufferLineClipped = clipBufferLine( strataGeom, &clippedBaseline, bufferDist );
196  if ( !bufferLineClipped )
197  {
198  continue;
199  }
200 
201  //save clipped baseline to file
202  QgsFeature blFeature( usedBaselineFields );
203  blFeature.setGeometry( clippedBaseline );
204  blFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
205  blFeature.setAttribute( QStringLiteral( "ok" ), "f" );
206  usedBaselineWriter.addFeature( blFeature );
207 
208  //start loop to create random points along the baseline
209  int nTransects = fet.attribute( mNPointsAttribute ).toInt();
210  int nCreatedTransects = 0;
211  int nIterations = 0;
212  int nMaxIterations = nTransects * 50;
213 
214  QgsSpatialIndex sIndex; //to check minimum distance
215  QMap< QgsFeatureId, QgsGeometry > lineFeatureMap;
216 
217  while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
218  {
219  double randomPosition = ( ( double )mt_rand() / MD_RAND_MAX ) * clippedBaseline.length();
220  QgsGeometry samplePoint = clippedBaseline.interpolate( randomPosition );
221  ++nIterations;
222  if ( samplePoint.isNull() )
223  {
224  continue;
225  }
226  QgsPointXY sampleQgsPointXY = samplePoint.asPoint();
227  QgsPointXY latLongSamplePoint = toLatLongTransform.transform( sampleQgsPointXY );
228 
229  QgsFeature samplePointFeature( outputPointFields );
230  samplePointFeature.setGeometry( samplePoint );
231  samplePointFeature.setAttribute( QStringLiteral( "id" ), nTotalTransects + 1 );
232  samplePointFeature.setAttribute( QStringLiteral( "station_id" ), nCreatedTransects + 1 );
233  samplePointFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
234  samplePointFeature.setAttribute( QStringLiteral( "station_code" ), strataId.toString() + '_' + QString::number( nCreatedTransects + 1 ) );
235  samplePointFeature.setAttribute( QStringLiteral( "start_lat" ), latLongSamplePoint.y() );
236  samplePointFeature.setAttribute( QStringLiteral( "start_long" ), latLongSamplePoint.x() );
237 
238  //find closest point on clipped buffer line
239  QgsPointXY minDistPoint;
240 
241  int afterVertex;
242  if ( bufferLineClipped->closestSegmentWithContext( sampleQgsPointXY, minDistPoint, afterVertex ) < 0 )
243  {
244  continue;
245  }
246 
247  //bearing between sample point and min dist point (transect direction)
248  double bearing = distanceArea.bearing( sampleQgsPointXY, minDistPoint ) / M_PI * 180.0;
249 
250  QgsPointXY ptFarAway( sampleQgsPointXY.x() + ( minDistPoint.x() - sampleQgsPointXY.x() ) * 1000000,
251  sampleQgsPointXY.y() + ( minDistPoint.y() - sampleQgsPointXY.y() ) * 1000000 );
252  QgsPolyline lineFarAway;
253  lineFarAway << sampleQgsPointXY << ptFarAway;
254  QgsGeometry lineFarAwayGeom = QgsGeometry::fromPolyline( lineFarAway );
255  QgsGeometry lineClipStratum = lineFarAwayGeom.intersection( strataGeom );
256  if ( lineClipStratum.isNull() )
257  {
258  continue;
259  }
260 
261  //cancel if distance between sample point and line is too large (line does not start at point)
262  if ( lineClipStratum.distance( samplePoint ) > 0.000001 )
263  {
264  continue;
265  }
266 
267  //if lineClipStratum is a multiline, take the part line closest to sampleQgsPoint
268  if ( lineClipStratum.wkbType() == QgsWkbTypes::MultiLineString
269  || lineClipStratum.wkbType() == QgsWkbTypes::MultiLineString25D )
270  {
271  QgsGeometry singleLine = closestMultilineElement( sampleQgsPointXY, lineClipStratum );
272  if ( !singleLine.isNull() )
273  {
274  lineClipStratum = singleLine;
275  }
276  }
277 
278  //cancel if length of lineClipStratum is too small
279  double transectLength = distanceArea.measureLength( lineClipStratum );
280  if ( transectLength < mMinTransectLength )
281  {
282  continue;
283  }
284 
285  //search closest existing profile. Cancel if dist < minDist
286  if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
287  {
288  continue;
289  }
290 
291  QgsFeatureId fid( nCreatedTransects );
292  QgsFeature sampleLineFeature( outputPointFields, fid );
293  sampleLineFeature.setGeometry( lineClipStratum );
294  sampleLineFeature.setAttribute( QStringLiteral( "id" ), nTotalTransects + 1 );
295  sampleLineFeature.setAttribute( QStringLiteral( "station_id" ), nCreatedTransects + 1 );
296  sampleLineFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
297  sampleLineFeature.setAttribute( QStringLiteral( "station_code" ), strataId.toString() + '_' + QString::number( nCreatedTransects + 1 ) );
298  sampleLineFeature.setAttribute( QStringLiteral( "start_lat" ), latLongSamplePoint.y() );
299  sampleLineFeature.setAttribute( QStringLiteral( "start_long" ), latLongSamplePoint.x() );
300  sampleLineFeature.setAttribute( QStringLiteral( "bearing" ), bearing );
301  outputLineWriter.addFeature( sampleLineFeature );
302 
303  //add point to file writer here.
304  //It can only be written if the corresponding transect has been as well
305  outputPointWriter.addFeature( samplePointFeature );
306 
307  sIndex.insertFeature( sampleLineFeature );
308  lineFeatureMap.insert( fid, lineClipStratum );
309 
310  ++nTotalTransects;
311  ++nCreatedTransects;
312  }
313 
314  QgsFeature bufferClipFeature;
315  bufferClipFeature.setGeometry( *bufferLineClipped );
316  delete bufferLineClipped;
317  bufferClipFeature.setAttribute( QStringLiteral( "id" ), strataId );
318  bufferClipLineWriter.addFeature( bufferClipFeature );
319  //delete bufferLineClipped;
320 
321  ++nFeatures;
322  }
323 
324  if ( pd )
325  {
326  pd->setValue( mStrataLayer->featureCount() );
327  }
328 
329  return 0;
330 }
331 
332 QgsGeometry QgsTransectSample::findBaselineGeometry( const QVariant &strataId )
333 {
334  if ( !mBaselineLayer )
335  {
336  return QgsGeometry();
337  }
338 
339  QgsFeatureIterator baseLineIt = mBaselineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QStringList( mBaselineStrataId ), mBaselineLayer->fields() ) );
340  QgsFeature fet;
341 
342  while ( baseLineIt.nextFeature( fet ) ) //todo: cache this in case there are many baslines
343  {
344  if ( strataId == fet.attribute( mBaselineStrataId ) || mShareBaseline )
345  {
346  return fet.geometry();
347  }
348  }
349  return QgsGeometry();
350 }
351 
352 bool QgsTransectSample::otherTransectWithinDistance( const QgsGeometry &geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex &sIndex,
353  const QMap< QgsFeatureId, QgsGeometry > &lineFeatureMap, QgsDistanceArea &da )
354 {
355  if ( geom.isNull() )
356  {
357  return false;
358  }
359 
360  QgsGeometry buffer = geom.buffer( minDistLayerUnit, 8 );
361  if ( buffer.isNull() )
362  {
363  return false;
364  }
365  QgsRectangle rect = buffer.boundingBox();
366  QList<QgsFeatureId> lineIdList = sIndex.intersects( rect );
367 
368  QList<QgsFeatureId>::const_iterator lineIdIt = lineIdList.constBegin();
369  for ( ; lineIdIt != lineIdList.constEnd(); ++lineIdIt )
370  {
371  const QMap< QgsFeatureId, QgsGeometry >::const_iterator idMapIt = lineFeatureMap.find( *lineIdIt );
372  if ( idMapIt != lineFeatureMap.constEnd() )
373  {
374  double dist = 0;
375  QgsPointXY pt1, pt2;
376  closestSegmentPoints( geom, idMapIt.value(), dist, pt1, pt2 );
377  dist = da.measureLine( pt1, pt2 ); //convert degrees to meters if necessary
378 
379  if ( dist < minDistance )
380  {
381  return true;
382  }
383  }
384  }
385 
386  return false;
387 }
388 
389 bool QgsTransectSample::closestSegmentPoints( const QgsGeometry &g1, const QgsGeometry &g2, double &dist, QgsPointXY &pt1, QgsPointXY &pt2 )
390 {
391  QgsWkbTypes::Type t1 = g1.wkbType();
393  {
394  return false;
395  }
396 
397  QgsWkbTypes::Type t2 = g2.wkbType();
399  {
400  return false;
401  }
402 
403  QgsPolyline pl1 = g1.asPolyline();
404  QgsPolyline pl2 = g2.asPolyline();
405 
406  if ( pl1.size() < 2 || pl2.size() < 2 )
407  {
408  return false;
409  }
410 
411  QgsPointXY p11 = pl1.at( 0 );
412  QgsPointXY p12 = pl1.at( 1 );
413  QgsPointXY p21 = pl2.at( 0 );
414  QgsPointXY p22 = pl2.at( 1 );
415 
416  double p1x = p11.x();
417  double p1y = p11.y();
418  double v1x = p12.x() - p11.x();
419  double v1y = p12.y() - p11.y();
420  double p2x = p21.x();
421  double p2y = p21.y();
422  double v2x = p22.x() - p21.x();
423  double v2y = p22.y() - p21.y();
424 
425  double denominatorU = v2x * v1y - v2y * v1x;
426  double denominatorT = v1x * v2y - v1y * v2x;
427 
428  if ( qgsDoubleNear( denominatorU, 0 ) || qgsDoubleNear( denominatorT, 0 ) )
429  {
430  //lines are parallel
431  //project all points on the other segment and take the one with the smallest distance
432  QgsPointXY minDistPoint1;
433  double d1 = p11.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint1 );
434  QgsPointXY minDistPoint2;
435  double d2 = p12.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint2 );
436  QgsPointXY minDistPoint3;
437  double d3 = p21.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint3 );
438  QgsPointXY minDistPoint4;
439  double d4 = p22.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint4 );
440 
441  if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
442  {
443  dist = sqrt( d1 );
444  pt1 = p11;
445  pt2 = minDistPoint1;
446  return true;
447  }
448  else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
449  {
450  dist = sqrt( d2 );
451  pt1 = p12;
452  pt2 = minDistPoint2;
453  return true;
454  }
455  else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
456  {
457  dist = sqrt( d3 );
458  pt1 = p21;
459  pt2 = minDistPoint3;
460  return true;
461  }
462  else
463  {
464  dist = sqrt( d4 );
465  pt1 = p21;
466  pt2 = minDistPoint4;
467  return true;
468  }
469  }
470 
471  double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
472  double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
473 
474  if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
475  {
476  dist = 0;
477  pt1.setX( p2x + u * v2x );
478  pt1.setY( p2y + u * v2y );
479  pt2 = pt1;
480  dist = 0;
481  return true;
482  }
483 
484  if ( t > 1.0 )
485  {
486  pt1.setX( p12.x() );
487  pt1.setY( p12.y() );
488  }
489  else if ( t < 0.0 )
490  {
491  pt1.setX( p11.x() );
492  pt1.setY( p11.y() );
493  }
494  if ( u > 1.0 )
495  {
496  pt2.setX( p22.x() );
497  pt2.setY( p22.y() );
498  }
499  if ( u < 0.0 )
500  {
501  pt2.setX( p21.x() );
502  pt2.setY( p21.y() );
503  }
504  if ( t >= 0.0 && t <= 1.0 )
505  {
506  //project pt2 onto g1
507  pt2.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), pt1 );
508  }
509  if ( u >= 0.0 && u <= 1.0 )
510  {
511  //project pt1 onto g2
512  pt1.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), pt2 );
513  }
514 
515  dist = sqrt( pt1.sqrDist( pt2 ) );
516  return true;
517 }
518 
519 QgsGeometry QgsTransectSample::closestMultilineElement( const QgsPointXY &pt, const QgsGeometry &multiLine )
520 {
521  if ( !multiLine || ( multiLine.wkbType() != QgsWkbTypes::MultiLineString
522  && multiLine.wkbType() != QgsWkbTypes::MultiLineString25D ) )
523  {
524  return QgsGeometry();
525  }
526 
527  double minDist = DBL_MAX;
528  double currentDist = 0;
529  QgsGeometry currentLine;
530  QgsGeometry closestLine;
531  QgsGeometry pointGeom = QgsGeometry::fromPoint( pt );
532 
533  QgsMultiPolyline multiPolyline = multiLine.asMultiPolyline();
534  QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
535  for ( ; it != multiPolyline.constEnd(); ++it )
536  {
537  currentLine = QgsGeometry::fromPolyline( *it );
538  currentDist = pointGeom.distance( currentLine );
539  if ( currentDist < minDist )
540  {
541  minDist = currentDist;
542  closestLine = currentLine;
543  }
544  }
545 
546  return closestLine;
547 }
548 
549 QgsGeometry *QgsTransectSample::clipBufferLine( const QgsGeometry &stratumGeom, QgsGeometry *clippedBaseline, double tolerance )
550 {
551  if ( !stratumGeom || !clippedBaseline || clippedBaseline->wkbType() == QgsWkbTypes::Unknown )
552  {
553  return nullptr;
554  }
555 
556  QgsGeometry usedBaseline = *clippedBaseline;
557  if ( mBaselineSimplificationTolerance >= 0 )
558  {
559  //int verticesBefore = usedBaseline->asMultiPolyline().count();
560  usedBaseline = clippedBaseline->simplify( mBaselineSimplificationTolerance );
561  if ( usedBaseline.isNull() )
562  {
563  return nullptr;
564  }
565  //int verticesAfter = usedBaseline->asMultiPolyline().count();
566 
567  //debug: write to file
568  /*QgsVectorFileWriter debugWriter( "/tmp/debug.shp", "utf-8", QgsFields(), QgsWkbTypes::LineString, &( mStrataLayer->crs() ) );
569  QgsFeature debugFeature; debugFeature.setGeometry( usedBaseline );
570  debugWriter.addFeature( debugFeature );*/
571  }
572 
573  double currentBufferDist = tolerance;
574  int maxLoops = 10;
575 
576  for ( int i = 0; i < maxLoops; ++i )
577  {
578  //loop with tolerance: create buffer, convert buffer to line, clip line by stratum, test if result is (single) line
579  QgsGeometry clipBaselineBuffer = usedBaseline.buffer( currentBufferDist, 8 );
580  if ( clipBaselineBuffer.isNull() )
581  {
582  continue;
583  }
584 
585  //it is also possible that clipBaselineBuffer is a multipolygon
586  QgsGeometry bufferLine; //buffer line or multiline
587  QgsGeometry bufferLineClipped;
588  QgsMultiPolyline mpl;
589  if ( clipBaselineBuffer.isMultipart() )
590  {
591  QgsMultiPolygon bufferMultiPolygon = clipBaselineBuffer.asMultiPolygon();
592  if ( bufferMultiPolygon.size() < 1 )
593  {
594  continue;
595  }
596 
597  for ( int j = 0; j < bufferMultiPolygon.size(); ++j )
598  {
599  int size = bufferMultiPolygon.at( j ).size();
600  for ( int k = 0; k < size; ++k )
601  {
602  mpl.append( bufferMultiPolygon.at( j ).at( k ) );
603  }
604  }
605  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
606  }
607  else
608  {
609  QgsPolygon bufferPolygon = clipBaselineBuffer.asPolygon();
610  if ( bufferPolygon.size() < 1 )
611  {
612  continue;
613  }
614 
615  int size = bufferPolygon.size();
616  mpl.reserve( size );
617  for ( int j = 0; j < size; ++j )
618  {
619  mpl.append( bufferPolygon[j] );
620  }
621  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
622  }
623  bufferLineClipped = bufferLine.intersection( stratumGeom );
624 
625  if ( bufferLineClipped.isNull() && bufferLineClipped.type() == QgsWkbTypes::LineGeometry )
626  {
627  //if stratumGeom is a multipolygon, bufferLineClipped must intersect each part
628  bool bufferLineClippedIntersectsStratum = true;
629  if ( stratumGeom.wkbType() == QgsWkbTypes::MultiPolygon || stratumGeom.wkbType() == QgsWkbTypes::MultiPolygon25D )
630  {
631  QVector<QgsPolygon> multiPoly = stratumGeom.asMultiPolygon();
632  QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
633  for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
634  {
635  QgsGeometry poly = QgsGeometry::fromPolygon( *multiIt );
636  if ( !poly.intersects( bufferLineClipped ) )
637  {
638  bufferLineClippedIntersectsStratum = false;
639  break;
640  }
641  }
642  }
643 
644  if ( bufferLineClippedIntersectsStratum )
645  {
646  return new QgsGeometry( bufferLineClipped );
647  }
648  }
649 
650  currentBufferDist /= 2;
651  }
652 
653  return nullptr; //no solution found even with reduced tolerances
654 }
655 
656 double QgsTransectSample::bufferDistance( double minDistanceFromAttribute ) const
657 {
658  double bufferDist = minDistanceFromAttribute;
659  if ( mBaselineBufferDistance >= 0 )
660  {
661  bufferDist = mBaselineBufferDistance;
662  }
663 
664  if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QgsUnitTypes::DistanceDegrees )
665  {
666  bufferDist /= 111319.9;
667  }
668 
669  return bufferDist;
670 }
QgsPolygon asPolygon() const
Return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list...
Wrapper for iterator of features from vector data provider or vector layer.
QgsMultiPolyline asMultiPolyline() const
Return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list.
A rectangle specified with double values.
Definition: qgsrectangle.h:38
int createSample(QProgressDialog *pd)
QgsTransectSample(QgsVectorLayer *strataLayer, const QString &strataIdAttribute, const QString &minDistanceAttribute, const QString &nPointsAttribute, DistanceUnits minDistUnits, QgsVectorLayer *baselineLayer, bool shareBaseline, const QString &baselineStrataId, const QString &outputPointLayer, const QString &outputLineLayer, const QString &usedBaselineLayer, double minTransectLength=0.0, double baselineBufferDistance=-1.0, double baselineSimplificationTolerance=-1.0)
static QgsGeometry fromPolyline(const QgsPolyline &polyline)
Creates a new geometry from a QgsPolyline object.
bool intersects(const QgsRectangle &r) const
Test for intersection with a rectangle (uses GEOS)
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
double distance(const QgsGeometry &geom) const
Returns the minimum distance between this geometry and another geometry, using GEOS.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double y
Definition: qgspointxy.h:47
A class to represent a 2D point.
Definition: qgspointxy.h:42
QList< QgsFeatureId > intersects(const QgsRectangle &rect) const
Returns features that intersect the specified rectangle.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QgsGeometry interpolate(double distance) const
Return interpolated point on line at distance.
Container of fields for a vector layer.
Definition: qgsfields.h:41
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:96
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:204
QgsGeometry buffer(double distance, int segments) const
Returns a buffer region around this geometry having the given width and with a specified number of se...
A convenience class for writing vector files to disk.
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:61
bool isValid() const
Return the status of the layer.
double closestSegmentWithContext(const QgsPointXY &point, QgsPointXY &minDistPoint, int &afterVertex, double *leftOf=nullptr, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Searches for the closest segment of geometry to the given point.
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:190
#define MD_RAND_MAX
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:203
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:71
QgsPolyline asPolyline() const
Return contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list...
QVector< QgsPointXY > QgsPolyline
Polyline is represented as a vector of points.
Definition: qgsgeometry.h:48
int mt_rand()
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.cpp:264
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:66
double sqrDistToSegment(double x1, double y1, double x2, double y2, QgsPointXY &minDistPoint, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Returns the minimum distance between this point and a segment.
Definition: qgspointxy.cpp:332
QgsFields fields() const override
Returns the list of fields of this layer.
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
void setY(double y)
Sets the y value of the point.
Definition: qgspointxy.h:108
#define M_PI
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QVector< QgsPolygon > QgsMultiPolygon
A collection of QgsPolygons that share a common collection of attributes.
Definition: qgsgeometry.h:76
Degrees, for planar geographic CRS distance measurements.
Definition: qgsunittypes.h:50
void mt_srand(unsigned value)
QgsCoordinateReferenceSystem crs() const
Returns the layer&#39;s spatial reference system.
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Append a field. The field must have unique name, otherwise it is rejected (returns false) ...
Definition: qgsfields.cpp:59
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:46
QVector< QgsPolyline > QgsPolygon
Polygon: first item of the list is outer ring, inner rings (if any) start from second item...
Definition: qgsgeometry.h:55
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
QgsWkbTypes::GeometryType type() const
Returns type of the geometry as a QgsWkbTypes::GeometryType.
double length() const
Returns the length of geometry using GEOS.
void setX(double x)
Sets the x value of the point.
Definition: qgspointxy.h:100
double x
Definition: qgspointxy.h:46
static QgsGeometry fromPolygon(const QgsPolygon &polygon)
Creates a new geometry from a QgsPolygon.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const override
Query the layer for features specified in request.
QVector< QgsPolyline > QgsMultiPolyline
A collection of QgsPolylines that share a common collection of attributes.
Definition: qgsgeometry.h:69
A general purpose distance and area calculator, capable of performing ellipsoid based calculations...
QgsGeometry simplify(double tolerance) const
Returns a simplified version of this geometry using a specified tolerance value.
bool insertFeature(const QgsFeature &f)
Add feature to index.
QgsPointXY asPoint() const
Return contents of the geometry as a point if wkbType is WKBPoint, otherwise returns [0...
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:377
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
This class represents a coordinate reference system (CRS).
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:137
Class for doing transforms between two map coordinate systems.
qint64 QgsFeatureId
Definition: qgsfeature.h:37
static QgsGeometry fromPoint(const QgsPointXY &point)
Creates a new geometry from a QgsPointXY object.
void setSourceCrs(const QgsCoordinateReferenceSystem &srcCRS)
Sets source spatial reference system.
bool nextFeature(QgsFeature &f)
Represents a vector layer which manages a vector based data sets.
QVariant::Type type() const
Gets variant type of the field as it will be retrieved from data source.
Definition: qgsfield.cpp:94
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:255
double measureLine(const QList< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
QgsField field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
Definition: qgsfields.cpp:140
static QgsGeometry fromMultiPolyline(const QgsMultiPolyline &multiline)
Creates a new geometry from a QgsMultiPolyline object.