22.01.2016, 13:49, Kai Muehlbauer kirjoitti:
This is very good news! If I understand correctly, I could use your
patch and self-compile gdal to have access to the improved
layer-intersection function via the python bindings. I'll give this a
try.
Attached is probably a bit better diff.
As I wrote, similar thing should be done for other layer methods.
Yes, you can use this with the trunk version.
I also created a ticket here: https://trac.osgeo.org/gdal/ticket/6323
to keep track of this enhancement.
ok
Ari
Cheers,
Kai
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Index: ogr/ogr_geometry.h
===================================================================
--- ogr/ogr_geometry.h (revision 33064)
+++ ogr/ogr_geometry.h (working copy)
@@ -1343,5 +1343,7 @@
void OGRDestroyPreparedGeometry( OGRPreparedGeometry* poPreparedGeom );
int OGRPreparedGeometryIntersects( const OGRPreparedGeometry* poPreparedGeom,
const OGRGeometry* poOtherGeom );
+int OGRPreparedGeometryContains( const OGRPreparedGeometry* poPreparedGeom,
+ const OGRGeometry* poOtherGeom );
#endif /* ndef OGR_GEOMETRY_H_INCLUDED */
Index: ogr/ogrgeometry.cpp
===================================================================
--- ogr/ogrgeometry.cpp (revision 33064)
+++ ogr/ogrgeometry.cpp (working copy)
@@ -4765,6 +4765,28 @@
#endif
}
+int OGRPreparedGeometryContains( UNUSED_IF_NO_GEOS const OGRPreparedGeometry* poPreparedGeom,
+ UNUSED_IF_NO_GEOS const OGRGeometry* poOtherGeom )
+{
+#ifdef HAVE_GEOS_PREPARED_GEOMETRY
+ if( poPreparedGeom == NULL || poOtherGeom == NULL )
+ return FALSE;
+
+ GEOSGeom hGEOSOtherGeom = poOtherGeom->exportToGEOS(poPreparedGeom->hGEOSCtxt);
+ if( hGEOSOtherGeom == NULL )
+ return FALSE;
+
+ int bRet = GEOSPreparedContains_r(poPreparedGeom->hGEOSCtxt,
+ poPreparedGeom->poPreparedGEOSGeom,
+ hGEOSOtherGeom);
+ GEOSGeom_destroy_r( poPreparedGeom->hGEOSCtxt, hGEOSOtherGeom );
+
+ return bRet;
+#else
+ return FALSE;
+#endif
+}
+
/************************************************************************/
/* OGRGeometryFromEWKB() */
/************************************************************************/
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,8 @@
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"));
+ int bPretestContains = CPLTestBool(CSLFetchNameValueDef(papszOptions, "PRETEST_CONTAINS", "NO"));
// check for GEOS
if (!OGRGeometryFactory::haveGEOS()) {
@@ -2093,42 +2095,61 @@
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;}
- OGRGeometry* poIntersection = x_geom->Intersection(y_geom);
- if( poIntersection == NULL || poIntersection->IsEmpty() ||
- (x_geom->getDimension() == 2 &&
- y_geom->getDimension() == 2 &&
- poIntersection->getDimension() < 2) )
- {
- delete poIntersection;
- delete y;
+ OGRGeometry *z_geom = NULL;
+
+ if (bUsePreparedGeometries && x_prepared_geom) {
+ if (bPretestContains && OGRPreparedGeometryContains(x_prepared_geom, y_geom))
+ {
+ z_geom = y_geom->clone();
+ }
+ else if (!(OGRPreparedGeometryIntersects(x_prepared_geom, y_geom)))
+ {
+ delete y;
+ continue;
+ }
}
- else
- {
- OGRFeature *z = new OGRFeature(poDefnResult);
- z->SetFieldsFrom(x, mapInput);
- z->SetFieldsFrom(y, mapMethod);
- if( bPromoteToMulti )
- poIntersection = promote_to_multi(poIntersection);
- z->SetGeometryDirectly(poIntersection);
- delete y;
- ret = pLayerResult->CreateFeature(z);
- delete z;
- if (ret != OGRERR_NONE) {
- if (!bSkipFailures) {
- delete x;
- goto done;
- } else {
- CPLErrorReset();
- ret = OGRERR_NONE;
- }
+ if (!z_geom) {
+ z_geom = x_geom->Intersection(y_geom);
+ if (z_geom == NULL || z_geom->IsEmpty() ||
+ (x_geom->getDimension() == 2 &&
+ y_geom->getDimension() == 2 &&
+ z_geom->getDimension() < 2))
+ {
+ delete z_geom;
+ delete y;
+ continue;
}
}
+ OGRFeature *z = new OGRFeature(poDefnResult);
+ z->SetFieldsFrom(x, mapInput);
+ z->SetFieldsFrom(y, mapMethod);
+ if( bPromoteToMulti )
+ z_geom = promote_to_multi(z_geom);
+ z->SetGeometryDirectly(z_geom);
+ delete y;
+ ret = pLayerResult->CreateFeature(z);
+ delete z;
+ if (ret != OGRERR_NONE) {
+ if (!bSkipFailures) {
+ OGRDestroyPreparedGeometry(x_prepared_geom);
+ delete x;
+ goto done;
+ } else {
+ CPLErrorReset();
+ ret = OGRERR_NONE;
+ }
+ }
}
+ 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
