OGR – How to Count Number of Points Within Multipolygon Shapefiles

countgdalogrpythonspatial-join

I have two shapefiles: 1) consists of points (hundreds of entries) and 2) multi polygons (tens of thousands of entries).

I simply need to count the total number of points contained within the polygon set (presumably using ogrinfo) and write out a shapefile containing these points.

Must use OGR/GDAL to do this. I feel that there is an ogr2ogr sql statement that allows me to select the points without performing a full on clip, which is expensive in terms of computation time (clipping the points by the polygon extent takes hours on hours and I cannot tell if the operation works). Or perhaps the best way would be a spatial join?

There is a similar question here:
https://stackoverflow.com/questions/29769991/python-ogr-select-and-write-points-which-intersect-with-polygon-layer/33399208#33399208

Best Answer

If pol.shp contains your polygons and pnt.shp - points:

Create input.vrt:

<OGRVRTDataSource>
<OGRVRTLayer name="pol">
    <SrcDataSource>pol.shp</SrcDataSource>
    <SrcLayer>pol</SrcLayer>
</OGRVRTLayer>
<OGRVRTLayer name="pnt">
    <SrcDataSource>pnt.shp</SrcDataSource>
    <SrcLayer>pnt</SrcLayer>
</OGRVRTLayer>
</OGRVRTDataSource>

Run the query:

ogrinfo -sql "SELECT COUNT(*) FROM pol, pnt WHERE ST_Intersects(pol.geometry,pnt.geometry) GROUP BY pol.id" -dialect SQLITE input.vrt

As usual, you need to have OGR compiled with SQLITE SQL support. For more, check ExecuteSQL page.

Related Question