Skip to content
Permalink
Browse files
add length() and area() methods to QgsGeometry
and use it in labeling and zonal statistics


git-svn-id: http://svn.osgeo.org/qgis/trunk@13588 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
jef committed May 29, 2010
1 parent 3b280fa commit c410f69bf924e6d58a3491be0a97b056ea268989
@@ -111,6 +111,16 @@ class QgsGeometry
// TODO: unsupported class... would be possible to use PyGEOS?
//void fromGeos(geos::Geometry* geos);

/** get area using GEOS
@note added in 1.5
*/
double area();

/** get length using GEOS
@note added in 1.5
*/
double length();

double distance(QgsGeometry& geom);

/**
@@ -23,13 +23,19 @@
#include "cpl_string.h"
#include <QProgressDialog>

QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand ) :
mRasterFilePath( rasterFile ), mRasterBand( rasterBand ), mPolygonLayer( polygonLayer ), mAttributePrefix( attributePrefix ), mInputNodataValue( -1 )
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand )
: mRasterFilePath( rasterFile )
, mRasterBand( rasterBand )
, mPolygonLayer( polygonLayer )
, mAttributePrefix( attributePrefix )
, mInputNodataValue( -1 )
{

}

QgsZonalStatistics::QgsZonalStatistics(): mRasterBand( 0 ), mPolygonLayer( 0 )
QgsZonalStatistics::QgsZonalStatistics()
: mRasterBand( 0 )
, mPolygonLayer( 0 )
{

}
@@ -162,13 +168,13 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
continue;
}

statisticsFromMiddlePointTest_improved( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, \
statisticsFromMiddlePointTest_improved( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
rasterBBox, sum, count );

if ( count <= 1 )
{
//the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, \
statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
rasterBBox, sum, count );
}

@@ -227,7 +233,7 @@ int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const Q
return 0;
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX, \
void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
{
double cellCenterX, cellCenterY;
@@ -260,20 +266,18 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry*
CPLFree( scanLine );
}

void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX, \
void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
{
sum = 0;
count = 0;
double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
QgsGeometry* pixelRectGeometry = 0;
QgsGeometry* intersectGeometry = 0;

double hCellSizeX = cellSizeX / 2.0;
double hCellSizeY = cellSizeY / 2.0;
double pixelArea = cellSizeX * cellSizeY;
double intersectionArea = 0;
double weight = 0;

for ( int row = 0; row < nCellsY; ++row )
@@ -286,10 +290,11 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeome
if ( pixelRectGeometry )
{
//intersection
intersectGeometry = pixelRectGeometry->intersection( poly );
QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly );
if ( intersectGeometry )
{
if ( GEOSArea( intersectGeometry->asGeos(), &intersectionArea ) )
double intersectionArea = intersectGeometry->area();
if ( intersectionArea >= 0.0 )
{
weight = intersectionArea / pixelArea;
count += weight;
@@ -305,7 +310,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeome
CPLFree( pixelData );
}

void QgsZonalStatistics::statisticsFromMiddlePointTest_improved( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, \
void QgsZonalStatistics::statisticsFromMiddlePointTest_improved( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
{
double cellCenterX, cellCenterY;
@@ -5859,6 +5859,43 @@ QgsMultiPolygon QgsGeometry::asMultiPolygon()
return polygons;
}

double QgsGeometry::area()
{
if ( !mGeos )
{
exportWkbToGeos();
}

double area;

try
{
if ( GEOSArea( mGeos, &area ) == 0 )
return -1.0;
}
CATCH_GEOS( -1.0 )

return area;
}

double QgsGeometry::length()
{
if ( !mGeos )
{
exportWkbToGeos();
}

double length;

try
{
if ( GEOSLength( mGeos, &length ) == 0 )
return -1.0;
}
CATCH_GEOS( -1.0 )

return length;
}
double QgsGeometry::distance( QgsGeometry& geom )
{
if ( !mGeos )
@@ -146,6 +146,16 @@ class CORE_EXPORT QgsGeometry
*/
bool isGeosEmpty();

/** get area of geometry using GEOS
@note added in 1.5
*/
double area();

/** get length of geometry using GEOS
@note added in 1.5
*/
double length();

double distance( QgsGeometry& geom );

/**
@@ -229,25 +229,19 @@ bool LayerSettings::checkMinimumSizeMM( const QgsRenderContext& ct, QgsGeometry*
return true;
}

GEOSGeometry* geosGeom = geom->asGeos();
if ( !geosGeom )
{
return true;
}

double mapUnitsPerMM = ct.mapToPixel().mapUnitsPerPixel() * ct.scaleFactor();
if ( featureType == QGis::Line )
{
double length;
if ( GEOSLength( geosGeom, &length ) )
double length = geom->length();
if ( length >= 0.0 )
{
return ( length >= ( minSize * mapUnitsPerMM ) );
}
}
else if ( featureType == QGis::Polygon )
{
double area;
if ( GEOSArea( geosGeom, &area ) )
double area = geom->area();
if ( area >= 0.0 )
{
return ( sqrt( area ) >= ( minSize * mapUnitsPerMM ) );
}

0 comments on commit c410f69

Please sign in to comment.