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