[gdal-dev] Tests if the geometry is within the other geometry

[Even Rouault]

On samedi 23 juin 2018 04:21:38 CEST Ian wrote:
> I have a shape file consisting of thousands of features/polygons.
> I have to extract only those polygons which are within the reference
> polygon.
> 
> I tried as follows:
> 
> refpoly = feat.GetGeometryRef()
> 
> fi =  'data.shp'
> ids = ogr.Open(fi,0)
> inlay = ids.GetLayer(0)
> for feat in inlay:
>     ctype = feat.GetField('Name')
>     geomi = feat.GetGeometryRef()
>        if geomi.Within(refpoly):
>            write ods....
> 
> But it is running for weeks and not yet finished.
> How to do it faster?

geomi.Within(refpoly) is equivalent to refpoly.Contains(geomi)

If you use Shapely and prepared geometry, and that refpoly is a rather complex 
geometry, then making it a prepared geometry could speed up the contains test.

http://toblerity.org/shapely/manual.html#prepared-geometry-operations

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com