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