2

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

Improve this question
4
  • Does must use mean that it is your home work? Myself I would import both shapefiles into PostGIS or Spatialite and do the analysis in the database. OpenJUMP has also dedicated tools for such tasks, probably QGIS as well.
    – user30184
    Oct 28, 2015 at 20:35
  • 1
    I would prefer to use those libraries because this analysis is part of a larger tool that is based largely on these modules. I lack the knowledge in PostGIS and Spatialite to successfully write this query. The closest I can think of is something like this: ogr2ogr -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM points.shp A, polygons.shp B WHERE ST_Intersects(A.geometry, B.geometry);
    – geoeye
    Oct 28, 2015 at 21:23
  • 1
    I was asking if you plan to test some theory or to do productive work because I am pretty sure that joining two datasets with GDAL can't offer support for using spatial index and therefore you will make Cartesian product en.wikipedia.org/wiki/Cartesian_product. With hundreds times tens of thousands operations you will get millions operations. Million times millisecond makes 1000 seconds = about fifteen minutes for each millisecond you need for one operation.
    – user30184
    Oct 28, 2015 at 22:19
  • I am trying to do productive work, not theory. I'd welcome your suggestions to improve processing speed when working with large sets of vector data. I'm quite experienced in ArcGIS but have only been scripting in python/gdal for ~6 months. I am building tools that can count points within polygons as well as intersect polygons and aggregate fields in the new intersection. However the data I am using is so large that these operations frequently take many many hours. Presumably due to the cartesian product problem that you have shed light on. VRT seems to be a way around this, do you know others?
    – geoeye
    Nov 9, 2015 at 21:43

1 Answer 1

Reset to default
5

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.

Improve this answer
7
  • Thanks for your response. I believe I have OGR compiled with SQLITE support, as I've successfully run other sqlite queries. However, I followed your script exactly and renamed my points/polygons accordingly. However, I get the error Parse error at EOF, not all elements have been closed, starting with OGRVRTDataSource All files are residing within the same directory and I have specified absolute paths to each of the datasources. Any ideas what i am missing here?!?
    – geoeye
    Nov 9, 2015 at 21:17
  • I also get Unable to open datasource /path/to/input.vrt' with the following drivers.' after the first error
    – geoeye
    Nov 9, 2015 at 21:22
  • 1
    EDIT: Works when you add the line </OGRCRTDataSource> at the end of input.vrt.
    – geoeye
    Nov 9, 2015 at 21:39
  • 1
    oops, I swear it was there, but got eaten because indentation was wrong) now it's fixed
    – Maxim Dubinin
    Nov 9, 2015 at 22:21
  • Can I ask what GROUP BY pol.id is doing in this instance and how it should be used?
    – geoeye
    Nov 10, 2015 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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