7 #include <QProgressDialog>
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 )
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 )
44 if ( !mStrataLayer || !mStrataLayer->
isValid() )
49 if ( !mBaselineLayer || !mBaselineLayer->
isValid() )
55 QVariant::Type stratumIdType = QVariant::Int;
56 if ( !mStrataIdAttribute.
isEmpty() )
58 stratumIdType = mStrataLayer->
fields().
field( mStrataIdAttribute ).
type();
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 ) );
71 &( mStrataLayer->
crs() ) );
77 outputPointFields.
append(
QgsField(
"bearing", QVariant::Double ) );
79 &( mStrataLayer->
crs() ) );
86 usedBaselineFields.
append(
QgsField(
"stratum_id", stratumIdType ) );
89 &( mStrataLayer->
crs() ) );
96 QFileInfo outputPointInfo( mOutputPointLayer );
97 QString bufferClipLineOutput = outputPointInfo.
absolutePath() +
"/out_buffer_clip_line.shp";
105 if ( mMinDistanceUnits ==
Meters )
125 int nTotalTransects = 0;
152 QgsGeometry* baselineGeom = findBaselineGeometry( strataId.
isValid() ? strataId : -1 );
159 double minDistanceLayerUnits = minDistance;
161 double bufferDist = bufferDistance( minDistance );
164 minDistanceLayerUnits = minDistance / 111319.9;
170 delete clippedBaseline;
173 QgsGeometry* bufferLineClipped = clipBufferLine( strataGeom, clippedBaseline, bufferDist );
174 if ( !bufferLineClipped )
176 delete clippedBaseline;
189 int nCreatedTransects = 0;
191 int nMaxIterations = nTransects * 50;
196 while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
206 QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );
210 samplePointFeature.
setAttribute(
"id", nTotalTransects + 1 );
211 samplePointFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
212 samplePointFeature.
setAttribute(
"stratum_id", strataId );
214 samplePointFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
215 samplePointFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
227 double bearing = distanceArea.
bearing( sampleQgsPoint, minDistPoint ) /
M_PI * 180.0;
230 QgsPoint ptFarAway( sampleQgsPoint.
x() + ( minDistPoint.
x() - sampleQgsPoint.
x() ) * 1000000,
231 sampleQgsPoint.
y() + ( minDistPoint.
y() - sampleQgsPoint.
y() ) * 1000000 );
233 lineFarAway << sampleQgsPoint << ptFarAway;
236 if ( !lineClipStratum )
238 delete lineFarAwayGeom;
delete lineClipStratum;
243 if ( lineClipStratum->
distance( *samplePoint ) > 0.000001 )
245 delete lineFarAwayGeom;
delete lineClipStratum;
253 QgsGeometry* singleLine = closestMultilineElement( sampleQgsPoint, lineClipStratum );
256 delete lineClipStratum;
257 lineClipStratum = singleLine;
262 double transectLength = distanceArea.
measureLength( lineClipStratum );
263 if ( transectLength < mMinTransectLength )
265 delete lineFarAwayGeom;
delete lineClipStratum;
270 if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
272 delete lineFarAwayGeom;
delete lineClipStratum;
279 sampleLineFeature.
setAttribute(
"id", nTotalTransects + 1 );
280 sampleLineFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
281 sampleLineFeature.
setAttribute(
"stratum_id", strataId );
283 sampleLineFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
284 sampleLineFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
286 outputLineWriter.
addFeature( sampleLineFeature );
290 outputPointWriter.
addFeature( samplePointFeature );
297 delete lineFarAwayGeom;
301 delete clippedBaseline;
304 bufferClipFeature.
setGeometry( bufferLineClipped );
306 bufferClipLineWriter.
addFeature( bufferClipFeature );
311 for ( ; featureMapIt != lineFeatureMap.
end(); ++featureMapIt )
313 delete( featureMapIt.
value() );
315 lineFeatureMap.
clear();
331 if ( !mBaselineLayer )
341 if ( strataId == fet.
attribute( mBaselineStrataId ) || mShareBaseline )
351 bool QgsTransectSample::otherTransectWithinDistance(
QgsGeometry* geom,
double minDistLayerUnit,
double minDistance,
QgsSpatialIndex& sIndex,
368 for ( ; lineIdIt != lineIdList.
constEnd(); ++lineIdIt )
371 if ( idMapIt != lineFeatureMap.
constEnd() )
375 closestSegmentPoints( *geom, *( idMapIt.
value() ), dist, pt1, pt2 );
378 if ( dist < minDistance )
407 if ( pl1.
size() < 2 || pl2.
size() < 2 )
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();
426 double denominatorU = v2x * v1y - v2y * v1x;
427 double denominatorT = v1x * v2y - v1y * v2x;
442 if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
444 dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
447 else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
449 dist = sqrt( d2 ); pt1 = p12; pt2 = minDistPoint2;
452 else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
454 dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
459 dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
464 double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
465 double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
467 if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
470 pt1.
setX( p2x + u * v2x );
471 pt1.
setY( p2y + u * v2y );
497 if ( t >= 0.0 && t <= 1.0 )
502 if ( u >= 0.0 && u <= 1.0 )
508 dist = sqrt( pt1.
sqrDist( pt2 ) );
520 double minDist = DBL_MAX;
521 double currentDist = 0;
528 for ( ; it != multiPolyline.
constEnd(); ++it )
531 currentDist = pointGeom->
distance( *currentLine );
532 if ( currentDist < minDist )
534 minDist = currentDist;
535 closestLine.
reset( currentLine );
544 return closestLine.
take();
555 if ( mBaselineSimplificationTolerance >= 0 )
558 usedBaseline = clippedBaseline->
simplify( mBaselineSimplificationTolerance );
571 double currentBufferDist = tolerance;
574 for (
int i = 0; i < maxLoops; ++i )
578 if ( !clipBaselineBuffer )
580 delete clipBaselineBuffer;
591 if ( bufferMultiPolygon.
size() < 1 )
593 delete clipBaselineBuffer;
597 for (
int j = 0; j < bufferMultiPolygon.
size(); ++j )
599 int size = bufferMultiPolygon.
at( j ).size();
600 for (
int k = 0; k < size; ++k )
602 mpl.
append( bufferMultiPolygon.
at( j ).at( k ) );
610 if ( bufferPolygon.
size() < 1 )
612 delete clipBaselineBuffer;
616 int size = bufferPolygon.
size();
618 for (
int j = 0; j < size; ++j )
620 mpl.
append( bufferPolygon[j] );
624 bufferLineClipped = bufferLine->
intersection( stratumGeom );
627 if ( bufferLineClipped && bufferLineClipped->
type() ==
QGis::Line )
630 bool bufferLineClippedIntersectsStratum =
true;
635 for ( ; multiIt != multiPoly.
constEnd(); ++multiIt )
640 bufferLineClippedIntersectsStratum =
false;
648 if ( bufferLineClippedIntersectsStratum )
650 delete clipBaselineBuffer;
651 if ( mBaselineSimplificationTolerance >= 0 )
655 return bufferLineClipped;
659 delete bufferLineClipped;
660 delete clipBaselineBuffer;
661 currentBufferDist /= 2;
664 if ( mBaselineSimplificationTolerance >= 0 )
671 double QgsTransectSample::bufferDistance(
double minDistanceFromAttribute )
const
673 double bufferDist = minDistanceFromAttribute;
674 if ( mBaselineBufferDistance >= 0 )
676 bufferDist = mBaselineBufferDistance;
681 bufferDist /= 111319.9;
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.
int createSample(QProgressDialog *pd)
const QgsField & field(int fieldIdx) const
Get field at particular index (must be in range 0..N-1)
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)
QList< QgsFeatureId > intersects(const QgsRectangle &rect) const
Returns features that intersect the specified rectangle.
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.
QgsFields fields() const
Returns the list of fields of this layer.
long srsid() const
Get the SrsId - if possible.
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
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.
A geometry is the spatial representation of a feature.
bool setAttribute(int field, const QVariant &attr)
Set an attribute's value by field index.
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.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
double x() const
Get the x value of the point.
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
long featureCount(QgsSymbolV2 *symbol)
Number of features rendered with specified symbol.
void setValue(int progress)
void setGeometry(const QgsGeometry &geom)
Set this feature's geometry from another QgsGeometry object.
QString number(int n, int base)
int toInt(bool *ok) 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.
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) ...
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.
bool isValid()
Return the status of the layer.
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.
static QgsGeometry * fromPoint(const QgsPoint &point)
Creates a new geometry from a QgsPoint object.
void setX(double x)
Sets the x value of the point.
void setY(double y)
Sets the y value of the point.
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
General purpose distance and area calculator.
QgsGeometry * intersection(const QgsGeometry *geometry) const
Returns a geometry representing the points shared by this geometry and other.
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
double measureLength(const QgsGeometry *geometry) const
Measures the length of a geometry.
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.
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.
static QgsGeometry * fromPolyline(const QgsPolyline &polyline)
Creates a new geometry from a QgsPolyline object.
double y() const
Get the y value of the point.
const QgsCoordinateReferenceSystem & crs() const
Returns layer's spatial reference system.
double toDouble(bool *ok) const
iterator insert(const Key &key, const T &value)
static QgsGeometry * fromPolygon(const QgsPolygon &polygon)
Creates a new geometry from a QgsPolygon.
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...
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.
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.
QGis::UnitType mapUnits() const
Get the units that the projection is in.