1

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?

1 Answer 1

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.

3
  • Jose, this looks like a great option but I think that I need more detail. How do I create the virtual OGR dataset? I see that you posted some script that I think is to be the virtual OGR dataset. I created a file called data.vrt in which I replaced the names of the .shp files with my own file names. Upon running the 2nd part of your solution in the command prompt, I get ERROR 3: Cannot open file 'sel_points.geojson' FAILURE: Unable to open datasource sel_points.geojson' with the following drivers.~long list of drivers~'
    – Casivio
    May 30, 2020 at 7:02
  • Clarified VRT file creation above
    – Jose
    May 31, 2020 at 15:44
  • What worked for me was ogr2ogr -f GeoJSON -dialect sqlite selectedPoints.geojson -sql "select points.* from points, grid where ST_Intersects(grid.geometry, points.geometry)" data.vrt I then subbed in ESRI Shapefile as the driver to get the desired .shp output. ogr2ogr -f "ESRI Shapefile" -dialect sqlite selectedPoints.shp -sql "select points.* from points, grid where ST_Intersects(grid.geometry, points.geometry)" data.vrt
    – Casivio
    Jul 1, 2020 at 23:34

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.