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?
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,
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
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.