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