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