7 #include <QProgressDialog>
16 QgsVectorLayer* baselineLayer,
bool shareBaseline, QString baselineStrataId,
const QString& outputPointLayer,
17 const QString& outputLineLayer,
const QString& usedBaselineLayer,
double minTransectLength ): mStrataLayer( strataLayer ),
18 mStrataIdAttribute( strataIdAttribute ), mMinDistanceAttribute( minDistanceAttribute ), mNPointsAttribute( nPointsAttribute ), mBaselineLayer( baselineLayer ), mShareBaseline( shareBaseline ),
19 mBaselineStrataId( baselineStrataId ), mOutputPointLayer( outputPointLayer ), mOutputLineLayer( outputLineLayer ), mUsedBaselineLayer( usedBaselineLayer ),
20 mMinDistanceUnits( minDistUnits ), mMinTransectLength( minTransectLength )
47 QVariant::Type stratumIdType = QVariant::Int;
56 outputPointFields.
append(
QgsField(
"station_id", QVariant::Int ) );
57 outputPointFields.
append(
QgsField(
"stratum_id", stratumIdType ) );
58 outputPointFields.
append(
QgsField(
"station_code", QVariant::String ) );
59 outputPointFields.
append(
QgsField(
"start_lat", QVariant::Double ) );
60 outputPointFields.
append(
QgsField(
"start_long", QVariant::Double ) );
69 outputPointFields.
append(
QgsField(
"bearing", QVariant::Double ) );
78 usedBaselineFields.
append(
QgsField(
"stratum_id", stratumIdType ) );
89 QString bufferClipLineOutput = outputPointInfo.absolutePath() +
"/out_buffer_clip_line.shp";
110 mt_srand( QTime::currentTime().msec() );
117 int nTotalTransects = 0;
129 pd->setValue( nFeatures );
131 if ( pd && pd->wasCanceled() )
144 bool strataIdOk =
true;
153 double minDistanceLayerUnits = minDistance;
155 double bufferDist = minDistance;
158 bufferDist = minDistance / 111319.9;
159 minDistanceLayerUnits = bufferDist;
165 delete clippedBaseline;
169 if ( !bufferLineClipped )
171 delete clippedBaseline;
184 int nCreatedTransects = 0;
186 int nMaxIterations = nTransects * 50;
189 QMap< QgsFeatureId, QgsGeometry* > lineFeatureMap;
191 while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
201 QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );
205 samplePointFeature.
setAttribute(
"id", nTotalTransects + 1 );
206 samplePointFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
207 samplePointFeature.
setAttribute(
"stratum_id", strataId );
208 samplePointFeature.
setAttribute(
"station_code", strataId.toString() +
"_" + QString::number( nCreatedTransects + 1 ) );
209 samplePointFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
210 samplePointFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
222 double bearing = distanceArea.
bearing( sampleQgsPoint, minDistPoint ) /
M_PI * 180.0;
225 QgsPoint ptFarAway( sampleQgsPoint.
x() + ( minDistPoint.
x() - sampleQgsPoint.
x() ) * 1000000,
226 sampleQgsPoint.
y() + ( minDistPoint.
y() - sampleQgsPoint.
y() ) * 1000000 );
228 lineFarAway << sampleQgsPoint << ptFarAway;
231 if ( !lineClipStratum )
233 delete lineFarAwayGeom;
delete lineClipStratum;
238 if ( lineClipStratum->
distance( *samplePoint ) > 0.000001 )
240 delete lineFarAwayGeom;
delete lineClipStratum;
251 delete lineClipStratum;
252 lineClipStratum = singleLine;
257 double transectLength = distanceArea.
measure( lineClipStratum );
260 delete lineFarAwayGeom;
delete lineClipStratum;
267 delete lineFarAwayGeom;
delete lineClipStratum;
274 sampleLineFeature.
setAttribute(
"id", nTotalTransects + 1 );
275 sampleLineFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
276 sampleLineFeature.
setAttribute(
"stratum_id", strataId );
277 sampleLineFeature.
setAttribute(
"station_code", strataId.toString() +
"_" + QString::number( nCreatedTransects + 1 ) );
278 sampleLineFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
279 sampleLineFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
281 outputLineWriter.
addFeature( sampleLineFeature );
285 outputPointWriter.
addFeature( samplePointFeature );
290 delete lineFarAwayGeom;
294 delete clippedBaseline;
297 bufferClipFeature.
setGeometry( bufferLineClipped );
299 bufferClipLineWriter.
addFeature( bufferClipFeature );
303 QMap< QgsFeatureId, QgsGeometry* >::iterator featureMapIt = lineFeatureMap.begin();
304 for ( ; featureMapIt != lineFeatureMap.end(); ++featureMapIt )
306 delete( featureMapIt.value() );
308 lineFeatureMap.clear();
343 const QMap< QgsFeatureId, QgsGeometry* >& lineFeatureMap,
QgsDistanceArea& da )
356 QList<QgsFeatureId> lineIdList = sIndex.
intersects( rect );
358 QList<QgsFeatureId>::const_iterator lineIdIt = lineIdList.constBegin();
359 for ( ; lineIdIt != lineIdList.constEnd(); ++lineIdIt )
361 const QMap< QgsFeatureId, QgsGeometry* >::const_iterator idMapIt = lineFeatureMap.find( *lineIdIt );
362 if ( idMapIt != lineFeatureMap.constEnd() )
369 if ( dist < minDistance )
398 if ( pl1.size() < 2 || pl2.size() < 2 )
408 double p1x = p11.
x();
409 double p1y = p11.
y();
410 double v1x = p12.
x() - p11.
x();
411 double v1y = p12.
y() - p11.
y();
412 double p2x = p21.
x();
413 double p2y = p21.
y();
414 double v2x = p22.
x() - p21.
x();
415 double v2y = p22.
y() - p21.
y();
417 double denominatorU = v2x * v1y - v2y * v1x;
418 double denominatorT = v1x * v2y - v1y * v2x;
433 if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
435 dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
438 else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
440 dist = sqrt( d2 ); pt1 = p12; pt2 = minDistPoint2;
443 else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
445 dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
450 dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
455 double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
456 double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
458 if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
461 pt1.
setX( p2x + u * v2x );
462 pt1.
setY( p2y + u * v2y );
488 if ( t >= 0.0 && t <= 1.0 )
493 if ( u >= 0.0 && u <= 1.0 )
499 dist = sqrt( pt1.
sqrDist( pt2 ) );
511 double minDist = DBL_MAX;
512 double currentDist = 0;
518 QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
519 for ( ; it != multiPolyline.constEnd(); ++it )
522 currentDist = pointGeom->
distance( *currentLine );
523 if ( currentDist < minDist )
525 minDist = currentDist;
526 closestLine = currentLine;
545 double currentBufferDist = tolerance;
548 for (
int i = 0; i < maxLoops; ++i )
551 QgsGeometry* clipBaselineBuffer = clippedBaseline->
buffer( currentBufferDist, 8 );
552 if ( !clipBaselineBuffer )
554 delete clipBaselineBuffer;
565 if ( bufferMultiPolygon.size() < 1 )
567 delete clipBaselineBuffer;
571 for (
int j = 0; j < bufferMultiPolygon.size(); ++j )
573 int size = bufferMultiPolygon.at( j ).size();
574 for (
int k = 0; k <
size; ++k )
576 mpl.append( bufferMultiPolygon.at( j ).at( k ) );
584 if ( bufferPolygon.size() < 1 )
586 delete clipBaselineBuffer;
590 int size = bufferPolygon.size();
591 for (
int j = 0; j <
size; ++j )
593 mpl.append( bufferPolygon[j] );
597 bufferLineClipped = bufferLine->
intersection( stratumGeom );
599 if ( bufferLineClipped && bufferLineClipped->
type() ==
QGis::Line )
602 bool bufferLineClippedIntersectsStratum =
true;
606 QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
607 for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
612 bufferLineClippedIntersectsStratum =
false;
620 if ( bufferLineClippedIntersectsStratum )
622 return bufferLineClipped;
626 delete bufferLineClipped;
delete clipBaselineBuffer;
delete bufferLine;
627 currentBufferDist /= 2;
static bool otherTransectWithinDistance(QgsGeometry *geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex &sIndex, const QMap< QgsFeatureId, QgsGeometry * > &lineFeatureMap, QgsDistanceArea &da)
Returns true if another transect is within the specified minimum distance.
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
int createSample(QProgressDialog *pd)
double length()
get length of geometry using GEOS
bool isMultipart()
Returns true if wkb of the geometry is of WKBMulti* type.
const QgsField & field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
QgsMultiPolyline asMultiPolyline() const
return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
QVector< QgsPoint > QgsPolyline
polyline is represented as a vector of points
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
QgsPolygon asPolygon() const
return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list ...
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON)
Searches for the closest segment of geometry to the given point.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QGis::GeometryType type()
Returns type of the vector.
Container of fields for a vector layer.
bool setAttribute(int field, const QVariant &attr)
Set an attribute by id.
WkbType
Used for symbology operations.
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...
double sqrDist(double x, double y) const
Returns the squared distance between this point and x,y.
DistanceUnits mMinDistanceUnits
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
QList< QgsFeatureId > intersects(QgsRectangle rect) const
returns features that intersect the specified rectangle
QString mBaselineStrataId
QgsMultiPolygon asMultiPolygon() const
return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
static bool closestSegmentPoints(QgsGeometry &g1, QgsGeometry &g2, double &dist, QgsPoint &pt1, QgsPoint &pt2)
Finds the closest points between two line segments.
void setGeometry(const QgsGeometry &geom)
Set this feature's geometry from another QgsGeometry object (deep copy)
QgsGeometry * buffer(double distance, int segments)
Returns a buffer region around this geometry having the given width and with a specified number of se...
double measure(QgsGeometry *geometry)
general measurement (line distance or polygon area)
QString mOutputPointLayer
double bearing(const QgsPoint &p1, const QgsPoint &p2)
compute bearing - in radians
double sqrDistToSegment(double x1, double y1, double x2, double y2, QgsPoint &minDistPoint, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Returns the minimum distance between this point and a segment.
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
void mt_srand(unsigned value)
QgsGeometry * findBaselineGeometry(QVariant strataId)
bool addFeature(QgsFeature &feature, QgsFeatureRendererV2 *renderer=0, QGis::UnitType outputUnit=QGis::Meters)
add feature to the currently opened shapefile
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Append a field. The field must have unique name, otherwise it is rejected (returns false) ...
QString mStrataIdAttribute
QGis::WkbType wkbType() const
Returns type of wkb (point / linestring / polygon etc.)
static QgsGeometry * closestMultilineElement(const QgsPoint &pt, QgsGeometry *multiLine)
Returns a copy of the multiline element closest to a point (caller takes ownership) ...
Encapsulate a field in an attribute table or data source.
QVector< QgsPolyline > QgsPolygon
polygon: first item of the list is outer ring, inner rings (if any) start from second item ...
QgsGeometry * intersection(QgsGeometry *geometry)
Returns a geometry representing the points shared by this geometry and other.
A class to represent a point geometry.
static QgsGeometry * fromPoint(const QgsPoint &point)
construct geometry from a point
QVector< QgsPolyline > QgsMultiPolyline
a collection of QgsPolylines that share a common collection of attributes
static QgsGeometry * fromMultiPolyline(const QgsMultiPolyline &multiline)
construct geometry from a multipolyline
QString mMinDistanceAttribute
General purpose distance and area calculator.
QgsPolyline asPolyline() const
return contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list ...
QgsVectorLayer * mBaselineLayer
QgsRectangle boundingBox()
Returns the bounding box of this feature.
virtual long featureCount() const
Number of features in the layer.
bool insertFeature(const QgsFeature &f)
add feature to index
WriterError hasError()
checks whether there were any errors in constructor
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
double mMinTransectLength
QString mNPointsAttribute
Class for storing a coordinate reference system (CRS)
double measureLine(const QList< QgsPoint > &points)
measures line
static QgsGeometry * fromPolyline(const QgsPolyline &polyline)
construct geometry from a polyline
const QgsCoordinateReferenceSystem & crs() const
Returns layer's spatial reference system.
QString mUsedBaselineLayer
static QgsGeometry * fromPolygon(const QgsPolygon &polygon)
construct geometry from a polygon
const QgsFields & pendingFields() const
returns field list in the to-be-committed state
QgsGeometry * interpolate(double distance)
QgsVectorLayer * mStrataLayer
static QgsGeometry * clipBufferLine(QgsGeometry *stratumGeom, QgsGeometry *clippedBaseline, double tolerance)
Returns clipped buffer line.
bool nextFeature(QgsFeature &f)
QgsGeometry * geometryAndOwnership()
Get the geometry object associated with this feature The caller assumes responsibility for the QgsGeo...
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.
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.
double distance(QgsGeometry &geom)
QGis::UnitType mapUnits() const