QGIS API Documentation  2.9.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 
142  QgsGeometry* strataGeom = fet.geometry();
143  if ( !strataGeom )
144  {
145  continue;
146  }
147 
148  //find baseline for strata
149  QVariant strataId = fet.attribute( mStrataIdAttribute );
150  QgsGeometry* baselineGeom = findBaselineGeometry( strataId.isValid() ? strataId : -1 );
151  if ( !baselineGeom )
152  {
153  continue;
154  }
155 
156  double minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
157  double minDistanceLayerUnits = minDistance;
158  //if minDistance is in meters and the data in degrees, we need to apply a rough conversion for the buffer distance
159  double bufferDist = minDistance;
160  if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QGis::DecimalDegrees )
161  {
162  bufferDist = minDistance / 111319.9;
163  minDistanceLayerUnits = bufferDist;
164  }
165 
166  QgsGeometry* clippedBaseline = strataGeom->intersection( baselineGeom );
167  if ( !clippedBaseline || clippedBaseline->wkbType() == QGis::WKBUnknown )
168  {
169  delete clippedBaseline;
170  continue;
171  }
172  QgsGeometry* bufferLineClipped = clipBufferLine( strataGeom, clippedBaseline, bufferDist );
173  if ( !bufferLineClipped )
174  {
175  delete clippedBaseline;
176  continue;
177  }
178 
179  //save clipped baseline to file
180  QgsFeature blFeature;
181  blFeature.setGeometry( *clippedBaseline );
182  blFeature.setAttribute( "stratum_id", strataId );
183  blFeature.setAttribute( "ok", "f" );
184  usedBaselineWriter.addFeature( blFeature );
185 
186  //start loop to create random points along the baseline
187  int nTransects = fet.attribute( mNPointsAttribute ).toInt();
188  int nCreatedTransects = 0;
189  int nIterations = 0;
190  int nMaxIterations = nTransects * 50;
191 
192  QgsSpatialIndex sIndex; //to check minimum distance
193  QMap< QgsFeatureId, QgsGeometry* > lineFeatureMap;
194 
195  while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
196  {
197  double randomPosition = (( double )mt_rand() / MD_RAND_MAX ) * clippedBaseline->length();
198  QgsGeometry* samplePoint = clippedBaseline->interpolate( randomPosition );
199  ++nIterations;
200  if ( !samplePoint )
201  {
202  continue;
203  }
204  QgsPoint sampleQgsPoint = samplePoint->asPoint();
205  QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );
206 
207  QgsFeature samplePointFeature;
208  samplePointFeature.setGeometry( samplePoint );
209  samplePointFeature.setAttribute( "id", nTotalTransects + 1 );
210  samplePointFeature.setAttribute( "station_id", nCreatedTransects + 1 );
211  samplePointFeature.setAttribute( "stratum_id", strataId );
212  samplePointFeature.setAttribute( "station_code", strataId.toString() + "_" + QString::number( nCreatedTransects + 1 ) );
213  samplePointFeature.setAttribute( "start_lat", latLongSamplePoint.y() );
214  samplePointFeature.setAttribute( "start_long", latLongSamplePoint.x() );
215 
216  //find closest point on clipped buffer line
217  QgsPoint minDistPoint;
218 
219  int afterVertex;
220  if ( bufferLineClipped->closestSegmentWithContext( sampleQgsPoint, minDistPoint, afterVertex ) < 0 )
221  {
222  continue;
223  }
224 
225  //bearing between sample point and min dist point (transect direction)
226  double bearing = distanceArea.bearing( sampleQgsPoint, minDistPoint ) / M_PI * 180.0;
227 
228  QgsPolyline sampleLinePolyline;
229  QgsPoint ptFarAway( sampleQgsPoint.x() + ( minDistPoint.x() - sampleQgsPoint.x() ) * 1000000,
230  sampleQgsPoint.y() + ( minDistPoint.y() - sampleQgsPoint.y() ) * 1000000 );
231  QgsPolyline lineFarAway;
232  lineFarAway << sampleQgsPoint << ptFarAway;
233  QgsGeometry* lineFarAwayGeom = QgsGeometry::fromPolyline( lineFarAway );
234  QgsGeometry* lineClipStratum = lineFarAwayGeom->intersection( strataGeom );
235  if ( !lineClipStratum )
236  {
237  delete lineFarAwayGeom; delete lineClipStratum;
238  continue;
239  }
240 
241  //cancel if distance between sample point and line is too large (line does not start at point
242  if ( lineClipStratum->distance( *samplePoint ) > 0.000001 )
243  {
244  delete lineFarAwayGeom; delete lineClipStratum;
245  continue;
246  }
247 
248  //if lineClipStratum is a multiline, take the part line closest to sampleQgsPoint
249  if ( lineClipStratum->wkbType() == QGis::WKBMultiLineString
250  || lineClipStratum->wkbType() == QGis::WKBMultiLineString25D )
251  {
252  QgsGeometry* singleLine = closestMultilineElement( sampleQgsPoint, lineClipStratum );
253  if ( singleLine )
254  {
255  delete lineClipStratum;
256  lineClipStratum = singleLine;
257  }
258  }
259 
260  //cancel if length of lineClipStratum is too small
261  double transectLength = distanceArea.measure( lineClipStratum );
262  if ( transectLength < mMinTransectLength )
263  {
264  delete lineFarAwayGeom; delete lineClipStratum;
265  continue;
266  }
267 
268  //search closest existing profile. Cancel if dist < minDist
269  if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
270  {
271  delete lineFarAwayGeom; delete lineClipStratum;
272  continue;
273  }
274 
275  QgsFeatureId fid( nCreatedTransects );
276  QgsFeature sampleLineFeature( fid );
277  sampleLineFeature.setGeometry( lineClipStratum );
278  sampleLineFeature.setAttribute( "id", nTotalTransects + 1 );
279  sampleLineFeature.setAttribute( "station_id", nCreatedTransects + 1 );
280  sampleLineFeature.setAttribute( "stratum_id", strataId );
281  sampleLineFeature.setAttribute( "station_code", strataId.toString() + "_" + QString::number( nCreatedTransects + 1 ) );
282  sampleLineFeature.setAttribute( "start_lat", latLongSamplePoint.y() );
283  sampleLineFeature.setAttribute( "start_long", latLongSamplePoint.x() );
284  sampleLineFeature.setAttribute( "bearing", bearing );
285  outputLineWriter.addFeature( sampleLineFeature );
286 
287  //add point to file writer here.
288  //It can only be written if the corresponding transect has been as well
289  outputPointWriter.addFeature( samplePointFeature );
290 
291  sIndex.insertFeature( sampleLineFeature );
292  lineFeatureMap.insert( fid, sampleLineFeature.geometryAndOwnership() );
293 
294  delete lineFarAwayGeom;
295  ++nTotalTransects;
296  ++nCreatedTransects;
297  }
298  delete clippedBaseline;
299 
300  QgsFeature bufferClipFeature;
301  bufferClipFeature.setGeometry( bufferLineClipped );
302  bufferClipFeature.setAttribute( "id", strataId );
303  bufferClipLineWriter.addFeature( bufferClipFeature );
304  //delete bufferLineClipped;
305 
306  //delete all line geometries in spatial index
307  QMap< QgsFeatureId, QgsGeometry* >::iterator featureMapIt = lineFeatureMap.begin();
308  for ( ; featureMapIt != lineFeatureMap.end(); ++featureMapIt )
309  {
310  delete( featureMapIt.value() );
311  }
312  lineFeatureMap.clear();
313  delete baselineGeom;
314 
315  ++nFeatures;
316  }
317 
318  if ( pd )
319  {
320  pd->setValue( mStrataLayer->featureCount() );
321  }
322 
323  return 0;
324 }
325 
326 QgsGeometry* QgsTransectSample::findBaselineGeometry( QVariant strataId )
327 {
328  if ( !mBaselineLayer )
329  {
330  return 0;
331  }
332 
333  QgsFeatureIterator baseLineIt = mBaselineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QStringList( mBaselineStrataId ), mBaselineLayer->pendingFields() ) );
334  QgsFeature fet;
335 
336  while ( baseLineIt.nextFeature( fet ) ) //todo: cache this in case there are many baslines
337  {
338  if ( strataId == fet.attribute( mBaselineStrataId ) || mShareBaseline )
339  {
340  return fet.geometryAndOwnership();
341  }
342  }
343  return 0;
344 }
345 
346 bool QgsTransectSample::otherTransectWithinDistance( QgsGeometry* geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex& sIndex,
347  const QMap< QgsFeatureId, QgsGeometry* >& lineFeatureMap, QgsDistanceArea& da )
348 {
349  if ( !geom )
350  {
351  return false;
352  }
353 
354  QgsGeometry* buffer = geom->buffer( minDistLayerUnit, 8 );
355  if ( !buffer )
356  {
357  return false;
358  }
359  QgsRectangle rect = buffer->boundingBox();
360  QList<QgsFeatureId> lineIdList = sIndex.intersects( rect );
361 
362  QList<QgsFeatureId>::const_iterator lineIdIt = lineIdList.constBegin();
363  for ( ; lineIdIt != lineIdList.constEnd(); ++lineIdIt )
364  {
365  const QMap< QgsFeatureId, QgsGeometry* >::const_iterator idMapIt = lineFeatureMap.find( *lineIdIt );
366  if ( idMapIt != lineFeatureMap.constEnd() )
367  {
368  double dist = 0;
369  QgsPoint pt1, pt2;
370  closestSegmentPoints( *geom, *( idMapIt.value() ), dist, pt1, pt2 );
371  dist = da.measureLine( pt1, pt2 ); //convert degrees to meters if necessary
372 
373  if ( dist < minDistance )
374  {
375  delete buffer;
376  return true;
377  }
378  }
379  }
380 
381  delete buffer;
382  return false;
383 }
384 
385 bool QgsTransectSample::closestSegmentPoints( QgsGeometry& g1, QgsGeometry& g2, double& dist, QgsPoint& pt1, QgsPoint& pt2 )
386 {
387  QGis::WkbType t1 = g1.wkbType();
388  if ( t1 != QGis::WKBLineString && t1 != QGis::WKBLineString25D )
389  {
390  return false;
391  }
392 
393  QGis::WkbType t2 = g2.wkbType();
394  if ( t2 != QGis::WKBLineString && t2 != QGis::WKBLineString25D )
395  {
396  return false;
397  }
398 
399  QgsPolyline pl1 = g1.asPolyline();
400  QgsPolyline pl2 = g2.asPolyline();
401 
402  if ( pl1.size() < 2 || pl2.size() < 2 )
403  {
404  return false;
405  }
406 
407  QgsPoint p11 = pl1.at( 0 );
408  QgsPoint p12 = pl1.at( 1 );
409  QgsPoint p21 = pl2.at( 0 );
410  QgsPoint p22 = pl2.at( 1 );
411 
412  double p1x = p11.x();
413  double p1y = p11.y();
414  double v1x = p12.x() - p11.x();
415  double v1y = p12.y() - p11.y();
416  double p2x = p21.x();
417  double p2y = p21.y();
418  double v2x = p22.x() - p21.x();
419  double v2y = p22.y() - p21.y();
420 
421  double denominatorU = v2x * v1y - v2y * v1x;
422  double denominatorT = v1x * v2y - v1y * v2x;
423 
424  if ( qgsDoubleNear( denominatorU, 0 ) || qgsDoubleNear( denominatorT, 0 ) )
425  {
426  //lines are parallel
427  //project all points on the other segment and take the one with the smallest distance
428  QgsPoint minDistPoint1;
429  double d1 = p11.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint1 );
430  QgsPoint minDistPoint2;
431  double d2 = p12.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint2 );
432  QgsPoint minDistPoint3;
433  double d3 = p21.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint3 );
434  QgsPoint minDistPoint4;
435  double d4 = p22.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint4 );
436 
437  if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
438  {
439  dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
440  return true;
441  }
442  else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
443  {
444  dist = sqrt( d2 ); pt1 = p12; pt2 = minDistPoint2;
445  return true;
446  }
447  else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
448  {
449  dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
450  return true;
451  }
452  else
453  {
454  dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
455  return true;
456  }
457  }
458 
459  double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
460  double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
461 
462  if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
463  {
464  dist = 0;
465  pt1.setX( p2x + u * v2x );
466  pt1.setY( p2y + u * v2y );
467  pt2 = pt1;
468  dist = 0;
469  return true;
470  }
471 
472  if ( t > 1.0 )
473  {
474  pt1.setX( p12.x() );
475  pt1.setY( p12.y() );
476  }
477  else if ( t < 0.0 )
478  {
479  pt1.setX( p11.x() );
480  pt1.setY( p11.y() );
481  }
482  if ( u > 1.0 )
483  {
484  pt2.setX( p22.x() );
485  pt2.setY( p22.y() );
486  }
487  if ( u < 0.0 )
488  {
489  pt2.setX( p21.x() );
490  pt2.setY( p21.y() );
491  }
492  if ( t >= 0.0 && t <= 1.0 )
493  {
494  //project pt2 onto g1
495  pt2.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), pt1 );
496  }
497  if ( u >= 0.0 && u <= 1.0 )
498  {
499  //project pt1 onto g2
500  pt1.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), pt2 );
501  }
502 
503  dist = sqrt( pt1.sqrDist( pt2 ) );
504  return true;
505 }
506 
507 QgsGeometry* QgsTransectSample::closestMultilineElement( const QgsPoint& pt, QgsGeometry* multiLine )
508 {
509  if ( !multiLine || ( multiLine->wkbType() != QGis::WKBMultiLineString
510  && multiLine->wkbType() != QGis::WKBMultiLineString25D ) )
511  {
512  return 0;
513  }
514 
515  double minDist = DBL_MAX;
516  double currentDist = 0;
517  QgsGeometry* currentLine = 0;
518  QgsGeometry* closestLine = 0;
519  QgsGeometry* pointGeom = QgsGeometry::fromPoint( pt );
520 
521  QgsMultiPolyline multiPolyline = multiLine->asMultiPolyline();
522  QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
523  for ( ; it != multiPolyline.constEnd(); ++it )
524  {
525  currentLine = QgsGeometry::fromPolyline( *it );
526  currentDist = pointGeom->distance( *currentLine );
527  if ( currentDist < minDist )
528  {
529  minDist = currentDist;
530  closestLine = currentLine;
531  }
532  else
533  {
534  delete currentLine;
535  }
536  }
537 
538  delete pointGeom;
539  return closestLine;
540 }
541 
542 QgsGeometry* QgsTransectSample::clipBufferLine( QgsGeometry* stratumGeom, QgsGeometry* clippedBaseline, double tolerance )
543 {
544  if ( !stratumGeom || !clippedBaseline || clippedBaseline->wkbType() == QGis::WKBUnknown )
545  {
546  return 0;
547  }
548 
549  double currentBufferDist = tolerance;
550  int maxLoops = 10;
551 
552  for ( int i = 0; i < maxLoops; ++i )
553  {
554  //loop with tolerance: create buffer, convert buffer to line, clip line by stratum, test if result is (single) line
555  QgsGeometry* clipBaselineBuffer = clippedBaseline->buffer( currentBufferDist, 8 );
556  if ( !clipBaselineBuffer )
557  {
558  delete clipBaselineBuffer;
559  continue;
560  }
561 
562  //it is also possible that clipBaselineBuffer is a multipolygon
563  QgsGeometry* bufferLine = 0; //buffer line or multiline
564  QgsGeometry* bufferLineClipped = 0;
565  QgsMultiPolyline mpl;
566  if ( clipBaselineBuffer->isMultipart() )
567  {
568  QgsMultiPolygon bufferMultiPolygon = clipBaselineBuffer->asMultiPolygon();
569  if ( bufferMultiPolygon.size() < 1 )
570  {
571  delete clipBaselineBuffer;
572  continue;
573  }
574 
575  for ( int j = 0; j < bufferMultiPolygon.size(); ++j )
576  {
577  int size = bufferMultiPolygon.at( j ).size();
578  for ( int k = 0; k < size; ++k )
579  {
580  mpl.append( bufferMultiPolygon.at( j ).at( k ) );
581  }
582  }
583  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
584  }
585  else
586  {
587  QgsPolygon bufferPolygon = clipBaselineBuffer->asPolygon();
588  if ( bufferPolygon.size() < 1 )
589  {
590  delete clipBaselineBuffer;
591  continue;
592  }
593 
594  int size = bufferPolygon.size();
595  for ( int j = 0; j < size; ++j )
596  {
597  mpl.append( bufferPolygon[j] );
598  }
599  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
600  }
601  bufferLineClipped = bufferLine->intersection( stratumGeom );
602  delete bufferLine;
603 
604  if ( bufferLineClipped && bufferLineClipped->type() == QGis::Line )
605  {
606  //if stratumGeom is a multipolygon, bufferLineClipped must intersect each part
607  bool bufferLineClippedIntersectsStratum = true;
608  if ( stratumGeom->wkbType() == QGis::WKBMultiPolygon || stratumGeom->wkbType() == QGis::WKBMultiPolygon25D )
609  {
610  QVector<QgsPolygon> multiPoly = stratumGeom->asMultiPolygon();
611  QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
612  for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
613  {
614  QgsGeometry* poly = QgsGeometry::fromPolygon( *multiIt );
615  if ( !poly->intersects( bufferLineClipped ) )
616  {
617  bufferLineClippedIntersectsStratum = false;
618  delete poly;
619  break;
620  }
621  delete poly;
622  }
623  }
624 
625  if ( bufferLineClippedIntersectsStratum )
626  {
627  delete clipBaselineBuffer;
628  return bufferLineClipped;
629  }
630  }
631 
632  delete bufferLineClipped;
633  delete clipBaselineBuffer;
634  currentBufferDist /= 2;
635  }
636 
637  return 0; //no solution found even with reduced tolerances
638 }
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)
double length()
get length of geometry using GEOS
const QgsField & field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
Definition: qgsfield.h:229
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.
QVector< QgsPoint > QgsPolyline
polyline is represented as a vector of points
Definition: qgsgeometry.h:33
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
Definition: qgsfeature.cpp:112
QgsPolygon asPolygon() const
return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list ...
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON)
Searches for the closest segment of geometry to the given point.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QGis::GeometryType type() const
Returns type of the vector.
Container of fields for a vector layer.
Definition: qgsfield.h:172
bool setAttribute(int field, const QVariant &attr)
Set an attribute by id.
Definition: qgsfeature.cpp:190
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:113
double sqrDist(double x, double y) const
Returns the squared distance between this point and x,y.
Definition: qgspoint.cpp:333
#define MD_RAND_MAX
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Definition: qgis.h:327
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 setGeometry(const QgsGeometry &geom)
Set this feature's geometry from another QgsGeometry object (deep copy)
Definition: qgsfeature.cpp:134
QgsGeometry * buffer(double distance, int segments)
Returns a buffer region around this geometry having the given width and with a specified number of se...
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 measure(QgsGeometry *geometry)
general measurement (line distance or polygon area)
double bearing(const QgsPoint &p1, const QgsPoint &p2)
compute bearing - in radians
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:414
#define M_PI
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QVector< QgsPolygon > QgsMultiPolygon
a collection of QgsPolygons that share a common collection of attributes
Definition: qgsgeometry.h:48
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:177
QGis::WkbType wkbType() const
Returns type of wkb (point / linestring / polygon etc.)
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:33
QVector< QgsPolyline > QgsPolygon
polygon: first item of the list is outer ring, inner rings (if any) start from second item ...
Definition: qgsgeometry.h:39
bool isValid()
QgsGeometry * intersection(QgsGeometry *geometry)
Returns a geometry representing the points shared by this geometry and other.
A class to represent a point.
Definition: qgspoint.h:63
static QgsGeometry * fromPoint(const QgsPoint &point)
construct geometry from a point
void setX(double x)
Definition: qgspoint.h:103
void setY(double y)
Definition: qgspoint.h:111
QVector< QgsPolyline > QgsMultiPolyline
a collection of QgsPolylines that share a common collection of attributes
Definition: qgsgeometry.h:45
static QgsGeometry * fromMultiPolyline(const QgsMultiPolyline &multiline)
construct geometry from a multipolyline
General purpose distance and area calculator.
QgsPolyline asPolyline() const
return contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list ...
QgsRectangle boundingBox()
Returns the bounding box of this feature.
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:230
bool isMultipart() const
Returns true if wkb of the geometry is of WKBMulti* type.
Class for storing a coordinate reference system (CRS)
double measureLine(const QList< QgsPoint > &points)
measures line
Class for doing transforms between two map coordinate systems.
static QgsGeometry * fromPolyline(const QgsPolyline &polyline)
construct geometry from a polyline
qint64 QgsFeatureId
Definition: qgsfeature.h:30
double y() const
Definition: qgspoint.h:134
const QgsCoordinateReferenceSystem & crs() const
Returns layer's spatial reference system.
static QgsGeometry * fromPolygon(const QgsPolygon &polygon)
construct geometry from a polygon
const QgsFields & pendingFields() const
returns field list in the to-be-committed state
QgsGeometry * interpolate(double distance)
bool nextFeature(QgsFeature &f)
QgsGeometry * geometryAndOwnership()
Get the geometry object associated with this feature The caller assumes responsibility for the QgsGeo...
Definition: qgsfeature.cpp:117
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)
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:64
double distance(QgsGeometry &geom)