Python (возможно, Shapely) для создания буфера и подсчета очков

Я пытаюсь провести довольно простой анализ, буферизуя точки (станции метро SEPTA [ниже]) и подсчитывая количество «происшествий» (также точек), которые находятся в буфере. Вот и все.

Я кое-что поискал в Интернете, но ничего конкретного не нашел. Может быть, это так просто, что никто не должен спрашивать. Мне действительно не помешала бы помощь.

Мне удалось исправить код и создать буфер для точек, но я не могу подсчитать точки, попадающие в буфер. Кроме того, возникла проблема с набором точек «инциденты», который я использовал, поэтому я переключил его на «Фермерские рынки». Ниже то, что у меня есть до сих пор. Опять же, мне просто нужно подсчитать очки.

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)

python buffer intersection shapely ogr
person Caleb Fritz    schedule 16.12.2015    source источник
comment
Я бы предложил использовать fiona и shapely для этого. (Вы также можете посмотреть на GeoPandas вместо fiona). Затем вы можете повторить точки и использовать метод Shapely intersect(), чтобы найти перекрытия.   -  person songololo    schedule 17.12.2015

Ответы (1)


arrow_upward
1
arrow_downward

Попробуйте это для подсчета точек после пересечения:

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

Или используя возможности OGR:

crime_multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for crimefeat in crimelyr:
    crime_multipoint.AddGeometry(crimefeat.geometry())
crime_multipoint.Intersection(multipoly).GetGeometryCount()
person Martin Valgur    schedule 16.12.2015
comment
Спасибо за помощь. Первый не сработал, а второй не дал мне никаких ошибок. Однако он дал мне нулевой результат, что не может быть правильным. - person Caleb Fritz; 16.12.2015