22.01.2016, 10:06, Kai Muehlbauer kirjoitti:
To quickly check if there is an improvement I used the python shapely
package. Converting the outer loop geometry to an shapely prepared
geometry and calling 'contains()' with the inner loop geometries to
find those fully contained geometries I get an improvement of nearly
40% to the former ogr only code.
As Even pointed out, there is only OGRPreparedGeometryIntersects()
(GEOSPreparedIntersects) exposed within ogrgeometry.cpp. In my use
case I would need GEOSPreparedContains to mimic the shapely behaviour.
I made a small test with the attached diff, which prepares the
geometries of the input layer.
The input layer had two polygons (boundaries of two Finnish
municipalities) and the method layer had 3091 points, all in either one
of the polygons.
The improvement is ~15%.
I think the prepared geometries could be used by default easily.
Prepared contains test would need addition of the respective GEOS method.
Ari
Next thing is that the ogr SWIG python bindings do not wrap those
available OGRPrepared-functions.
I'm no expert when it comes to c/c++-python interaction to create a
patch to have this functionality within python ogr/gdal. So hoping for
help in this matter, I'm asking for suggestions how to get the
GEOSPreparedGeometry stuff called from OGR.
Cheers,
Kai
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Index: ogr/ogrsf_frmts/generic/ogrlayer.cpp
===================================================================
--- ogr/ogrsf_frmts/generic/ogrlayer.cpp (revision 33064)
+++ ogr/ogrsf_frmts/generic/ogrlayer.cpp (working copy)
@@ -2033,6 +2033,7 @@
double progress_ticker = 0;
int bSkipFailures = CPLTestBool(CSLFetchNameValueDef(papszOptions, "SKIP_FAILURES", "NO"));
int bPromoteToMulti = CPLTestBool(CSLFetchNameValueDef(papszOptions, "PROMOTE_TO_MULTI", "NO"));
+ int bUsePreparedGeometries = CPLTestBool(CSLFetchNameValueDef(papszOptions, "USE_PREPARED_GEOMETRIES", "NO"));
// check for GEOS
if (!OGRGeometryFactory::haveGEOS()) {
@@ -2093,10 +2094,23 @@
continue;
}
+ //
+ OGRPreparedGeometry* x_prepared_geom = NULL;
+ if (bUsePreparedGeometries) x_prepared_geom = OGRCreatePreparedGeometry(x_geom);
+
pLayerMethod->ResetReading();
while (OGRFeature *y = pLayerMethod->GetNextFeature()) {
OGRGeometry *y_geom = y->GetGeometryRef();
if (!y_geom) {delete y; continue;}
+
+ if (bUsePreparedGeometries &&
+ !((x_prepared_geom && OGRPreparedGeometryIntersects(x_prepared_geom, y_geom)) ||
+ x_geom->Intersects(y_geom)))
+ {
+ delete y;
+ continue;
+ }
+
OGRGeometry* poIntersection = x_geom->Intersection(y_geom);
if( poIntersection == NULL || poIntersection->IsEmpty() ||
(x_geom->getDimension() == 2 &&
@@ -2119,6 +2133,7 @@
delete z;
if (ret != OGRERR_NONE) {
if (!bSkipFailures) {
+ OGRDestroyPreparedGeometry(x_prepared_geom);
delete x;
goto done;
} else {
@@ -2129,6 +2144,7 @@
}
}
+ OGRDestroyPreparedGeometry(x_prepared_geom);
delete x;
}
if (pfnProgress && !pfnProgress(1.0, "", pProgressArg)) {
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
