1

I'm trying to do a fairly simple analysis by buffering points (SEPTA subway stations [below]) and counting the number of 'incidents' (also points) that lie within the buffer. That's it.

I've done some looking online, but I can't find anything specific. Maybe it's so easy that nobody needs to ask. I could really use some help.

I was able to fix the code and create a buffer for the points, but I can't seem to be able to count the points that come within the buffer. Also, there was an issue with the 'incidents' point set that I was using, so I switched it with 'Farmers Markets'. Below is what I've got so far. Again, I just need to count the points.

from osgeo import ogr

septaclip = ogr.Open(r'/home/user/Downloads/SEPTAclip.shp')
septalyr = septaclip.GetLayer(0)
citylimits = ogr.Open(r'/home/user/Downloads/City_Limits.shp')
citylyr = citylimits.GetLayer(0)
crimestat = ogr.Open(r'/home/user/Downloads/Farmers_Markets.shp')
crimelyr = crimestat.GetLayer(0)

memory_driver = ogr.GetDriverByName('Memory')
memory_ds = memory_driver.CreateDataSource('Temp')
buff_lyr = memory_ds.CreateLayer('Buffer')
buff_feat = ogr.Feature(buff_lyr.GetLayerDefn())

multipoly = ogr.Geometry(ogr.wkbMultiPolygon)
for septafeat in septalyr:             
    buff_geo = septafeat.geometry().Buffer(3000) 
    multipoly.AddGeometry(buff_geo)
#multipoly = (multipoly.UnionCascaded())  

for crimefeat in crimelyr:
    buffcrime = crimefeat.geometry().Intersection(multipoly)
Mogsdad
  • 44,709
  • 21
  • 151
  • 275
Caleb Fritz
  • 11
  • 3
  • 1
    I would suggest using fiona and shapely for this. (You could also look at GeoPandas instead of fiona). You can then iterate the points and use Shapely's `intersect()` method to find the overlaps. – songololo Dec 17 '15 at 12:17

1 Answers1

1

Try this for counting the points after intersection:

count = 0
for crimefeat in crimelyr:
    if not crimefeat.geometry().Intersection(multipoly).IsEmpty():
        count += 1

Or using OGR's capabilities:

crime_multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for crimefeat in crimelyr:
    crime_multipoint.AddGeometry(crimefeat.geometry())
crime_multipoint.Intersection(multipoly).GetGeometryCount()
Martin Valgur
  • 5,793
  • 1
  • 33
  • 45
  • Thanks for the help. The first one didn't work, and the second one didn't give me any errors However, it gave me an output of zero, which can't be right. – Caleb Fritz Dec 16 '15 at 03:48