QGIS API Documentation  2.99.0-Master (c558d51)
qgstransectsample.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstransectsample.cpp
3  ---------------------
4  begin : July 2013
5  copyright : (C) 2013 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 #include "qgstransectsample.h"
16 #include "qgsdistancearea.h"
17 #include "qgsfeatureiterator.h"
18 #include "qgsgeometry.h"
19 #include "qgsspatialindex.h"
20 #include "qgsvectorfilewriter.h"
21 #include "qgsvectorlayer.h"
22 #include <QProgressDialog>
23 #include <QFileInfo>
24 #ifndef _MSC_VER
25 #include <stdint.h>
26 #endif
27 #include "mersenne-twister.h"
28 #include <limits>
29 
30 QgsTransectSample::QgsTransectSample( QgsVectorLayer* strataLayer, const QString& strataIdAttribute, const QString& minDistanceAttribute, const QString& nPointsAttribute, DistanceUnits minDistUnits,
31  QgsVectorLayer* baselineLayer, bool shareBaseline, const QString& baselineStrataId, const QString& outputPointLayer,
32  const QString& outputLineLayer, const QString& usedBaselineLayer, double minTransectLength,
33  double baselineBufferDistance, double baselineSimplificationTolerance )
34  : mStrataLayer( strataLayer )
35  , mStrataIdAttribute( strataIdAttribute )
36  , mMinDistanceAttribute( minDistanceAttribute )
37  , mNPointsAttribute( nPointsAttribute )
38  , mBaselineLayer( baselineLayer )
39  , mShareBaseline( shareBaseline )
40  , mBaselineStrataId( baselineStrataId )
41  , mOutputPointLayer( outputPointLayer )
42  , mOutputLineLayer( outputLineLayer )
43  , mUsedBaselineLayer( usedBaselineLayer )
44  , mMinDistanceUnits( minDistUnits )
45  , mMinTransectLength( minTransectLength )
46  , mBaselineBufferDistance( baselineBufferDistance )
47  , mBaselineSimplificationTolerance( baselineSimplificationTolerance )
48 {
49 }
50 
52  : mStrataLayer( nullptr )
53  , mBaselineLayer( nullptr )
54  , mShareBaseline( false )
55  , mMinDistanceUnits( Meters )
56  , mMinTransectLength( 0.0 )
57  , mBaselineBufferDistance( -1.0 )
58  , mBaselineSimplificationTolerance( -1.0 )
59 {
60 }
61 
62 int QgsTransectSample::createSample( QProgressDialog* pd )
63 {
64  Q_UNUSED( pd );
65 
66  if ( !mStrataLayer || !mStrataLayer->isValid() )
67  {
68  return 1;
69  }
70 
71  if ( !mBaselineLayer || !mBaselineLayer->isValid() )
72  {
73  return 2;
74  }
75 
76  //stratum id is not necessarily an integer
77  QVariant::Type stratumIdType = QVariant::Int;
78  if ( !mStrataIdAttribute.isEmpty() )
79  {
80  stratumIdType = mStrataLayer->fields().field( mStrataIdAttribute ).type();
81  }
82 
83  //create vector file writers for output
84  QgsFields outputPointFields;
85  outputPointFields.append( QgsField( QStringLiteral( "id" ), stratumIdType ) );
86  outputPointFields.append( QgsField( QStringLiteral( "station_id" ), QVariant::Int ) );
87  outputPointFields.append( QgsField( QStringLiteral( "stratum_id" ), stratumIdType ) );
88  outputPointFields.append( QgsField( QStringLiteral( "station_code" ), QVariant::String ) );
89  outputPointFields.append( QgsField( QStringLiteral( "start_lat" ), QVariant::Double ) );
90  outputPointFields.append( QgsField( QStringLiteral( "start_long" ), QVariant::Double ) );
91 
92  QgsVectorFileWriter outputPointWriter( mOutputPointLayer, QStringLiteral( "utf-8" ), outputPointFields, QgsWkbTypes::Point,
93  mStrataLayer->crs() );
94  if ( outputPointWriter.hasError() != QgsVectorFileWriter::NoError )
95  {
96  return 3;
97  }
98 
99  outputPointFields.append( QgsField( QStringLiteral( "bearing" ), QVariant::Double ) ); //add bearing attribute for lines
100  QgsVectorFileWriter outputLineWriter( mOutputLineLayer, QStringLiteral( "utf-8" ), outputPointFields, QgsWkbTypes::LineString,
101  mStrataLayer->crs() );
102  if ( outputLineWriter.hasError() != QgsVectorFileWriter::NoError )
103  {
104  return 4;
105  }
106 
107  QgsFields usedBaselineFields;
108  usedBaselineFields.append( QgsField( QStringLiteral( "stratum_id" ), stratumIdType ) );
109  usedBaselineFields.append( QgsField( QStringLiteral( "ok" ), QVariant::String ) );
110  QgsVectorFileWriter usedBaselineWriter( mUsedBaselineLayer, QStringLiteral( "utf-8" ), usedBaselineFields, QgsWkbTypes::LineString,
111  mStrataLayer->crs() );
112  if ( usedBaselineWriter.hasError() != QgsVectorFileWriter::NoError )
113  {
114  return 5;
115  }
116 
117  //debug: write clipped buffer bounds with stratum id to same directory as out_point
118  QFileInfo outputPointInfo( mOutputPointLayer );
119  QString bufferClipLineOutput = outputPointInfo.absolutePath() + "/out_buffer_clip_line.shp";
120  QgsFields bufferClipLineFields;
121  bufferClipLineFields.append( QgsField( QStringLiteral( "id" ), stratumIdType ) );
122  QgsVectorFileWriter bufferClipLineWriter( bufferClipLineOutput, QStringLiteral( "utf-8" ), bufferClipLineFields, QgsWkbTypes::LineString, mStrataLayer->crs() );
123 
124  //configure distanceArea depending on minDistance units and output CRS
125  QgsDistanceArea distanceArea;
126  distanceArea.setSourceCrs( mStrataLayer->crs().srsid() );
127  if ( mMinDistanceUnits == Meters )
128  {
129  distanceArea.setEllipsoidalMode( true );
130  }
131  else
132  {
133  distanceArea.setEllipsoidalMode( false );
134  }
135 
136  //possibility to transform output points to lat/long
138 
139  //init random number generator
140  mt_srand( QTime::currentTime().msec() );
141 
143  fr.setSubsetOfAttributes( QStringList() << mStrataIdAttribute << mMinDistanceAttribute << mNPointsAttribute, mStrataLayer->fields() );
144  QgsFeatureIterator strataIt = mStrataLayer->getFeatures( fr );
145 
146  QgsFeature fet;
147  int nTotalTransects = 0;
148  int nFeatures = 0;
149 
150  if ( pd )
151  {
152  pd->setMaximum( mStrataLayer->featureCount() );
153  }
154 
155  while ( strataIt.nextFeature( fet ) )
156  {
157  if ( pd )
158  {
159  pd->setValue( nFeatures );
160  }
161  if ( pd && pd->wasCanceled() )
162  {
163  break;
164  }
165 
166  if ( !fet.hasGeometry() )
167  {
168  continue;
169  }
170  QgsGeometry strataGeom = fet.geometry();
171 
172  //find baseline for strata
173  QVariant strataId = fet.attribute( mStrataIdAttribute );
174  QgsGeometry baselineGeom = findBaselineGeometry( strataId.isValid() ? strataId : -1 );
175  if ( baselineGeom.isEmpty() )
176  {
177  continue;
178  }
179 
180  double minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
181  double minDistanceLayerUnits = minDistance;
182  //if minDistance is in meters and the data in degrees, we need to apply a rough conversion for the buffer distance
183  double bufferDist = bufferDistance( minDistance );
184  if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QgsUnitTypes::DistanceDegrees )
185  {
186  minDistanceLayerUnits = minDistance / 111319.9;
187  }
188 
189  QgsGeometry clippedBaseline = strataGeom.intersection( baselineGeom );
190  if ( !clippedBaseline || clippedBaseline.wkbType() == QgsWkbTypes::Unknown )
191  {
192  continue;
193  }
194  QgsGeometry* bufferLineClipped = clipBufferLine( strataGeom, &clippedBaseline, bufferDist );
195  if ( !bufferLineClipped )
196  {
197  continue;
198  }
199 
200  //save clipped baseline to file
201  QgsFeature blFeature( usedBaselineFields );
202  blFeature.setGeometry( clippedBaseline );
203  blFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
204  blFeature.setAttribute( QStringLiteral( "ok" ), "f" );
205  usedBaselineWriter.addFeature( blFeature );
206 
207  //start loop to create random points along the baseline
208  int nTransects = fet.attribute( mNPointsAttribute ).toInt();
209  int nCreatedTransects = 0;
210  int nIterations = 0;
211  int nMaxIterations = nTransects * 50;
212 
213  QgsSpatialIndex sIndex; //to check minimum distance
214  QMap< QgsFeatureId, QgsGeometry > lineFeatureMap;
215 
216  while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
217  {
218  double randomPosition = (( double )mt_rand() / MD_RAND_MAX ) * clippedBaseline.length();
219  QgsGeometry samplePoint = clippedBaseline.interpolate( randomPosition );
220  ++nIterations;
221  if ( samplePoint.isEmpty() )
222  {
223  continue;
224  }
225  QgsPoint sampleQgsPoint = samplePoint.asPoint();
226  QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );
227 
228  QgsFeature samplePointFeature( outputPointFields );
229  samplePointFeature.setGeometry( samplePoint );
230  samplePointFeature.setAttribute( QStringLiteral( "id" ), nTotalTransects + 1 );
231  samplePointFeature.setAttribute( QStringLiteral( "station_id" ), nCreatedTransects + 1 );
232  samplePointFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
233  samplePointFeature.setAttribute( QStringLiteral( "station_code" ), strataId.toString() + '_' + QString::number( nCreatedTransects + 1 ) );
234  samplePointFeature.setAttribute( QStringLiteral( "start_lat" ), latLongSamplePoint.y() );
235  samplePointFeature.setAttribute( QStringLiteral( "start_long" ), latLongSamplePoint.x() );
236 
237  //find closest point on clipped buffer line
238  QgsPoint minDistPoint;
239 
240  int afterVertex;
241  if ( bufferLineClipped->closestSegmentWithContext( sampleQgsPoint, minDistPoint, afterVertex ) < 0 )
242  {
243  continue;
244  }
245 
246  //bearing between sample point and min dist point (transect direction)
247  double bearing = distanceArea.bearing( sampleQgsPoint, minDistPoint ) / M_PI * 180.0;
248 
249  QgsPoint ptFarAway( sampleQgsPoint.x() + ( minDistPoint.x() - sampleQgsPoint.x() ) * 1000000,
250  sampleQgsPoint.y() + ( minDistPoint.y() - sampleQgsPoint.y() ) * 1000000 );
251  QgsPolyline lineFarAway;
252  lineFarAway << sampleQgsPoint << ptFarAway;
253  QgsGeometry lineFarAwayGeom = QgsGeometry::fromPolyline( lineFarAway );
254  QgsGeometry lineClipStratum = lineFarAwayGeom.intersection( strataGeom );
255  if ( lineClipStratum.isEmpty() )
256  {
257  continue;
258  }
259 
260  //cancel if distance between sample point and line is too large (line does not start at point
261  if ( lineClipStratum.distance( samplePoint ) > 0.000001 )
262  {
263  continue;
264  }
265 
266  //if lineClipStratum is a multiline, take the part line closest to sampleQgsPoint
267  if ( lineClipStratum.wkbType() == QgsWkbTypes::MultiLineString
268  || lineClipStratum.wkbType() == QgsWkbTypes::MultiLineString25D )
269  {
270  QgsGeometry singleLine = closestMultilineElement( sampleQgsPoint, lineClipStratum );
271  if ( !singleLine.isEmpty() )
272  {
273  lineClipStratum = singleLine;
274  }
275  }
276 
277  //cancel if length of lineClipStratum is too small
278  double transectLength = distanceArea.measureLength( lineClipStratum );
279  if ( transectLength < mMinTransectLength )
280  {
281  continue;
282  }
283 
284  //search closest existing profile. Cancel if dist < minDist
285  if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
286  {
287  continue;
288  }
289 
290  QgsFeatureId fid( nCreatedTransects );
291  QgsFeature sampleLineFeature( outputPointFields, fid );
292  sampleLineFeature.setGeometry( lineClipStratum );
293  sampleLineFeature.setAttribute( QStringLiteral( "id" ), nTotalTransects + 1 );
294  sampleLineFeature.setAttribute( QStringLiteral( "station_id" ), nCreatedTransects + 1 );
295  sampleLineFeature.setAttribute( QStringLiteral( "stratum_id" ), strataId );
296  sampleLineFeature.setAttribute( QStringLiteral( "station_code" ), strataId.toString() + '_' + QString::number( nCreatedTransects + 1 ) );
297  sampleLineFeature.setAttribute( QStringLiteral( "start_lat" ), latLongSamplePoint.y() );
298  sampleLineFeature.setAttribute( QStringLiteral( "start_long" ), latLongSamplePoint.x() );
299  sampleLineFeature.setAttribute( QStringLiteral( "bearing" ), bearing );
300  outputLineWriter.addFeature( sampleLineFeature );
301 
302  //add point to file writer here.
303  //It can only be written if the corresponding transect has been as well
304  outputPointWriter.addFeature( samplePointFeature );
305 
306  sIndex.insertFeature( sampleLineFeature );
307  lineFeatureMap.insert( fid, lineClipStratum );
308 
309  ++nTotalTransects;
310  ++nCreatedTransects;
311  }
312 
313  QgsFeature bufferClipFeature;
314  bufferClipFeature.setGeometry( *bufferLineClipped );
315  delete bufferLineClipped;
316  bufferClipFeature.setAttribute( QStringLiteral( "id" ), strataId );
317  bufferClipLineWriter.addFeature( bufferClipFeature );
318  //delete bufferLineClipped;
319 
320  ++nFeatures;
321  }
322 
323  if ( pd )
324  {
325  pd->setValue( mStrataLayer->featureCount() );
326  }
327 
328  return 0;
329 }
330 
331 QgsGeometry QgsTransectSample::findBaselineGeometry( const QVariant& strataId )
332 {
333  if ( !mBaselineLayer )
334  {
335  return QgsGeometry();
336  }
337 
338  QgsFeatureIterator baseLineIt = mBaselineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QStringList( mBaselineStrataId ), mBaselineLayer->fields() ) );
339  QgsFeature fet;
340 
341  while ( baseLineIt.nextFeature( fet ) ) //todo: cache this in case there are many baslines
342  {
343  if ( strataId == fet.attribute( mBaselineStrataId ) || mShareBaseline )
344  {
345  return fet.geometry();
346  }
347  }
348  return QgsGeometry();
349 }
350 
351 bool QgsTransectSample::otherTransectWithinDistance( const QgsGeometry& geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex& sIndex,
352  const QMap< QgsFeatureId, QgsGeometry >& lineFeatureMap, QgsDistanceArea& da )
353 {
354  if ( geom.isEmpty() )
355  {
356  return false;
357  }
358 
359  QgsGeometry buffer = geom.buffer( minDistLayerUnit, 8 );
360  if ( buffer.isEmpty() )
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  return true;
381  }
382  }
383  }
384 
385  return false;
386 }
387 
388 bool QgsTransectSample::closestSegmentPoints( const QgsGeometry& g1, const QgsGeometry& g2, double& dist, QgsPoint& pt1, QgsPoint& pt2 )
389 {
390  QgsWkbTypes::Type t1 = g1.wkbType();
392  {
393  return false;
394  }
395 
396  QgsWkbTypes::Type t2 = g2.wkbType();
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 );
443  pt1 = p11;
444  pt2 = minDistPoint1;
445  return true;
446  }
447  else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
448  {
449  dist = sqrt( d2 );
450  pt1 = p12;
451  pt2 = minDistPoint2;
452  return true;
453  }
454  else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
455  {
456  dist = sqrt( d3 );
457  pt1 = p21;
458  pt2 = minDistPoint3;
459  return true;
460  }
461  else
462  {
463  dist = sqrt( d4 );
464  pt1 = p21;
465  pt2 = minDistPoint4;
466  return true;
467  }
468  }
469 
470  double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
471  double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
472 
473  if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
474  {
475  dist = 0;
476  pt1.setX( p2x + u * v2x );
477  pt1.setY( p2y + u * v2y );
478  pt2 = pt1;
479  dist = 0;
480  return true;
481  }
482 
483  if ( t > 1.0 )
484  {
485  pt1.setX( p12.x() );
486  pt1.setY( p12.y() );
487  }
488  else if ( t < 0.0 )
489  {
490  pt1.setX( p11.x() );
491  pt1.setY( p11.y() );
492  }
493  if ( u > 1.0 )
494  {
495  pt2.setX( p22.x() );
496  pt2.setY( p22.y() );
497  }
498  if ( u < 0.0 )
499  {
500  pt2.setX( p21.x() );
501  pt2.setY( p21.y() );
502  }
503  if ( t >= 0.0 && t <= 1.0 )
504  {
505  //project pt2 onto g1
506  pt2.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), pt1 );
507  }
508  if ( u >= 0.0 && u <= 1.0 )
509  {
510  //project pt1 onto g2
511  pt1.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), pt2 );
512  }
513 
514  dist = sqrt( pt1.sqrDist( pt2 ) );
515  return true;
516 }
517 
518 QgsGeometry QgsTransectSample::closestMultilineElement( const QgsPoint& pt, const QgsGeometry& multiLine )
519 {
520  if ( !multiLine || ( multiLine.wkbType() != QgsWkbTypes::MultiLineString
521  && multiLine.wkbType() != QgsWkbTypes::MultiLineString25D ) )
522  {
523  return QgsGeometry();
524  }
525 
526  double minDist = DBL_MAX;
527  double currentDist = 0;
528  QgsGeometry currentLine;
529  QgsGeometry closestLine;
530  QgsGeometry pointGeom = QgsGeometry::fromPoint( pt );
531 
532  QgsMultiPolyline multiPolyline = multiLine.asMultiPolyline();
533  QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
534  for ( ; it != multiPolyline.constEnd(); ++it )
535  {
536  currentLine = QgsGeometry::fromPolyline( *it );
537  currentDist = pointGeom.distance( currentLine );
538  if ( currentDist < minDist )
539  {
540  minDist = currentDist;
541  closestLine = currentLine;
542  }
543  }
544 
545  return closestLine;
546 }
547 
548 QgsGeometry* QgsTransectSample::clipBufferLine( const QgsGeometry& stratumGeom, QgsGeometry* clippedBaseline, double tolerance )
549 {
550  if ( !stratumGeom || !clippedBaseline || clippedBaseline->wkbType() == QgsWkbTypes::Unknown )
551  {
552  return nullptr;
553  }
554 
555  QgsGeometry usedBaseline = *clippedBaseline;
556  if ( mBaselineSimplificationTolerance >= 0 )
557  {
558  //int verticesBefore = usedBaseline->asMultiPolyline().count();
559  usedBaseline = clippedBaseline->simplify( mBaselineSimplificationTolerance );
560  if ( usedBaseline.isEmpty() )
561  {
562  return nullptr;
563  }
564  //int verticesAfter = usedBaseline->asMultiPolyline().count();
565 
566  //debug: write to file
567  /*QgsVectorFileWriter debugWriter( "/tmp/debug.shp", "utf-8", QgsFields(), QgsWkbTypes::LineString, &( mStrataLayer->crs() ) );
568  QgsFeature debugFeature; debugFeature.setGeometry( usedBaseline );
569  debugWriter.addFeature( debugFeature );*/
570  }
571 
572  double currentBufferDist = tolerance;
573  int maxLoops = 10;
574 
575  for ( int i = 0; i < maxLoops; ++i )
576  {
577  //loop with tolerance: create buffer, convert buffer to line, clip line by stratum, test if result is (single) line
578  QgsGeometry clipBaselineBuffer = usedBaseline.buffer( currentBufferDist, 8 );
579  if ( clipBaselineBuffer.isEmpty() )
580  {
581  continue;
582  }
583 
584  //it is also possible that clipBaselineBuffer is a multipolygon
585  QgsGeometry bufferLine; //buffer line or multiline
586  QgsGeometry bufferLineClipped;
587  QgsMultiPolyline mpl;
588  if ( clipBaselineBuffer.isMultipart() )
589  {
590  QgsMultiPolygon bufferMultiPolygon = clipBaselineBuffer.asMultiPolygon();
591  if ( bufferMultiPolygon.size() < 1 )
592  {
593  continue;
594  }
595 
596  for ( int j = 0; j < bufferMultiPolygon.size(); ++j )
597  {
598  int size = bufferMultiPolygon.at( j ).size();
599  for ( int k = 0; k < size; ++k )
600  {
601  mpl.append( bufferMultiPolygon.at( j ).at( k ) );
602  }
603  }
604  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
605  }
606  else
607  {
608  QgsPolygon bufferPolygon = clipBaselineBuffer.asPolygon();
609  if ( bufferPolygon.size() < 1 )
610  {
611  continue;
612  }
613 
614  int size = bufferPolygon.size();
615  mpl.reserve( size );
616  for ( int j = 0; j < size; ++j )
617  {
618  mpl.append( bufferPolygon[j] );
619  }
620  bufferLine = QgsGeometry::fromMultiPolyline( mpl );
621  }
622  bufferLineClipped = bufferLine.intersection( stratumGeom );
623 
624  if ( bufferLineClipped.isEmpty() && bufferLineClipped.type() == QgsWkbTypes::LineGeometry )
625  {
626  //if stratumGeom is a multipolygon, bufferLineClipped must intersect each part
627  bool bufferLineClippedIntersectsStratum = true;
628  if ( stratumGeom.wkbType() == QgsWkbTypes::MultiPolygon || stratumGeom.wkbType() == QgsWkbTypes::MultiPolygon25D )
629  {
630  QVector<QgsPolygon> multiPoly = stratumGeom.asMultiPolygon();
631  QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
632  for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
633  {
634  QgsGeometry poly = QgsGeometry::fromPolygon( *multiIt );
635  if ( !poly.intersects( bufferLineClipped ) )
636  {
637  bufferLineClippedIntersectsStratum = false;
638  break;
639  }
640  }
641  }
642 
643  if ( bufferLineClippedIntersectsStratum )
644  {
645  return new QgsGeometry( bufferLineClipped );
646  }
647  }
648 
649  currentBufferDist /= 2;
650  }
651 
652  return nullptr; //no solution found even with reduced tolerances
653 }
654 
655 double QgsTransectSample::bufferDistance( double minDistanceFromAttribute ) const
656 {
657  double bufferDist = minDistanceFromAttribute;
658  if ( mBaselineBufferDistance >= 0 )
659  {
660  bufferDist = mBaselineBufferDistance;
661  }
662 
663  if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QgsUnitTypes::DistanceDegrees )
664  {
665  bufferDist /= 111319.9;
666  }
667 
668  return bufferDist;
669 }
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 y
Definition: qgspoint.h:116
QgsGeometry interpolate(double distance) const
Return interpolated point on line at distance.
QgsTransectSample(QgsVectorLayer *strataLayer, const QString &strataIdAttribute, const QString &minDistanceAttribute, const QString &nPointsAttribute, DistanceUnits minDistUnits, QgsVectorLayer *baselineLayer, bool shareBaseline, const QString &baselineStrataId, const QString &outputPointLayer, const QString &outputLineLayer, const QString &usedBaselineLayer, double minTransectLength=0.0, double baselineBufferDistance=-1.0, double baselineSimplificationTolerance=-1.0)
bool isValid() const
Return the status of the layer.
static QgsGeometry fromPolyline(const QgsPolyline &polyline)
Creates a new geometry from a QgsPolyline object.
QList< QgsFeatureId > intersects(const QgsRectangle &rect) const
Returns features that intersect the specified rectangle.
double distance(const QgsGeometry &geom) const
Returns the minimum distance between this geometry and another geometry, using GEOS.
QgsFields fields() const
Returns the list of fields of this layer.
long srsid() const
Returns the internal CRS ID, if available.
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)
QVector< QgsPoint > QgsPolyline
Polyline is represented as a vector of points.
Definition: qgsgeometry.h:46
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=nullptr, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Searches for the closest segment of geometry to the given point.
QgsPolygon asPolygon() const
Return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list...
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
Container of fields for a vector layer.
Definition: qgsfields.h:36
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:78
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:228
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:135
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...
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:348
#define MD_RAND_MAX
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:196
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
int mt_rand()
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const
Query the layer for features specified in request.
static QgsGeometry fromPoint(const QgsPoint &point)
Creates a new geometry from a QgsPoint object.
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:452
#define M_PI
This class wraps a request for features to a vector layer (or directly its vector data provider)...
double measureLine(const QList< QgsPoint > &points) const
Measures the length of a line with multiple segments.
QVector< QgsPolygon > QgsMultiPolygon
A collection of QgsPolygons that share a common collection of attributes.
Definition: qgsgeometry.h:61
Degrees, for planar geographic CRS distance measurements.
Definition: qgsunittypes.h:50
void mt_srand(unsigned value)
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: qgsfields.cpp:61
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:47
QVector< QgsPolyline > QgsPolygon
Polygon: first item of the list is outer ring, inner rings (if any) start from second item...
Definition: qgsgeometry.h:52
QgsField field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
Definition: qgsfields.cpp:142
A class to represent a point.
Definition: qgspoint.h:111
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:113
static QgsGeometry fromPolygon(const QgsPolygon &polygon)
Creates a new geometry from a QgsPolygon.
void setX(double x)
Sets the x value of the point.
Definition: qgspoint.h:161
void setY(double y)
Sets the y value of the point.
Definition: qgspoint.h:169
double length() const
Returns the length of geometry using GEOS.
QVector< QgsPolyline > QgsMultiPolyline
A collection of QgsPolylines that share a common collection of attributes.
Definition: qgsgeometry.h:58
QgsGeometry simplify(double tolerance) const
Returns a simplified version of this geometry using a specified tolerance value.
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...
bool insertFeature(const QgsFeature &f)
Add feature to index.
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:277
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
This class represents a coordinate reference system (CRS).
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:162
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:214
Class for doing transforms between two map coordinate systems.
qint64 QgsFeatureId
Definition: qgsfeature.h:32
bool isEmpty() const
Returns true if the geometry is empty (ie, contains no underlying geometry accessible via geometry)...
QgsCoordinateReferenceSystem crs() const
Returns the layer&#39;s spatial reference system.
bool nextFeature(QgsFeature &f)
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.
QgsWkbTypes::GeometryType type() const
Returns type of the geometry as a QgsWkbTypes::GeometryType.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
QVariant::Type type() const
Gets variant type of the field as it will be retrieved from data source.
Definition: qgsfield.cpp:98
static QgsGeometry fromMultiPolyline(const QgsMultiPolyline &multiline)
Creates a new geometry from a QgsMultiPolyline object.
double x
Definition: qgspoint.h:115