Fix crashes due to geometry pointers in QgsGeometryAnalyzer; Add QgsOverlayAnalyzer to analysis library (currently only supports intersections)

git-svn-id: http://svn.osgeo.org/qgis/trunk@12040 c8812cc2-4d05-0410-92ff-de0c093fc19c
This commit is contained in:
cfarmer 2009-11-08 23:10:48 +00:00
parent 3c8107349b
commit c1b360bd93
5 changed files with 449 additions and 56 deletions

View File

@ -8,4 +8,5 @@
%Import core/core.sip
%Include qgsgeometryanalyzer.sip
%Include qgsoverlayanalyzer.sip

View File

@ -31,6 +31,7 @@ SET(QGIS_ANALYSIS_SRCS
raster/qgstotalcurvaturefilter.cpp
vector/qgsgeometryanalyzer.cpp
vector/qgszonalstatistics.cpp
vector/qgsoverlayanalyzer.cpp
)
SET(QGIS_ANALYSIS_MOC_HDRS
@ -43,6 +44,7 @@ INCLUDE_DIRECTORIES(
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/../core/
${CMAKE_CURRENT_SOURCE_DIR}/../core/renderer
${CMAKE_CURRENT_SOURCE_DIR}/../core/spatialindex
interpolation
${PROJ_INCLUDE_DIR}
${GEOS_INCLUDE_DIR}
@ -99,7 +101,7 @@ INSTALL(TARGETS qgis_analysis
# Added by Tim to install headers
SET(QGIS_ANALYSIS_HDRS vector/qgsgeometryanalyzer.h vector/qgszonalstatistics.h
SET(QGIS_ANALYSIS_HDRS vector/qgsgeometryanalyzer.h vector/qgszonalstatistics.h vector/qgsgeometryanalyzer.h
)
INSTALL(CODE "MESSAGE(\"Installing ANALYSIS headers...\")")

View File

@ -329,6 +329,7 @@ bool QgsGeometryAnalyzer::extent( QgsVectorLayer* layer, const QString& shapefil
QList<double> QgsGeometryAnalyzer::simpleMeasure( QgsGeometry* mpGeometry )
{
QList<double> list;
double perim;
if ( mpGeometry->wkbType() == QGis::WKBPoint )
{
QgsPoint pt = mpGeometry->asPoint();
@ -341,11 +342,11 @@ QList<double> QgsGeometryAnalyzer::simpleMeasure( QgsGeometry* mpGeometry )
list.append( measure.measure( mpGeometry ) );
if ( mpGeometry->type() == QGis::Polygon )
{
list.append( perimeterMeasure( mpGeometry, measure ) );
perim = perimeterMeasure( mpGeometry, measure );
list.append( perim );
}
}
return list;
}
double QgsGeometryAnalyzer::perimeterMeasure( QgsGeometry* geometry, QgsDistanceArea& measure )
@ -383,22 +384,12 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
{
return false;
}
QgsVectorDataProvider* dp = layer->dataProvider();
if ( !dp )
{
return false;
}
QGis::WkbType outputType = QGis::WKBPolygon;
const QgsCoordinateReferenceSystem crs = layer->srs();
QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), dp->fields(), outputType, &crs );
QgsFeature currentFeature;
QgsGeometry* dissolveGeometry; //dissolve geometry
QMultiMap<QString, int> map;
bool useField = false;
if ( uniqueIdField == -1 )
{
uniqueIdField = 0;
@ -406,14 +397,63 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
else
{
useField = true;
}
QgsFieldMap fields;
fields.insert( 0 , QgsField( QString( "UID" ), QVariant::String ) );
fields.insert( 1 , QgsField( QString( "AREA" ), QVariant::Double ) );
fields.insert( 2 , QgsField( QString( "PERIM" ), QVariant::Double ) );
QGis::WkbType outputType = QGis::WKBPolygon;
const QgsCoordinateReferenceSystem crs = layer->srs();
QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), fields, outputType, &crs );
QgsFeature currentFeature;
QgsGeometry* dissolveGeometry; //dissolve geometry
QMultiMap<QString, int> map;
if ( onlySelectedFeatures )
{
//use QgsVectorLayer::featureAtId
const QgsFeatureIds selection = layer->selectedFeaturesIds();
QgsFeatureIds::const_iterator it = selection.constBegin();
for ( ; it != selection.constEnd(); ++it )
{
// if ( p )
// {
// p->setValue( processedFeatures );
// }
// if ( p && p->wasCanceled() )
// {
// // break; // it may be better to do something else here?
// return false;
// }
if ( !layer->featureAtId( *it, currentFeature, true, true ) )
{
continue;
}
map.insert( currentFeature.attributeMap()[ uniqueIdField ].toString(), currentFeature.id() );
}
}
else
{
layer->select( layer->pendingAllAttributesList(), QgsRectangle(), true, false );
while ( layer->nextFeature( currentFeature ) )
{
map.insert( currentFeature.attributeMap()[uniqueIdField].toString(), currentFeature.id() );
// if ( p )
// {
// p->setValue( processedFeatures );
// }
// if ( p && p->wasCanceled() )
// {
// // break; // it may be better to do something else here?
// return false;
// }
map.insert( currentFeature.attributeMap()[ uniqueIdField ].toString(), currentFeature.id() );
}
}
QMultiMap<QString, int>::const_iterator jt;
for (jt = map.constBegin(); jt != map.constEnd(); ++jt)
QMultiMap<QString, int>::const_iterator jt = map.constBegin();
while ( jt != map.constEnd() )
{
QString currentKey = jt.key();
int processedFeatures = 0;
@ -421,13 +461,13 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
if ( onlySelectedFeatures )
{
//use QgsVectorLayer::featureAtId
const QgsFeatureIds selection = layer->selectedFeaturesIds();
const QgsFeatureIds selection = layer->selectedFeaturesIds();
if ( p )
{
p->setMaximum( selection.size() );
p->setMaximum( selection.size() );
}
processedFeatures = 0;
while ( jt != map.end() && ( jt.key() == currentKey || !useField ) )
while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
{
if ( p && p->wasCanceled() )
{
@ -437,18 +477,26 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
{
if ( p )
{
p->setValue( processedFeatures );
p->setValue( processedFeatures );
}
if ( !layer->featureAtId( jt.value(), currentFeature, true, true ) )
{
continue;
continue;
}
convexFeature( currentFeature, processedFeatures, &dissolveGeometry );
++processedFeatures;
}
++jt;
}
QList<double> values;
dissolveGeometry = dissolveGeometry->convexHull();
values = simpleMeasure( dissolveGeometry );
QgsAttributeMap attributeMap;
attributeMap.insert( 0 , QVariant( currentKey ) );
attributeMap.insert( 1 , values[ 0 ] );
attributeMap.insert( 2 , values[ 1 ] );
QgsFeature dissolveFeature;
dissolveFeature.setAttributeMap( attributeMap );
dissolveFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( dissolveFeature );
}
@ -461,7 +509,7 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
p->setMaximum( featureCount );
}
processedFeatures = 0;
while ( jt != map.end() && ( jt.key() == currentKey || !useField ) )
while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
{
if ( p )
{
@ -480,9 +528,19 @@ bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shap
++processedFeatures;
++jt;
}
QgsFeature dissolveFeature;
dissolveFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( dissolveFeature );
QList<double> values;
// QgsGeometry* tmpGeometry = 0;
dissolveGeometry = dissolveGeometry->convexHull();
// values = simpleMeasure( tmpGeometry );
values = simpleMeasure( dissolveGeometry );
QgsAttributeMap attributeMap;
attributeMap.insert( 0 , QVariant( currentKey ) );
attributeMap.insert( 1 , QVariant( values[ 0 ] ) );
attributeMap.insert( 2 , QVariant( values[ 1 ] ) );
QgsFeature dissolveFeature;
dissolveFeature.setAttributeMap( attributeMap );
dissolveFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( dissolveFeature );
}
}
return true;
@ -522,22 +580,12 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
{
return false;
}
QgsVectorDataProvider* dp = layer->dataProvider();
if ( !dp )
{
return false;
}
QGis::WkbType outputType = dp->geometryType();
const QgsCoordinateReferenceSystem crs = layer->srs();
QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), dp->fields(), outputType, &crs );
QgsFeature currentFeature;
QgsGeometry* dissolveGeometry; //dissolve geometry
QMultiMap<QString, int> map;
bool useField = false;
if ( uniqueIdField == -1 )
{
uniqueIdField = 0;
@ -545,28 +593,56 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
else
{
useField = true;
}
QGis::WkbType outputType = dp->geometryType();
const QgsCoordinateReferenceSystem crs = layer->srs();
QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), dp->fields(), outputType, &crs );
QgsFeature currentFeature;
QMultiMap<QString, int> map;
if ( onlySelectedFeatures )
{
//use QgsVectorLayer::featureAtId
const QgsFeatureIds selection = layer->selectedFeaturesIds();
QgsFeatureIds::const_iterator it = selection.constBegin();
for ( ; it != selection.constEnd(); ++it )
{
if ( !layer->featureAtId( *it, currentFeature, true, true ) )
{
continue;
}
map.insert( currentFeature.attributeMap()[ uniqueIdField ].toString(), currentFeature.id() );
}
}
else
{
layer->select( layer->pendingAllAttributesList(), QgsRectangle(), true, false );
while ( layer->nextFeature( currentFeature ) )
{
map.insert( currentFeature.attributeMap()[uniqueIdField].toString(), currentFeature.id() );
map.insert( currentFeature.attributeMap()[ uniqueIdField ].toString(), currentFeature.id() );
}
}
QMultiMap<QString, int>::const_iterator jt;
for (jt = map.constBegin(); jt != map.constEnd(); ++jt)
QgsGeometry* dissolveGeometry; //dissolve geometry
QMultiMap<QString, int>::const_iterator jt = map.constBegin();
QgsFeature outputFeature;
while ( jt != map.constEnd() )
{
QString currentKey = jt.key();
int processedFeatures = 0;
bool first = true;
//take only selection
if ( onlySelectedFeatures )
{
//use QgsVectorLayer::featureAtId
const QgsFeatureIds selection = layer->selectedFeaturesIds();
const QgsFeatureIds selection = layer->selectedFeaturesIds();
if ( p )
{
p->setMaximum( selection.size() );
p->setMaximum( selection.size() );
}
processedFeatures = 0;
while ( jt != map.end() && ( jt.key() == currentKey || !useField ) )
while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
{
if ( p && p->wasCanceled() )
{
@ -576,20 +652,22 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
{
if ( p )
{
p->setValue( processedFeatures );
p->setValue( processedFeatures );
}
if ( !layer->featureAtId( jt.value(), currentFeature, true, true ) )
{
continue;
continue;
}
if ( first )
{
outputFeature.setAttributeMap( currentFeature.attributeMap() );
first = false;
}
dissolveFeature( currentFeature, processedFeatures, &dissolveGeometry );
++processedFeatures;
}
++jt;
}
QgsFeature dissolveFeature;
dissolveFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( dissolveFeature );
}
//take all features
else
@ -599,8 +677,7 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
{
p->setMaximum( featureCount );
}
processedFeatures = 0;
while ( jt != map.end() && ( jt.key() == currentKey || !useField ) )
while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
{
if ( p )
{
@ -615,14 +692,17 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
{
continue;
}
{
outputFeature.setAttributeMap( currentFeature.attributeMap() );
first = false;
}
dissolveFeature( currentFeature, processedFeatures, &dissolveGeometry );
++processedFeatures;
++jt;
}
QgsFeature dissolveFeature;
dissolveFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( dissolveFeature );
}
outputFeature.setGeometry( dissolveGeometry );
vWriter.addFeature( outputFeature );
}
return true;
}
@ -630,7 +710,6 @@ bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapef
void QgsGeometryAnalyzer::dissolveFeature( QgsFeature& f, int nProcessedFeatures, QgsGeometry** dissolveGeometry )
{
QgsGeometry* featureGeometry = f.geometry();
QgsGeometry* tmpGeometry = 0;
if ( !featureGeometry )
{
@ -639,13 +718,15 @@ void QgsGeometryAnalyzer::dissolveFeature( QgsFeature& f, int nProcessedFeatures
if ( nProcessedFeatures == 0 )
{
*dissolveGeometry = featureGeometry;
int geomSize = featureGeometry->wkbSize();
*dissolveGeometry = new QgsGeometry();
unsigned char* wkb = new unsigned char[geomSize];
memcpy(wkb, featureGeometry->asWkb(), geomSize);
(*dissolveGeometry)->fromWkb(wkb, geomSize);
}
else
{
tmpGeometry = *dissolveGeometry;
*dissolveGeometry = ( *dissolveGeometry )->combine( featureGeometry );
delete tmpGeometry;
}
}

View File

@ -0,0 +1,204 @@
/***************************************************************************
qgsoverlayanalyzer.cpp - QGIS Tools for vector geometry analysis
-------------------
begin : 8 Nov 2009
copyright : (C) Carson J. Q. Farmer
email : carson.farmer@gmail.com
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
/* $Id: qgis.h 9774 2008-12-12 05:41:24Z timlinux $ */
#include "qgsoverlayanalyzer.h"
#include "qgsapplication.h"
#include "qgsfield.h"
#include "qgsfeature.h"
#include "qgslogger.h"
#include "qgscoordinatereferencesystem.h"
#include "qgsvectorfilewriter.h"
#include "qgsvectordataprovider.h"
#include "qgsdistancearea.h"
#include <QProgressDialog>
bool QgsOverlayAnalyzer::intersection( QgsVectorLayer* layerA, QgsVectorLayer* layerB,
const QString& shapefileName, bool onlySelectedFeatures,
QProgressDialog* p )
{
if ( !layerA && !layerB )
{
return false;
}
QgsVectorDataProvider* dpA = layerA->dataProvider();
QgsVectorDataProvider* dpB = layerB->dataProvider();
if ( !dpA && !dpB )
{
return false;
}
QGis::WkbType outputType = dpA->geometryType();
const QgsCoordinateReferenceSystem crs = layerA->srs();
QgsFieldMap fields;
//fields = combineFieldLists( dpA->fields(), dpB->fields() );
QgsVectorFileWriter vWriter( shapefileName, dpA->encoding(), fields, outputType, &crs );
QgsFeature currentFeature;
QgsGeometry* dissolveGeometry; //dissolve geometry (if dissolve enabled)
QgsSpatialIndex index;
//take only selection
if ( onlySelectedFeatures )
{
const QgsFeatureIds selectionB = layerB->selectedFeaturesIds();
QgsFeatureIds::const_iterator it = selectionB.constBegin();
for ( ; it != selectionB.constEnd(); ++it )
{
if ( !layerB->featureAtId( *it, currentFeature, true, true ) )
{
continue;
}
index.insertFeature( currentFeature );
}
//use QgsVectorLayer::featureAtId
const QgsFeatureIds selectionA = layerA->selectedFeaturesIds();
if ( p )
{
p->setMaximum( selectionA.size() );
}
int processedFeatures = 0;
it = selectionA.constBegin();
for ( ; it != selectionA.constEnd(); ++it )
{
if ( p )
{
p->setValue( processedFeatures );
}
if ( p && p->wasCanceled() )
{
break;
}
if ( !layerA->featureAtId( *it, currentFeature, true, true ) )
{
continue;
}
intersectFeature( currentFeature, &vWriter, layerB, &index );
++processedFeatures;
}
if ( p )
{
p->setValue( selectionA.size() );
}
}
//take all features
else
{
layerB->select( layerB->pendingAllAttributesList(), QgsRectangle(), true, false );
while ( dpB->nextFeature( currentFeature ) )
{
index.insertFeature( currentFeature );
}
layerA->select( layerA->pendingAllAttributesList(), QgsRectangle(), true, false );
int featureCount = layerA->featureCount();
if ( p )
{
p->setMaximum( featureCount );
}
int processedFeatures = 0;
while ( layerA->nextFeature( currentFeature ) )
{
if ( p )
{
p->setValue( processedFeatures );
}
if ( p && p->wasCanceled() )
{
break;
}
intersectFeature( currentFeature, &vWriter, layerB, &index );
++processedFeatures;
}
if ( p )
{
p->setValue( featureCount );
}
}
return true;
}
void QgsOverlayAnalyzer::intersectFeature( QgsFeature& f, QgsVectorFileWriter* vfw,
QgsVectorLayer* vl, QgsSpatialIndex* index )
{
QgsGeometry* featureGeometry = f.geometry();
QgsGeometry* intersectGeometry = 0;
QgsFeature currentFeature;
if ( !featureGeometry )
{
return;
}
QList<int> intersects;
intersects = index->intersects( featureGeometry->boundingBox() );
QList<int>::const_iterator it = intersects.constBegin();
for ( ; it != intersects.constEnd(); ++it )
{
if ( !vl->featureAtId( *it, currentFeature, true, true ) )
{
continue;
}
if ( featureGeometry->intersects( currentFeature.geometry() ) )
{
intersectGeometry = featureGeometry->intersection( currentFeature.geometry() );
QgsFeature outFeature;
outFeature.setGeometry( intersectGeometry );
outFeature.setAttributeMap( f.attributeMap() );
//add it to vector file writer
if ( vfw )
{
vfw->addFeature( outFeature );
}
}
}
}
void QgsOverlayAnalyzer::combineFieldLists( QgsFieldMap fieldListA, QgsFieldMap fieldListB )
{
QMap<int, QgsField>::const_iterator i = fieldListB.constBegin();
int count = 0;
while ( i != fieldListB.constEnd() )
{
if ( !fieldListA.contains( i.key() ) )
{
fieldListA.insert( fieldListA.size()-1, i.value() );
count = 0;
}
else
{
QgsField field;
field = i.value();
QString name = field.name();
name.append( "_" ).append( QString( count ) );
fieldListA.insert( fieldListA.size()-1 , QgsField( name, field.type() ) );
++count;
continue;
}
++i;
}
}

View File

@ -0,0 +1,105 @@
/***************************************************************************
qgsoverlayanalyzer.h - QGIS Tools for vector geometry analysis
-------------------
begin : 19 March 2009
copyright : (C) Carson Farmer
email : carson.farmer@gmail.com
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
/* $Id: qgis.h 9774 2008-12-12 05:41:24Z timlinux $ */
#ifndef QGSOVERLAYANALYZERH
#define QGSOVERLAYANALYZERH
#include "qgsvectorlayer.h"
#include "qgsfield.h"
#include "qgsspatialindex.h"
#include "qgsfeature.h"
#include "qgsgeometry.h"
#include "qgsfield.h"
#include "qgsdistancearea.h"
class QgsVectorFileWriter;
class QProgressDialog;
/** \ingroup analysis
* The QGis class provides vector overlay analysis functions
*/
class ANALYSIS_EXPORT QgsOverlayAnalyzer
{
public:
/**Perform an intersection on two input vector layers and write output to a new shape file
@param layerA input vector layer
@param layerB input vector layer
@param shapefileName path to the output shp
@param onlySelectedFeatures if true, only selected features are considered, else all the features
@param p progress dialog (or 0 if no progress dialog is to be shown)
@note: added in version 1.4*/
bool intersection( QgsVectorLayer* layerA, QgsVectorLayer* layerB, \
const QString& shapefileName, bool onlySelectedFeatures = false, \
QProgressDialog* p = 0 );
// /**Perform a union of two input vector layers and write output to a new shape file
// @param layerA input vector layer
// @param layerB input vector layer
// @param shapefileName path to the output shp
// @param onlySelectedFeatures if true, only selected features are considered, else all the features
// @param p progress dialog (or 0 if no progress dialog is to be shown)
// @note: added in version 1.4*/
// bool combine( QgsVectorLayer* layerA, QgsVectorLayer* layerB,
// const QString& shapefileName, bool onlySelectedFeatures = false,
// QProgressDialog* p = 0 );
//
// /**Clip a vector layer based on the boundary of another vector layer and
// write output to a new shape file
// @param layerA input vector layer
// @param layerB input vector layer
// @param shapefileName path to the output shp
// @param onlySelectedFeatures if true, only selected features are considered, else all the features
// @param p progress dialog (or 0 if no progress dialog is to be shown)
// @note: added in version 1.4*/
// bool clip( QgsVectorLayer* layerA, QgsVectorLayer* layerB,
// const QString& shapefileName, bool onlySelectedFeatures = false,
// QProgressDialog* p = 0 );
//
// /**Difference a vector layer based on the geometries of another vector layer
// and write the output to a new shape file
// @param layerA input vector layer
// @param layerB input vector layer
// @param shapefileName path to the output shp
// @param onlySelectedFeatures if true, only selected features are considered, else all the features
// @param p progress dialog (or 0 if no progress dialog is to be shown)
// @note: added in version 1.4*/
// bool difference( QgsVectorLayer* layerA, QgsVectorLayer* layerB,
// const QString& shapefileName, bool onlySelectedFeatures = false,
// QProgressDialog* p = 0 );
//
// /**Intersect two vector layers and write the geometries of each layer that
// do not intersect with the other layer to a new shape file (Symmetrical difference)
// @param layerA input vector layer
// @param layerB input vector layer
// @param shapefileName path to the output shp
// @param onlySelectedFeatures if true, only selected features are considered, else all the features
// @param p progress dialog (or 0 if no progress dialog is to be shown)
// @note: added in version 1.4*/
// bool symDifference( QgsVectorLayer* layerA, QgsVectorLayer* layerB,
// const QString& shapefileName, bool onlySelectedFeatures = false,
// QProgressDialog* p = 0 );
private:
void combineFieldLists( QgsFieldMap fieldListA, QgsFieldMap fieldListB );
void intersectFeature( QgsFeature& f, QgsVectorFileWriter* vfw, QgsVectorLayer* dp, QgsSpatialIndex* index );
};
#endif //QGSVECTORANALYZER