OGR/Python- Check whether or not points are in a polygon

Question

I have two shapefiles, one points file and one polygon. They are in the same projection (Lat/Long). I would like to run a quick command-line command or a Python script with the OGR module that creates a new points shapefile of only the points that are within that polygon. Is there a good way to do this?

Answer 1

You can do that with ogr2ogr, both on Python and on the command line. You first need to create a virtual OGR datset with your two shapefiles, points.shp and grid.shp, which you'll describe as a file, let's call it data.vrt. The contents of data.vrt should just be (change filenames etc as required):

<OGRVRTDataSource>
    <OGRVRTLayer name="grid">
        <SrcDataSource>grid.shp</SrcDataSource>
        <SrcLayer>grid</SrcLayer>
    </OGRVRTLayer>
    <OGRVRTLayer name="points">
        <SrcDataSource>points.shp</SrcDataSource>
        <SrcLayer>points</SrcLayer>
    </OGRVRTLayer>
</OGRVRTDataSource>

In the example above, the name in the OGRVRTLayer is the name of the output virtual layer (you can set it to whatever you want), the SrcDataSource points to the shapefile (you can add the full path to the shapefile here). SrcLayer is just the shapefile name minus the .shp extension.

You can check your file is sensible using

ogrinfo -al data.vrt

This will spit out tons of information to the shell. Then, happy with the VRT file being OK, you can write some simple SQL to look for all points that intercept with that polygon:

ogr2ogr -f GeoJSON -dialect sqlite \
 -sql "select points.* from points, grid 
where ST_Intersects(grid.geometry, points.geometry)" sel_points.geojson data.vrt 

I have done this with GeoJSON files (slightly simpler VRT file) an you can see it here.