2

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.

enter image description here

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'

  • If you've got Esri, and it's working in Esri, why OGR2OGR? If you intention is to go open source why not select by location in QGIS and export selected records... – Michael Stimson Dec 9 '14 at 21:34
  • The Esri processing is taking way too long (weeks vs. ogr2ogr hours/days). I have 5 buffer sizes for 80,000+ grid points and I need to select only the exploded buffer features that can be selected by the grid point (buffer centroid). – GIS Danny Dec 10 '14 at 2:34
  • Is it one polygon (buffer) for each size per point? It seems to me that the buffers are created by the points then have been erased by some other data. OGR2OGR is somewhat limited in this respect however it is quicker, would you consider using Near in Esri on the points to find the closest polygon then selecting by attributes (FID = Near_FID)? The other option is moving the shapefile to spatial lite or PostGis which do support spatial operators. – Michael Stimson Dec 10 '14 at 2:41
  • That's an interesting approach. I have some stuff processing so I'll try and do some tests on that approach. I'm down this rabbit hole sort of b/c I know it works. Can MakeFeatureLayer_management and SelectLayerByLocation_management be improved using in_memory? – GIS Danny Dec 10 '14 at 16:52
  • This is a similar processing function - saying to use Spatial Join - gis.stackexchange.com/questions/12448/… – GIS Danny Dec 10 '14 at 17:08
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']
                        })
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.