Finding an ogr2ogr equivalent of SelectLayerByLocation

Question

So I'm trying to take individual grid point buffers and only select the regions that the points intersect with (via SelectByLocation). Generally I run the following code to:

• Explode multipart features (not inlcuded in code below)
• create a Feature Layer
• Make the Selection and
• then turn the selected features in the Layer into a new Feature Class

Below is the code for ArcGIS's arcpy:

for a in geogList:
    for b in radbufListFn:
        for c in idList:
            layerName = a+b+"_erase_exp_"+c+"_Layer"
            arcpy.MakeFeatureLayer_management(wp+"geogs/split_geogs_explode.gdb/"+a+b+"_erase_exp_"+c,layerName)
            arcpy.SelectLayerByLocation_management(layerName,"INTERSECT",wp+"geogs/split_points.gdb/"+a+"_"+c,"#","NEW_SELECTION")
            arcpy.FeatureClassToFeatureClass_conversion(layerName,wp+"geogs/split_geogs_explode_select.gdb",a+b+'_'+c)

I've seen that there's a ST_Intersection and ST_Intersect but I'm not really sure these do what I need ogr2ogr. I need my output to be just the A part of the polygon.

Is there an ogr2ogr equivalent of SelectLayerByLocation?

edit: I found a blogpost here http://www.gisdoctor.com/site/2012/07/29/statement-spatial-join/

I'm not totally sure if this is the code I want, but I modified mine to this for testing purposes:

ogr2ogr -f "ESRI Shapefile" -sql "select PP.*, SP.* from point PP inner join buffer_clip_exp SP with(Index(geom_sidx)) on PP.geom.STIntersects(SP.geom) = 1" test_out.shp buffer_clip_exp.shp point.shp

And I get this error:

layer names ignored in combination with -sql. ERROR 1: SQL Expression Parsing Error: syntax error, unexpected reserved keyword, expecting end of string. Occured around : select PP., SP. from point PP inner join buffer_clip_exp SP with(Index

a caret (^) is below the 'i' in 'inner'

Answer 1

I found some code here Is there a python option to "join attributes by location"? that uses Shapely and Fiona that outputted the points and then rearranged the code to output the polygon instead.

This seems to be working on one of the 5 grid point buffer geographies and running through the ~80,000 grid points.

import fiona
from shapely.geometry import shape
from copy import deepcopy

with fiona.open("buffer_clip_exp.shp", "r") as n: 

    with fiona.open("point.shp", "r") as s:

        # create a schema for the attributes
        outSchema =  deepcopy(n.schema)
        outSchema['properties'].update(s.schema['properties'])

        with fiona.open ("TEST7.shp", "w", s.driver, outSchema, s.crs) as output:

            for school in s: 
                for neighborhood in n:
                    # check if point is in polygon and set attribute
                    if shape(school['geometry']).within(shape(neighborhood['geometry'])):  
                        #school['properties'] = neighborhood['properties'] 
                    # write out
                        output.write({                                 
                            'properties': neighborhood['properties'], 
                            'geometry': neighborhood['geometry']
                        })