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