QGIS API Documentation  2.17.0-Master (0497e4a)
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 }
QgsPolygon asPolygon() const
Return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list...
Wrapper for iterator of features from vector data provider or vector layer.
QgsMultiPolyline asMultiPolyline() const
Return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
int createSample(QProgressDialog *pd)
QgsTransectSample(QgsVectorLayer *strataLayer, const QString &strataIdAttribute, const QString &minDistanceAttribute, const QString &nPointsAttribute, DistanceUnits minDistUnits, QgsVectorLayer *baselineLayer, bool shareBaseline, const QString &baselineStrataId, const QString &outputPointLayer, const QString &outputLineLayer, const QString &usedBaselineLayer, double minTransectLength=0.0, double baselineBufferDistance=-1.0, double baselineSimplificationTolerance=-1.0)
QGis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
QgsPoint asPoint() const
Return contents of the geometry as a point if wkbType is WKBPoint, otherwise returns [0...
bool intersects(const QgsRectangle &r) const
Test for intersection with a rectangle (uses GEOS)
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
double distance(const QgsGeometry &geom) const
Returns the minimum distance between this geometry and another geometry, using GEOS.
void append(const T &value)
void setMaximum(int maximum)
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.
QList< QgsFeatureId > intersects(const QgsRectangle &rect) const
Returns features that intersect the specified rectangle.
const_iterator constEnd() const
const QgsCoordinateReferenceSystem & crs() const
Returns layer&#39;s spatial reference system.
#define Q_NOWARN_DEPRECATED_PUSH
Definition: qgis.h:515
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
Container of fields for a vector layer.
Definition: qgsfield.h:252
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
const QgsGeometry * constGeometry() const
Gets a const pointer to the geometry object associated with this feature.
Definition: qgsfeature.cpp:82
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()
#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
QgsPolyline asPolyline() const
Return contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list...
double y() const
Get the y value of the point.
Definition: qgspoint.h:193
int mt_rand()
void reset(T *other)
QgsFields fields() const
Returns the list of fields of this layer.
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
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
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...
bool isEmpty() const
const_iterator constEnd() const
QGis::UnitType mapUnits() const
Returns the units for the projection used by the CRS.
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
compute bearing - in radians
QgsGeometry * interpolate(double distance) const
Return interpolated point on line at distance.
#define M_PI
This class wraps a request for features to a vector layer (or directly its vector data provider)...
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:346
const T & value() const
double measureLine(const QList< QgsPoint > &points) const
Measures the length of a line with multiple segments.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:44
iterator end()
QGis::GeometryType type() const
Returns type of the geometry as a QGis::GeometryType.
bool isValid()
Return the status of the layer.
T & value() const
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.
double length() const
Returns the length of geometry using GEOS.
QgsGeometry * simplify(double tolerance) const
Returns a simplified version of this geometry using a specified tolerance value.
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
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.
QTime currentTime()
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.
const T & at(int i) const
const_iterator constBegin() const
bool insertFeature(const QgsFeature &f)
Add feature to index.
WriterError hasError()
Checks whether there were any errors in constructor.
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
Class for storing a coordinate reference system (CRS)
QgsGeometry * intersection(const QgsGeometry *geometry) const
Returns a geometry representing the points shared by this geometry and other.
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 measureLength(const QgsGeometry *geometry) const
Measures the length of a geometry.
bool isValid() const
double toDouble(bool *ok) const
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:353
iterator insert(const Key &key, const T &value)
const QgsField & field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
Definition: qgsfield.cpp:427
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
long srsid() const
Returns the SrsId, if available.
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
int size() const
Represents a vector layer which manages a vector based data sets.
QVariant::Type type() const
Gets variant type of the field as it will be retrieved from data source.
Definition: qgsfield.cpp:97
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:271
QString toString() const
iterator find(const Key &key)
double x() const
Get the x value of the point.
Definition: qgspoint.h:185
void setEllipsoidalMode(bool flag)
Sets whether coordinates must be projected to ellipsoid before measuring.