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