[gdal-dev] Intersection of two Layers

[Kai Muehlbauer]

Am 20.01.2016 um 07:35 schrieb Kai Muehlbauer:
> Am 19.01.2016 um 13:09 schrieb Even Rouault:
>> Le mardi 19 janvier 2016 08:40:30, Ari Jolma a écrit :
>>> 18.01.2016, 11:56, Kai Muehlbauer kirjoitti:
>>>> Hi,
>>>>
>>>> at the moment I calculate the intersection of two OGR Layers (A, B) by
>>>> iterating over features/geometries (GA) of layer A and then calling
>>>> GA.Intersection(GB) for every Feature/Geometry (GB) of B. I set a
>>>> spatial filter for layer B with every GA. The intersections
>>>> features/geometries are written to a result layer (C).
>>>>
>>>> In my use case the features/geometries of layer B are much smaller
>>>> than the features of layer A. Means that quite a big number of B's
>>>> geometries are fully contained in A's geometries.
>>>>
>>>> To speed things up, I check for containment GA.contains(GB). If GB is
>>>> contained in GA, I just copy GB to C, otherwise I calculate
>>>> GA.Intersection(GB) and write this to C.
>>>>
>>>> I found that I also could use the layer intersection function
>>>> (OGRLayer::Intersection). Unfortunately the intersection is calculated
>>>> for every feature/geometry.
>>>>
>>>> Would it be possible to extend the layer intersection function to
>>>> allow for containment check (eg. as an option) and just copy the
>>>> contained geometries without calculation of expensive intersection?
>>>> Suggestions, opinions?
>>>
>>> There's the Options argument, which could be used. Something like
>>> "PRETEST=CONTAINMENT" could indicate a wish to test containment before
>>> intersection.
>>
>> The use of GEOS prepared geometry operations could also potentially
>> help. OGR
>> has already mapped GEOSPreparedIntersects as
>> OGRPreparedGeometryIntersects and
>> that could probably be unconditionnaly used, with the geometry of the
>> outer
>> loop being the prepared geometry with OGRCreatePreparedGeometry()
>> And perhaps GEOSPreparedContains could be used for containement (if
>> the above
>> suggestion desn't speed up things enough).
>
> Thank you Ari and Even for the comments and suggestions. I'll try to
> test the "prepared" stuff first, hoping that those functions are also
> available within the python bindings (which I use).
>
> I'll give an update here, when I have results on this.

OK, here we go.

JFTR, I'am using the SWIG python bindings for ogr.

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.

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