QGIS API Documentation  2.11.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 
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 );
292  lineFeatureMap.insert( fid, sampleLineFeature.geometryAndOwnership() );
294 
295  delete lineFarAwayGeom;
296  ++nTotalTransects;
297  ++nCreatedTransects;
298  }
299  delete clippedBaseline;
300 
301  QgsFeature bufferClipFeature;
302  bufferClipFeature.setGeometry( bufferLineClipped );
303  bufferClipFeature.setAttribute( "id", strataId );
304  bufferClipLineWriter.addFeature( bufferClipFeature );
305  //delete bufferLineClipped;
306 
307  //delete all line geometries in spatial index
308  QMap< QgsFeatureId, QgsGeometry* >::iterator featureMapIt = lineFeatureMap.begin();
309  for ( ; featureMapIt != lineFeatureMap.end(); ++featureMapIt )
310  {
311  delete( featureMapIt.value() );
312  }
313  lineFeatureMap.clear();
314  delete baselineGeom;
315 
316  ++nFeatures;
317  }
318 
319  if ( pd )
320  {
321  pd->setValue( mStrataLayer->featureCount() );
322  }
323 
324  return 0;
325 }
326 
327 QgsGeometry* QgsTransectSample::findBaselineGeometry( QVariant strataId )
328 {
329  if ( !mBaselineLayer )
330  {
331  return 0;
332  }
333 
334  QgsFeatureIterator baseLineIt = mBaselineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QStringList( mBaselineStrataId ), mBaselineLayer->pendingFields() ) );
335  QgsFeature fet;
336 
337  while ( baseLineIt.nextFeature( fet ) ) //todo: cache this in case there are many baslines
338  {
339  if ( strataId == fet.attribute( mBaselineStrataId ) || mShareBaseline )
340  {
342  return fet.geometryAndOwnership();
344  }
345  }
346  return 0;
347 }
348 
349 bool QgsTransectSample::otherTransectWithinDistance( QgsGeometry* geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex& sIndex,
350  const QMap< QgsFeatureId, QgsGeometry* >& lineFeatureMap, QgsDistanceArea& da )
351 {
352  if ( !geom )
353  {
354  return false;
355  }
356 
357  QgsGeometry* buffer = geom->buffer( minDistLayerUnit, 8 );
358  if ( !buffer )
359  {
360  return false;
361  }
362  QgsRectangle rect = buffer->boundingBox();
363  QList<QgsFeatureId> lineIdList = sIndex.intersects( rect );
364 
365  QList<QgsFeatureId>::const_iterator lineIdIt = lineIdList.constBegin();
366  for ( ; lineIdIt != lineIdList.constEnd(); ++lineIdIt )
367  {
368  const QMap< QgsFeatureId, QgsGeometry* >::const_iterator idMapIt = lineFeatureMap.find( *lineIdIt );
369  if ( idMapIt != lineFeatureMap.constEnd() )
370  {
371  double dist = 0;
372  QgsPoint pt1, pt2;
373  closestSegmentPoints( *geom, *( idMapIt.value() ), dist, pt1, pt2 );
374  dist = da.measureLine( pt1, pt2 ); //convert degrees to meters if necessary
375 
376  if ( dist < minDistance )
377  {
378  delete buffer;
379  return true;
380  }
381  }
382  }
383 
384  delete buffer;
385  return false;
386 }
387 
388 bool QgsTransectSample::closestSegmentPoints( QgsGeometry& g1, QgsGeometry& g2, double& dist, QgsPoint& pt1, QgsPoint& pt2 )
389 {
390  QGis::WkbType t1 = g1.wkbType();
391  if ( t1 != QGis::WKBLineString && t1 != QGis::WKBLineString25D )
392  {
393  return false;
394  }
395 
396  QGis::WkbType t2 = g2.wkbType();
397  if ( t2 != QGis::WKBLineString && t2 != QGis::WKBLineString25D )
398  {
399  return false;
400  }
401 
402  QgsPolyline pl1 = g1.asPolyline();
403  QgsPolyline pl2 = g2.asPolyline();
404 
405  if ( pl1.size() < 2 || pl2.size() < 2 )
406  {
407  return false;
408  }
409 
410  QgsPoint p11 = pl1.at( 0 );
411  QgsPoint p12 = pl1.at( 1 );
412  QgsPoint p21 = pl2.at( 0 );
413  QgsPoint p22 = pl2.at( 1 );
414 
415  double p1x = p11.x();
416  double p1y = p11.y();
417  double v1x = p12.x() - p11.x();
418  double v1y = p12.y() - p11.y();
419  double p2x = p21.x();
420  double p2y = p21.y();
421  double v2x = p22.x() - p21.x();
422  double v2y = p22.y() - p21.y();
423 
424  double denominatorU = v2x * v1y - v2y * v1x;
425  double denominatorT = v1x * v2y - v1y * v2x;
426 
427  if ( qgsDoubleNear( denominatorU, 0 ) || qgsDoubleNear( denominatorT, 0 ) )
428  {
429  //lines are parallel
430  //project all points on the other segment and take the one with the smallest distance
431  QgsPoint minDistPoint1;
432  double d1 = p11.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint1 );
433  QgsPoint minDistPoint2;
434  double d2 = p12.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint2 );
435  QgsPoint minDistPoint3;
436  double d3 = p21.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint3 );
437  QgsPoint minDistPoint4;
438  double d4 = p22.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint4 );
439 
440  if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
441  {
442  dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
443  return true;
444  }
445  else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
446  {
447  dist = sqrt( d2 ); pt1 = p12; pt2 = minDistPoint2;
448  return true;
449  }
450  else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
451  {
452  dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
453  return true;
454  }
455  else
456  {
457  dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
458  return true;
459  }
460  }
461 
462  double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
463  double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
464 
465  if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
466  {
467  dist = 0;
468  pt1.setX( p2x + u * v2x );
469  pt1.setY( p2y + u * v2y );
470  pt2 = pt1;
471  dist = 0;
472  return true;
473  }
474 
475  if ( t > 1.0 )
476  {
477  pt1.setX( p12.x() );
478  pt1.setY( p12.y() );
479  }
480  else if ( t < 0.0 )
481  {
482  pt1.setX( p11.x() );
483  pt1.setY( p11.y() );
484  }
485  if ( u > 1.0 )
486  {
487  pt2.setX( p22.x() );
488  pt2.setY( p22.y() );
489  }
490  if ( u < 0.0 )
491  {
492  pt2.setX( p21.x() );
493  pt2.setY( p21.y() );
494  }
495  if ( t >= 0.0 && t <= 1.0 )
496  {
497  //project pt2 onto g1
498  pt2.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), pt1 );
499  }
500  if ( u >= 0.0 && u <= 1.0 )
501  {
502  //project pt1 onto g2
503  pt1.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), pt2 );
504  }
505 
506  dist = sqrt( pt1.sqrDist( pt2 ) );
507  return true;
508 }
509 
510 QgsGeometry* QgsTransectSample::closestMultilineElement( const QgsPoint& pt, QgsGeometry* multiLine )
511 {
512  if ( !multiLine || ( multiLine->wkbType() != QGis::WKBMultiLineString
513  && multiLine->wkbType() != QGis::WKBMultiLineString25D ) )
514  {
515  return 0;
516  }
517 
518  double minDist = DBL_MAX;
519  double currentDist = 0;
520  QgsGeometry* currentLine = 0;
521  QScopedPointer<QgsGeometry> closestLine;
522  QgsGeometry* pointGeom = QgsGeometry::fromPoint( pt );
523 
524  QgsMultiPolyline multiPolyline = multiLine->asMultiPolyline();
525  QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
526  for ( ; it != multiPolyline.constEnd(); ++it )
527  {
528  currentLine = QgsGeometry::fromPolyline( *it );
529  currentDist = pointGeom->distance( *currentLine );
530  if ( currentDist < minDist )
531  {
532  minDist = currentDist;
533  closestLine.reset( currentLine );
534  }
535  else
536  {
537  delete currentLine;
538  }
539  }
540 
541  delete pointGeom;
542  return closestLine.take();
543 }
544 
545 QgsGeometry* QgsTransectSample::clipBufferLine( const QgsGeometry* stratumGeom, QgsGeometry* clippedBaseline, double tolerance )
546 {
547  if ( !stratumGeom || !clippedBaseline || clippedBaseline->wkbType() == QGis::WKBUnknown )
548  {
549  return 0;
550  }
551 
552  double currentBufferDist = tolerance;
553  int maxLoops = 10;
554 
555  for ( int i = 0; i < maxLoops; ++i )
556  {
557  //loop with tolerance: create buffer, convert buffer to line, clip line by stratum, test if result is (single) line
558  QgsGeometry* clipBaselineBuffer = clippedBaseline->buffer( currentBufferDist, 8 );
559  if ( !clipBaselineBuffer )
560  {
561  delete clipBaselineBuffer;
562  continue;
563  }
564 
565  //it is also possible that clipBaselineBuffer is a multipolygon
566  QgsGeometry* bufferLine = 0; //buffer line or multiline
567  QgsGeometry* bufferLineClipped = 0;
568  QgsMultiPolyline mpl;
569  if ( clipBaselineBuffer->isMultipart() )
570  {
571  QgsMultiPolygon bufferMultiPolygon = clipBaselineBuffer->asMultiPolygon();
572  if ( bufferMultiPolygon.size() < 1 )
573  {
574  delete clipBaselineBuffer;
575  continue;
576  }
577 
578  for ( int j = 0; j < bufferMultiPolygon.size(); ++j )
579  {
580  int size = bufferMultiPolygon.at( j ).size();
581  for ( int k = 0; k < size; ++k )
582  {
583  mpl.append( bufferMultiPolygon.at( j ).at( k ) );
584  }
585  }
586  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
587  }
588  else
589  {
590  QgsPolygon bufferPolygon = clipBaselineBuffer->asPolygon();
591  if ( bufferPolygon.size() < 1 )
592  {
593  delete clipBaselineBuffer;
594  continue;
595  }
596 
597  int size = bufferPolygon.size();
598  for ( int j = 0; j < size; ++j )
599  {
600  mpl.append( bufferPolygon[j] );
601  }
602  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
603  }
604  bufferLineClipped = bufferLine->intersection( stratumGeom );
605  delete bufferLine;
606 
607  if ( bufferLineClipped && bufferLineClipped->type() == QGis::Line )
608  {
609  //if stratumGeom is a multipolygon, bufferLineClipped must intersect each part
610  bool bufferLineClippedIntersectsStratum = true;
611  if ( stratumGeom->wkbType() == QGis::WKBMultiPolygon || stratumGeom->wkbType() == QGis::WKBMultiPolygon25D )
612  {
613  QVector<QgsPolygon> multiPoly = stratumGeom->asMultiPolygon();
614  QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
615  for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
616  {
617  QgsGeometry* poly = QgsGeometry::fromPolygon( *multiIt );
618  if ( !poly->intersects( bufferLineClipped ) )
619  {
620  bufferLineClippedIntersectsStratum = false;
621  delete poly;
622  break;
623  }
624  delete poly;
625  }
626  }
627 
628  if ( bufferLineClippedIntersectsStratum )
629  {
630  delete clipBaselineBuffer;
631  return bufferLineClipped;
632  }
633  }
634 
635  delete bufferLineClipped;
636  delete clipBaselineBuffer;
637  currentBufferDist /= 2;
638  }
639 
640  return 0; //no solution found even with reduced tolerances
641 }
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 append(const T &value)
void setMaximum(int maximum)
double distance(const QgsGeometry &geom) const
Returns the minimum distanace between this geometry and another geometry, using GEOS.
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.
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:464
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:173
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:75
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
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:350
QList< QgsFeatureId > intersects(QgsRectangle rect) const
returns features that intersect the specified rectangle
double x() const
Definition: qgspoint.h:126
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)
void setValue(int progress)
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
QString number(int n, int base)
int toInt(bool *ok) const
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)
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
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
const T & value() const
QGis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Searches for the closest segment of geometry to the given point.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:38
iterator end()
bool isValid()
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:63
iterator begin()
static QgsGeometry * fromPoint(const QgsPoint &point)
Creates a new geometry from a QgsPoint object.
void setX(double x)
Definition: qgspoint.h:103
void setY(double y)
Definition: qgspoint.h:111
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:465
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
virtual long featureCount() const
Number of features in the layer.
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:236
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:68
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
Definition: qgspoint.h:134
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 QgsFields & pendingFields() const
returns field list in the to-be-committed state
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:73
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:74