[gdal-dev] Zonal statistics with GDAL/StarSpan
Brent Fraser
bfraser at geoanalytic.com
Fri Dec 19 18:13:11 EST 2008
> Hi,
>
> I performed the OpenJump procedure. But, I still get:
>
> Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none
> --vector afadmn2n.shp --fid 3 --raster rfe_2006_04_pct.tif --stats
> results.txt avg --out-prefix bar --out-type table
> Number of features: 550
> starspan_csv: 1: Extracting from rfe_2006_04_pct.tif
> 0% terminate called after throwing an instance of
> 'geos::util::TopologyException'
> what(): TopologyException: side location conflict -8.25 23.1536
> Abort trap
> So:
>
> a) Is there some way to get a dump of FIDs with their corresponding
> country/province names?
You could use dbfdump. I think it is delivered with FWTools.
> b) What other steps might I take to detect/fix any side location
> conflict(s)
I've attached the source to a GDAL application I wrote (it simply accesses
the GDAL/Geos buffer function) to repair geometry of polygons. Use the -r
option (the -c option doesn't work). You'll need to link with a version
of GDAL that uses GEOS.
> c) The docs on the starspan wiki are out of sync with the current
> version (1.2.03). Are there updated docs for starspan? Should I fall
> back to a previous version?
>
I used version 1.2.01, but I don't recall why...
Brent
-------------- next part --------------
/* ***************** SOURCE HEADER: **********************
Name:
Purpose:
Invocation:
Arguments:
Data Files Accessed:
Program History:
v1.0 19 Jun 2007 BWF initial coding
Note:
Algorithm Description:
Future Enhancements:
*************** End of SOURCE HEADER ************** */
#include <assert.h>
#include "GDAL.h"
#include "cpl_string.h"
#include "ogrsf_frmts.h"
/************************************************************************/
/************************************************************************/
/************************************************************************/
static void Usage()
{
printf(
"ogrClean v1.0\n"
"Create clean planar polygon layer from overlapping polygons.\n"
"Usage: ogrClean -a <input> Assigned unique integer to FID attribute\n"
" ogrClean -r <input> Repair topology\n"
" ogrClean -c <input> Create non-overlapping topology\n"
" ogrClean -d <input> <output> Delete records with NULL geometry.\n"
);
exit( 1 );
}
/************************************************************************/
/************************************************************************/
void notice(const char *fmt, ...)
{
va_list ap;
fprintf( stdout, "NOTICE: ");
va_start (ap, fmt);
vfprintf( stdout, fmt, ap);
va_end(ap);
fprintf( stdout, "\n" );
}
void log_and_exit(const char *fmt, ...)
{
va_list ap;
fprintf( stdout, "ERROR: ");
va_start (ap, fmt);
vfprintf( stdout, fmt, ap);
va_end(ap);
fprintf( stdout, "\n" );
exit(1);
}
/************************************************************************/
/************************************************************************/
int repairPolygons(const char *pszSrcFilename, bool bOpAssignFID, bool bOpRepairTopo )
{
OGRDataSource *poDS;
OGRLayer *poLayer;
OGRFeature *poFeature;
OGRGeometry *poGeometry;
OGRwkbGeometryType poGeomType;
char *pszGeom=NULL;
OGRErr oErr;
long nLayers=0;
long nFeatures=0;
bool bFoundFID;
int iField;
long i=0, id=0;
long iCountFixed=0;
bool bGeomUpdated = false;
/*-----------------------------------------------------------------------*/
/* Open input: */
/*-----------------------------------------------------------------------*/
OGRRegisterAll();
poDS = OGRSFDriverRegistrar::Open( pszSrcFilename, TRUE ); // we will be updating the dataset
if( poDS == NULL )
{
return(1); //can't open dataset
}
nLayers = poDS->GetLayerCount();
if (nLayers == 0)
{
return(2); //no layers in dataset
}
// TBD: allow mult layers:
poLayer = poDS->GetLayer( 0 );
poLayer->ResetReading();
nFeatures = poLayer->GetFeatureCount();
if ( (poFeature = poLayer->GetNextFeature()) != NULL )
{
//--- Check Geometry type: ---//
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL )
{
poGeomType = poGeometry->getGeometryType(); // TBD: check if polgon
if ( poGeomType != wkbMultiPolygon && poGeomType != wkbPolygon )
{
printf("Data must be polygons.\n");
}
}else
{
return(4); //no geometry in feature
}
//--- Get Field of unique ID: ---//
bFoundFID = false;
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if ( EQUAL(poFieldDefn->GetNameRef(),"FID") )
{
bFoundFID = true;
break;
}
}
if (!bFoundFID)
{
printf("Data must have a unique ID field called FID.\n");
return(5);
}
}else
{
printf("No features in dataset.\n");
return(3); //no features in dataset
}
/*-----------------------------------------------------------------------*/
/* Process Features: */
/*-----------------------------------------------------------------------*/
poLayer->ResetReading();
i=0;
while( (poFeature = poLayer->GetNextFeature()) != NULL )
{
GDALTermProgress( ((double)i)/nFeatures, "", NULL);
if ( bOpAssignFID )
{
id = i;
}else
{
id = poFeature->GetFieldAsInteger(iField);
}
// printf("%d. FID %d\n", i, id);
if ( bOpRepairTopo && poFeature->GetGeometryRef() ) //not null
{
bGeomUpdated = false;
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef()->clone(); // make a copy of geometry, otherwise there's a problem re-assigning it later.
poGeomType = poGeometry->getGeometryType();
if ( poGeomType != wkbMultiPolygon && poGeomType != wkbPolygon ) // not a polygon!
{
printf("** Warning. FID %d has geometry type of %s.\n", id, poGeometry->getGeometryName() );
poGeometry->exportToWkt(&pszGeom);
printf("%s\n",pszGeom );
}
if ( !poGeometry->IsValid() ) // need to fix, so...
{
// printf("Warning. FID %d has invalid topology.\n", id);
OGRGeometry *poTempGeometry = poGeometry->Buffer(0.0); // Let Geos fix the problem.
if (poTempGeometry) // not sure how it could be null...
{
if ( poTempGeometry->IsValid() )
{
// printf(" Can correct by buffering.\n");
poGeometry = poTempGeometry->clone(); //replace with this new geom (need to destroy old?)
oErr = poFeature->SetGeometry(poGeometry);
bGeomUpdated = true;
iCountFixed++;
}else
{
printf("** Warning. FID %d still has invalid topology after buffering (so not updated).\n", id);
}
}else
{
printf("** Warning. FID %d has NULL geometry after buffering (so not updated).\n", id);
}
if (poTempGeometry != NULL)
{
OGRGeometryFactory::destroyGeometry( poTempGeometry);
poTempGeometry = NULL;
}
}else // valid
{
if ( poGeometry->IsEmpty() )
{
// printf("Warning. FID %d has empty geometry.\n", id);
}
}
if (poGeometry != NULL) // can delete since SetGeometry makes a copy
{
OGRGeometryFactory::destroyGeometry( poGeometry);
poGeometry = NULL;
}
}else
{
// printf("Warning. FID %d has NULL geometry.\n", id);
}
if ( bOpAssignFID )
{
poFeature->SetField( iField, id );
}
if ( bGeomUpdated || bOpAssignFID )
{
oErr = poLayer->SetFeature(poFeature);
// poDS->SyncToDisk();
}
OGRFeature::DestroyFeature( poFeature );
i++; //first record is 0
}
GDALTermProgress( ((double)nFeatures)/nFeatures, "", NULL);
if ( bOpRepairTopo )
printf("\n-- Info: Fixed %d polygons.\n", iCountFixed );
OGRDataSource::DestroyDataSource( poDS );
return(0);
}
/************************************************************************/
/************************************************************************/
int cleanPolygons(const char *pszSrcFilename, const char *pszAdjFilename ) //, const char *pszDstFilename, const char *pzsDstFormat)
{
OGRDataSource *poDS;
OGRLayer *poLayer;
OGRFeature *poFeature;
OGRDataSource *poDSadjacent;
OGRLayer *poLayerAdjacent;
OGRFeature *poFeatureAdjacent;
long nLayers;
OGRwkbGeometryType poGeomType;
OGRGeometry *poGeometry;
char *pszGeom=NULL;
OGRErr oErr;
long i=0;
long nFeatures=0;
unsigned int nCoords=0;
// double x, y;
int iField;
int id, idAdjacent;
bool bFoundFID;
bool bDone=false;
bool bGeomUpdated = false;
/*-----------------------------------------------------------------------*/
/* Open input: */
/*-----------------------------------------------------------------------*/
OGRRegisterAll();
poDS = OGRSFDriverRegistrar::Open( pszSrcFilename, TRUE ); // we will be updating the dataset
if( poDS == NULL )
{
return(1); //can't open dataset
}
nLayers = poDS->GetLayerCount();
if (nLayers == 0)
{
return(2); //no layers in dataset
}
// TBD: allow mult layers:
poLayer = poDS->GetLayer( 0 );
poLayer->ResetReading();
nFeatures = poLayer->GetFeatureCount();
if ( (poFeature = poLayer->GetNextFeature()) != NULL )
{
//--- Check Geometry type: ---//
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL )
{
poGeomType = poGeometry->getGeometryType(); // TBD: check if polgon
if ( poGeomType != wkbMultiPolygon && poGeomType != wkbPolygon )
{
printf("Data must be polygons.\n");
}
}else
{
return(4); //no geometry in feature
}
//--- Get Field of unique ID: ---//
bFoundFID = false;
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if ( EQUAL(poFieldDefn->GetNameRef(),"FID") )
{
bFoundFID = true;
break;
}
}
if (!bFoundFID)
{
printf("Data must have a unique ID field called FID.\n");
return(5);
}
}else
{
printf("No features in dataset.\n");
return(3); //no features in dataset
}
poLayer->ResetReading();
// poDSadjacent = OGRSFDriverRegistrar::Open( pszSrcFilename, FALSE ); // read only
// if( poDSadjacent == NULL ) return(1); //can't open dataset
// poLayerAdjacent = poDSadjacent->GetLayer( 0 );
/*-----------------------------------------------------------------------*/
/* Process Features: */
/*-----------------------------------------------------------------------*/
i=0;
while( (poFeature = poLayer->GetNextFeature()) != NULL )
{
i++;
GDALTermProgress( ((double)i)/nFeatures, "", NULL);
OGRGeometry *poGeometry;
id = poFeature->GetFieldAsInteger(iField);
// printf("%d. FID %d\n", i, id);
bGeomUpdated = false;
if ( poFeature->GetGeometryRef() ) //not null
{
poGeometry = poFeature->GetGeometryRef()->clone(); // make a copy of geometry, otherwise there's a problem re-assigning it later.
if (!poGeometry->IsValid() )
{
printf("Warning. FID %d has invalid geometry, buffering to fix.\n", id);
poGeometry = poGeometry->Buffer(0.0);
}
//---- Set Spatial Filter to only this Poly-of-interest ---//
if ( pszAdjFilename )
{
poDSadjacent = OGRSFDriverRegistrar::Open( pszAdjFilename, FALSE ); // read only
}else
{
poDSadjacent = OGRSFDriverRegistrar::Open( pszSrcFilename, FALSE ); // read only
}
if( poDSadjacent == NULL ) return(1); //can't open dataset
poLayerAdjacent = poDSadjacent->GetLayer( 0 );
poLayerAdjacent->SetSpatialFilter(poGeometry);
poLayerAdjacent->ResetReading();
//--- Get Adjacent Polys: ---//
bDone = false;
while( ((poFeatureAdjacent = poLayerAdjacent->GetNextFeature()) != NULL) && ! bDone ) // more polys nearby and some poly geom left
{
OGRGeometry *poGeometryAdjacent = poFeatureAdjacent->GetGeometryRef(); // no need to clone since need read only
idAdjacent = poFeatureAdjacent->GetFieldAsInteger(iField);
if ( poGeometryAdjacent ) //might have null geom? YES (due to this processing)
{
if (!poGeometryAdjacent->IsValid() )
{
poGeometryAdjacent = poGeometryAdjacent->Buffer(0.0);
}
if ( ( id!=idAdjacent || pszAdjFilename ) && (
!poGeometryAdjacent->IsEmpty() &&
poGeometryAdjacent->IsValid() &&
poGeometryAdjacent->Intersects(poGeometry)) ) //has some geom and not the same poly
{
OGRGeometry *poTempGeometry = poGeometry->Difference(poGeometryAdjacent); // subtract the adj poly
if (poTempGeometry) // not sure how it could be null...
{
poGeomType = poTempGeometry->getGeometryType(); // make sure a poly is the result (will be linestring for shared boundary)
if ( poGeomType == wkbMultiPolygon || poGeomType == wkbPolygon )
{
if (!poTempGeometry->IsValid() )
{
poTempGeometry = poTempGeometry->Buffer(0.0);
if (!poTempGeometry->IsValid() )
{
printf("Warning. FID %d - FID %d resulted in invalid geometry, even after buffering.\n", id, idAdjacent);
}
}else
{
poGeometry = poTempGeometry->clone(); //replace with this new geom (need to destroy old?)
//poGeometry->exportToWkt(&pszGeom);
//printf("%d - %d: %s\n", id, idAdjacent, pszGeom);
bGeomUpdated = true;
}
}else // not a polygon
{
if ( !poTempGeometry->IsEmpty())
{
printf("Warning. FID %d now has geometry type of %s.\n", id, poTempGeometry->getGeometryName() );
poTempGeometry->exportToWkt(&pszGeom);
printf("%s\n",pszGeom );
}else // empty "other" (GEOMETERYCOLLECTION) geometry
{
// printf("Warning. FID %d now has empty geometry.\n", id );
poGeometry = poTempGeometry->clone(); //replace with this new geom (need to destroy old?)
bGeomUpdated = true;
bDone = true;
}
}
}else // null polygon geom
{
poGeometry = poTempGeometry;
printf("Warning. FID %d now has NULL geometry.\n", id );
bGeomUpdated = true;
bDone = true;
}
}else
{
if (poGeometryAdjacent->IsEmpty())
{
printf(" Warning. Adjacent FID %d has empty geometry.\n", idAdjacent);
}
if (!poGeometryAdjacent->IsValid()) // attempted to buffer above, so don't bother again
{
printf(" Warning. Adjacent FID %d has invalid geometry.\n", idAdjacent);
}
}
OGRGeometryFactory::destroyGeometry( poGeometryAdjacent ); // done with this nearby poly
poGeometryAdjacent = NULL;
}else // null geom
{
//printf(" Warning. Adjacent FID %d has NULL geometry.\n", idAdjacent);
}
//OGRFeature::DestroyFeature( poFeatureAdjacent );
} //no more nearby polys
OGRDataSource::DestroyDataSource( poDSadjacent );
if ( bGeomUpdated )
{
if ( poGeometry )
{
//--- Set the new clipped geometry: ---//
if (poGeometry->getGeometryType() == wkbMultiPolygon || poGeometry->getGeometryType() == wkbPolygon )
{
if (!poGeometry->IsValid() )
{
printf("Warning. FID %d has invalid geometry, buffering to fix.\n", id);
poGeometry = poGeometry->Buffer(0.0);
}
oErr = poFeature->SetGeometry(poGeometry);
}else
{
oErr = poFeature->SetGeometry(NULL); //tbd: delete feature?
}
if (poGeometry->IsEmpty())
{
oErr = poFeature->SetGeometry(NULL); //tbd: delete feature?
}
} else
{
oErr = poFeature->SetGeometry(NULL); //tbd: delete feature?
}
oErr = poLayer->SetFeature(poFeature);
poDS->SyncToDisk();
}
if (poGeometry != NULL)
{
OGRGeometryFactory::destroyGeometry( poGeometry);
poGeometry = NULL;
}
}else // null geom
{
// printf("Warning. FID %d has NULL geometry.\n", id );
}
OGRFeature::DestroyFeature( poFeature );
} // end of polys
/*-----------------------------------------------------------------------*/
/* close Datasources: */
/*-----------------------------------------------------------------------*/
// OGRDataSource::DestroyDataSource( poDSadjacent );
OGRDataSource::DestroyDataSource( poDS );
return(0);
}
/************************************************************************/
/************************************************************************/
int deleteNullPolygons(const char *pszSrcFilename, const char *pszDstFilename, const char *pzsDstFormat)
{
OGRDataSource *poDS;
OGRLayer *poLayer;
OGRFeature *poFeature;
OGRFeatureDefn *poFDefn;
OGRwkbGeometryType poGeomType;
OGRGeometry *poGeometry;
OGRDataSource *poDSout;
OGRLayer *poLayerOut;
OGRFeature *poFeatureOut;
long nLayers, nFeatures;
int i;
int iFiltered=0;
bool bFoundFID=false;
bool bAssignFID=false;
int iField=0;
/*-----------------------------------------------------------------------*/
/* Open input: */
/*-----------------------------------------------------------------------*/
OGRRegisterAll();
poDS = OGRSFDriverRegistrar::Open( pszSrcFilename, FALSE ); // we will NOT be updating the dataset
if( poDS == NULL )
{
printf("Failed to open dataset: %s\n",pszSrcFilename);
return(1); //can't open dataset
}
nLayers = poDS->GetLayerCount();
if (nLayers == 0)
{
return(2); //no layers in dataset
}
// TBD: allow mult layers:
poLayer = poDS->GetLayer( 0 );
poLayer->ResetReading();
nFeatures = poLayer->GetFeatureCount();
if ( (poFeature = poLayer->GetNextFeature()) != NULL )
{
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL )
{
poGeomType = poGeometry->getGeometryType();
}else
{
return(4); //no geometry in feature
}
}else
{
return(3); //no features in dataset
}
//--- Get Field of unique ID: ---//
bFoundFID = false;
poFDefn = poLayer->GetLayerDefn();
for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if ( EQUAL(poFieldDefn->GetNameRef(),"FID") )
{
bFoundFID = true;
break;
}
}
/*-----------------------------------------------------------------------*/
/* Create output: */
/*-----------------------------------------------------------------------*/
OGRSFDriver *poDriver;
poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName( pzsDstFormat );
if( poDriver == NULL )
{
// printf( "%s driver not available.\n", pszDriverName );
return(101);
}
poDSout = poDriver->CreateDataSource( pszDstFilename, NULL );
if( poDSout == NULL )
{
// printf( "Creation of output file failed.\n" );
return(102);
}
poLayerOut = poDSout->CreateLayer( "new_layer", NULL, poGeomType, NULL );
if( poLayerOut == NULL )
{
// printf( "Layer creation failed.\n" );
return(103);
}
poFDefn = poLayer->GetLayerDefn();
for(iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if( poLayerOut->CreateField( poFieldDefn ) != OGRERR_NONE )
{
// printf( "Creation of field failed.\n" );
return(104);
}
}
if (!bFoundFID)
{
OGRFieldDefn oField( "FID", OFTInteger );
bAssignFID = true;
if( poLayerOut->CreateField( &oField ) != OGRERR_NONE )
{
printf( "Creating FID field failed.\n" );
return(5);
}
}
bFoundFID = false;
poFDefn = poLayerOut->GetLayerDefn();
for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if ( EQUAL(poFieldDefn->GetNameRef(),"FID") )
{
bFoundFID = true;
break;
}
}
/*-----------------------------------------------------------------------*/
/* Process Features: */
/*-----------------------------------------------------------------------*/
poLayer->ResetReading();
poFDefn = poLayer->GetLayerDefn();
i = 0;
while( (poFeature = poLayer->GetNextFeature()) != NULL )
{
i++;
GDALTermProgress( ((double)i)/nFeatures, "", NULL);
// OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL && !poGeometry->IsEmpty() )
{
poFeatureOut = OGRFeature::CreateFeature( poLayerOut->GetLayerDefn() );
poFeatureOut->SetFrom(poFeature);
if (bAssignFID)
{
poFeatureOut->SetField( iField, i );
}
if( poLayerOut->CreateFeature( poFeatureOut ) != OGRERR_NONE )
{
printf( "Failed to create feature.\n" );
return( 200 );
}
}else
{
iFiltered++;
}
OGRFeature::DestroyFeature( poFeature );
}
GDALTermProgress( ((double)nFeatures)/nFeatures, "", NULL);
printf("\n-- Info: Filtered %d NULL or Empty polygons from a total of %d polygons.\n", iFiltered,i );
OGRDataSource::DestroyDataSource( poDS );
OGRDataSource::DestroyDataSource( poDSout );
return(0);
}
/************************************************************************/
/* Main() */
/************************************************************************/
int main( int argc, char *argv[] )
{
bool bOpAssignFID = false;
bool bOpRepairTopo = false;
bool bOpCleanPolys = false;
bool bOpDelNullPolys = false;
int i;
int rc=0;
const char *pszSrcFilename = NULL;
const char *pszAdjFilename = NULL;
const char *pszDstFilename = NULL;
char *pzsDstFormat = NULL;
char *pzsDstFormatDefault = "ESRI Shapefile";
/* -------------------------------------------------------------------- */
/* Parse arguments. */
/* -------------------------------------------------------------------- */
argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
if( argc < 1 )
return( -argc );
if( argc == 1 )
Usage();
for( i = 1; i < argc; i++ )
{
if( EQUAL(argv[i],"-a") )
bOpAssignFID = true;
else if( EQUAL(argv[i],"-r") )
bOpRepairTopo = true;
else if( EQUAL(argv[i],"-c") )
bOpCleanPolys = true;
else if( EQUAL(argv[i],"-d") )
bOpDelNullPolys = true;
else if( EQUAL(argv[i],"-o") && i < argc-1 )
pszAdjFilename = argv[++i];
else if( EQUAL(argv[i],"-of") && i < argc-1 )
pzsDstFormat = argv[++i];
else if( pszSrcFilename == NULL )
{
pszSrcFilename = argv[i];
}
else if( pszDstFilename == NULL )
{
pszDstFilename = argv[i];
}
else
Usage();
}
if( pszSrcFilename == NULL )
{
CPLError( CE_Failure, CPLE_AppDefined, "An input file is required." );
return(1);
}
if( pszDstFilename == NULL && bOpDelNullPolys )
{
CPLError( CE_Failure, CPLE_AppDefined, "An output file is required when deleting null geometry records." );
return(1);
}
if ( pzsDstFormat == NULL )
{
pzsDstFormat = pzsDstFormatDefault;
}
/* -------------------------------------------------------------------- */
/* Do requested Operations: */
/* -------------------------------------------------------------------- */
if ( bOpAssignFID || bOpRepairTopo )
rc = repairPolygons(pszSrcFilename, bOpAssignFID, bOpRepairTopo );
if ( bOpCleanPolys )
rc = cleanPolygons(pszSrcFilename, pszAdjFilename);
if ( bOpDelNullPolys )
rc = deleteNullPolygons(pszSrcFilename, pszDstFilename, pzsDstFormat);
/* -------------------------------------------------------------------- */
/* Exit: */
/* -------------------------------------------------------------------- */
return(0);
}
More information about the gdal-dev
mailing list