7 #include <QProgressDialog>
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 )
24 QgsTransectSample::QgsTransectSample()
25 : mStrataLayer( NULL )
26 , mBaselineLayer( NULL )
27 , mShareBaseline( false )
28 , mMinDistanceUnits( Meters )
29 , mMinTransectLength( 0.0 )
41 if ( !mStrataLayer || !mStrataLayer->
isValid() )
46 if ( !mBaselineLayer || !mBaselineLayer->
isValid() )
52 QVariant::Type stratumIdType = QVariant::Int;
53 if ( !mStrataIdAttribute.
isEmpty() )
61 outputPointFields.
append(
QgsField(
"station_id", QVariant::Int ) );
62 outputPointFields.
append(
QgsField(
"stratum_id", stratumIdType ) );
63 outputPointFields.
append(
QgsField(
"station_code", QVariant::String ) );
64 outputPointFields.
append(
QgsField(
"start_lat", QVariant::Double ) );
65 outputPointFields.
append(
QgsField(
"start_long", QVariant::Double ) );
68 &( mStrataLayer->
crs() ) );
74 outputPointFields.
append(
QgsField(
"bearing", QVariant::Double ) );
76 &( mStrataLayer->
crs() ) );
83 usedBaselineFields.
append(
QgsField(
"stratum_id", stratumIdType ) );
86 &( mStrataLayer->
crs() ) );
93 QFileInfo outputPointInfo( mOutputPointLayer );
94 QString bufferClipLineOutput = outputPointInfo.
absolutePath() +
"/out_buffer_clip_line.shp";
102 if ( mMinDistanceUnits ==
Meters )
122 int nTotalTransects = 0;
149 QgsGeometry* baselineGeom = findBaselineGeometry( strataId.
isValid() ? strataId : -1 );
156 double minDistanceLayerUnits = minDistance;
158 double bufferDist = minDistance;
161 bufferDist = minDistance / 111319.9;
162 minDistanceLayerUnits = bufferDist;
168 delete clippedBaseline;
171 QgsGeometry* bufferLineClipped = clipBufferLine( strataGeom, clippedBaseline, bufferDist );
172 if ( !bufferLineClipped )
174 delete clippedBaseline;
187 int nCreatedTransects = 0;
189 int nMaxIterations = nTransects * 50;
194 while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
204 QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );
208 samplePointFeature.
setAttribute(
"id", nTotalTransects + 1 );
209 samplePointFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
210 samplePointFeature.
setAttribute(
"stratum_id", strataId );
212 samplePointFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
213 samplePointFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
225 double bearing = distanceArea.
bearing( sampleQgsPoint, minDistPoint ) /
M_PI * 180.0;
228 QgsPoint ptFarAway( sampleQgsPoint.
x() + ( minDistPoint.
x() - sampleQgsPoint.
x() ) * 1000000,
229 sampleQgsPoint.
y() + ( minDistPoint.
y() - sampleQgsPoint.
y() ) * 1000000 );
231 lineFarAway << sampleQgsPoint << ptFarAway;
234 if ( !lineClipStratum )
236 delete lineFarAwayGeom;
delete lineClipStratum;
241 if ( lineClipStratum->
distance( *samplePoint ) > 0.000001 )
243 delete lineFarAwayGeom;
delete lineClipStratum;
251 QgsGeometry* singleLine = closestMultilineElement( sampleQgsPoint, lineClipStratum );
254 delete lineClipStratum;
255 lineClipStratum = singleLine;
260 double transectLength = distanceArea.
measure( lineClipStratum );
261 if ( transectLength < mMinTransectLength )
263 delete lineFarAwayGeom;
delete lineClipStratum;
268 if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
270 delete lineFarAwayGeom;
delete lineClipStratum;
277 sampleLineFeature.
setAttribute(
"id", nTotalTransects + 1 );
278 sampleLineFeature.
setAttribute(
"station_id", nCreatedTransects + 1 );
279 sampleLineFeature.
setAttribute(
"stratum_id", strataId );
281 sampleLineFeature.
setAttribute(
"start_lat", latLongSamplePoint.
y() );
282 sampleLineFeature.
setAttribute(
"start_long", latLongSamplePoint.
x() );
284 outputLineWriter.
addFeature( sampleLineFeature );
288 outputPointWriter.
addFeature( samplePointFeature );
295 delete lineFarAwayGeom;
299 delete clippedBaseline;
302 bufferClipFeature.
setGeometry( bufferLineClipped );
304 bufferClipLineWriter.
addFeature( bufferClipFeature );
309 for ( ; featureMapIt != lineFeatureMap.
end(); ++featureMapIt )
311 delete( featureMapIt.
value() );
313 lineFeatureMap.
clear();
329 if ( !mBaselineLayer )
339 if ( strataId == fet.
attribute( mBaselineStrataId ) || mShareBaseline )
349 bool QgsTransectSample::otherTransectWithinDistance(
QgsGeometry* geom,
double minDistLayerUnit,
double minDistance,
QgsSpatialIndex& sIndex,
366 for ( ; lineIdIt != lineIdList.
constEnd(); ++lineIdIt )
369 if ( idMapIt != lineFeatureMap.
constEnd() )
373 closestSegmentPoints( *geom, *( idMapIt.
value() ), dist, pt1, pt2 );
376 if ( dist < minDistance )
405 if ( pl1.
size() < 2 || pl2.
size() < 2 )
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();
424 double denominatorU = v2x * v1y - v2y * v1x;
425 double denominatorT = v1x * v2y - v1y * v2x;
440 if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
442 dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
445 else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
447 dist = sqrt( d2 ); pt1 = p12; pt2 = minDistPoint2;
450 else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
452 dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
457 dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
462 double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
463 double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;
465 if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
468 pt1.
setX( p2x + u * v2x );
469 pt1.
setY( p2y + u * v2y );
495 if ( t >= 0.0 && t <= 1.0 )
500 if ( u >= 0.0 && u <= 1.0 )
506 dist = sqrt( pt1.
sqrDist( pt2 ) );
518 double minDist = DBL_MAX;
519 double currentDist = 0;
526 for ( ; it != multiPolyline.
constEnd(); ++it )
529 currentDist = pointGeom->
distance( *currentLine );
530 if ( currentDist < minDist )
532 minDist = currentDist;
533 closestLine.
reset( currentLine );
542 return closestLine.
take();
552 double currentBufferDist = tolerance;
555 for (
int i = 0; i < maxLoops; ++i )
558 QgsGeometry* clipBaselineBuffer = clippedBaseline->
buffer( currentBufferDist, 8 );
559 if ( !clipBaselineBuffer )
561 delete clipBaselineBuffer;
572 if ( bufferMultiPolygon.
size() < 1 )
574 delete clipBaselineBuffer;
578 for (
int j = 0; j < bufferMultiPolygon.
size(); ++j )
580 int size = bufferMultiPolygon.
at( j ).size();
581 for (
int k = 0; k < size; ++k )
583 mpl.
append( bufferMultiPolygon.
at( j ).at( k ) );
591 if ( bufferPolygon.
size() < 1 )
593 delete clipBaselineBuffer;
597 int size = bufferPolygon.
size();
598 for (
int j = 0; j < size; ++j )
600 mpl.
append( bufferPolygon[j] );
604 bufferLineClipped = bufferLine->
intersection( stratumGeom );
607 if ( bufferLineClipped && bufferLineClipped->
type() ==
QGis::Line )
610 bool bufferLineClippedIntersectsStratum =
true;
615 for ( ; multiIt != multiPoly.
constEnd(); ++multiIt )
620 bufferLineClippedIntersectsStratum =
false;
628 if ( bufferLineClippedIntersectsStratum )
630 delete clipBaselineBuffer;
631 return bufferLineClipped;
635 delete bufferLineClipped;
636 delete clipBaselineBuffer;
637 currentBufferDist /= 2;
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)
void append(const T &value)
void setMaximum(int maximum)
double distance(const QgsGeometry &geom) const
Returns the minimum distanace between this geometry and another geometry, using GEOS.
QgsMultiPolyline asMultiPolyline() const
Return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list.
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
QgsPolygon asPolygon() const
Return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list...
const_iterator constEnd() const
#define Q_NOWARN_DEPRECATED_PUSH
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)
QList< QgsFeatureId > intersects(QgsRectangle rect) const
returns features that intersect the specified rectangle
QgsMultiPolygon asMultiPolygon() const
Return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
void setValue(int progress)
double measure(const QgsGeometry *geometry) const
general measurement (line distance or polygon area)
void setGeometry(const QgsGeometry &geom)
Set this feature's geometry from another QgsGeometry object.
QString number(int n, int base)
int toInt(bool *ok) const
QgsTransectSample(QgsVectorLayer *strataLayer, QString strataIdAttribute, QString minDistanceAttribute, QString nPointsAttribute, DistanceUnits minDistUnits, QgsVectorLayer *baselineLayer, bool shareBaseline, QString baselineStrataId, const QString &outputPointLayer, const QString &outputLineLayer, const QString &usedBaselineLayer, double minTransectLength=0.0)
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.
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.
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
virtual long featureCount() const
Number of features in the layer.
const_iterator constBegin() const
bool insertFeature(const QgsFeature &f)
add feature to index
WriterError hasError()
checks whether there were any errors in constructor
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
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.
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 QgsFields & pendingFields() const
returns field list in the to-be-committed state
const_iterator constEnd() const
bool nextFeature(QgsFeature &f)
const_iterator constBegin() const
QString absolutePath() const
Q_DECL_DEPRECATED QgsGeometry * geometryAndOwnership()
Get the geometry object associated with this feature, and transfer ownership of the geometry to the c...
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