Working with Shapefiles containing over 200K, and stepping through code that checked for point containment within polygons, I found a couple optimizations.
The primary optimization is in the ReadGeosGeometries method. This original implementation used to iterate the shapefile, generating GEOS geometries for every shape, even though only a subset of them were required.
An earlier optimization created a list of ID's, being those shapes returned from the QTree operation of the reference shape (i.e. those within proximity of the Extents of the reference shape). It then generated GEOS geometries only for those shapes within proximity of the reference shape.
This optimization changes the order of iteration. Instead of iterating the entire shapefile (which may have hundreds of thousands of shapes), looking for the ID's returned by the QTree operation, it now iterates only the subset of ID's (which is often just a handful), only generating the GEOS geometries that are necessary for the geospatial operation.
Finally, if the QTree operation returns zero shapes, the code can simply bypass the remainder of the function, rather than enter multiple loops, of which no iterations will take place.