QGIS API Documentation  2.99.0-Master (0a63d1f)
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 }
QgsPolygon asPolygon() const
Return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list...
Wrapper for iterator of features from vector data provider or vector layer.
QgsMultiPolyline asMultiPolyline() const
Return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list.
A rectangle specified with double values.
Definition: qgsrectangle.h:36
int createSample(QProgressDialog *pd)
double y
Definition: qgspoint.h:148
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)
QgsPoint asPoint() const
Return contents of the geometry as a point if wkbType is WKBPoint, otherwise returns [0...
static QgsGeometry fromPolyline(const QgsPolyline &polyline)
Creates a new geometry from a QgsPolyline object.
bool intersects(const QgsRectangle &r) const
Test for intersection with a rectangle (uses GEOS)
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
double distance(const QgsGeometry &geom) const
Returns the minimum distance between this geometry and another geometry, using GEOS.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
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:47
QList< QgsFeatureId > intersects(const QgsRectangle &rect) const
Returns features that intersect the specified rectangle.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QgsGeometry interpolate(double distance) const
Return interpolated point on line at distance.
Container of fields for a vector layer.
Definition: qgsfields.h:39
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:79
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:213
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 convenience class for writing vector files to disk.
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:136
bool isValid() const
Return the status of the layer.
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:199
#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:193
QgsPolyline asPolyline() const
Return contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list...
int mt_rand()
QgsFields fields() const
Returns the list of fields of this layer.
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:65
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:486
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
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.
#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:62
Degrees, for planar geographic CRS distance measurements.
Definition: qgsunittypes.h:51
void mt_srand(unsigned value)
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 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:59
double measureLine(const QList< QgsPoint > &points) const
Measures the length of a line with multiple segments.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:45
QVector< QgsPolyline > QgsPolygon
Polygon: first item of the list is outer ring, inner rings (if any) start from second item...
Definition: qgsgeometry.h:53
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:113
A class to represent a point.
Definition: qgspoint.h:143
QgsWkbTypes::GeometryType type() const
Returns type of the geometry as a QgsWkbTypes::GeometryType.
double length() const
Returns the length of geometry using GEOS.
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:193
void setY(double y)
Sets the y value of the point.
Definition: qgspoint.h:201
QVector< QgsPolyline > QgsMultiPolyline
A collection of QgsPolylines that share a common collection of attributes.
Definition: qgsgeometry.h:59
General purpose distance and area calculator.
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.
QgsGeometry simplify(double tolerance) const
Returns a simplified version of this geometry using a specified tolerance value.
bool insertFeature(const QgsFeature &f)
Add feature to index.
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
This class represents a coordinate reference system (CRS).
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:147
Class for doing transforms between two map coordinate systems.
qint64 QgsFeatureId
Definition: qgsfeature.h:33
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:382
bool nextFeature(QgsFeature &f)
long srsid() const
Returns the internal CRS ID, if available.
Represents a vector layer which manages a vector based data sets.
QVariant::Type type() const
Gets variant type of the field as it will be retrieved from data source.
Definition: qgsfield.cpp:94
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:262
QgsField field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
Definition: qgsfields.cpp:140
static QgsGeometry fromMultiPolyline(const QgsMultiPolyline &multiline)
Creates a new geometry from a QgsMultiPolyline object.
double x
Definition: qgspoint.h:147