Как подсчитать количество полигонов, которые пересекает фигура? - PullRequest
0 голосов
/ 10 января 2019

У меня очень большой набор данных с полигонами и точками с буферами вокруг них. Я хотел бы создать новый столбец в данных точек, который включает количество полигонов, которые пересекает буфер точек.

Вот упрощенный пример:

import pandas as pd
import geopandas as gp
from shapely.geometry import Polygon
from shapely.geometry import Point
import matplotlib.pyplot as plt

## Create polygons and points ##

df = gp.GeoDataFrame([['a',Polygon([(1, 0), (1, 1), (2,2), (1,2)])],
                     ['b',Polygon([(1, 0.25), (2,1.25), (3,0.25)])]],
                     columns = ['name','geometry'])
df = gp.GeoDataFrame(df, geometry = 'geometry')

points = gp.GeoDataFrame( [['box', Point(1.5, 1.115), 4],
                        ['triangle', Point(2.5,1.25), 8]],
                     columns=['name', 'geometry', 'value'], 
                     geometry='geometry')

##Set a buffer around the points##

buf = points.buffer(0.5)
points['buffer'] = buf
points = points.drop(['geometry'], axis = 1) 
points = points.rename(columns = {'buffer': 'geometry'})

Эти данные выглядят так: enter image description here Что я хотел бы сделать, так это создать еще один столбец в кадре данных точек, содержащий количество полигонов, которые пересекает точка.

Я попытался использовать цикл for следующим образом:

points['intersect'] = []
for geo1 in points['geometry']:
    for geo2 in df['geometry']:
        if geo1.intersects(geo2):
            points['intersect'].append('1')

Который я бы затем суммировал, чтобы получить общее количество пересечений. Однако я получаю ошибку: «Длина значений не соответствует длине индекса». Я знаю, что это потому, что он пытается назначить три строки данных для кадра только с двумя строками.

Как я могу объединить подсчеты, чтобы первой точке было присвоено значение 2, а второй - 1?

1 Ответ

0 голосов
/ 10 января 2019

Если у вас большой набор данных, я бы выбрал решение с использованием пространственного индекса rtree, что-то вроде этого.

import pandas as pd
import geopandas as gp
from shapely.geometry import Polygon
from shapely.geometry import Point
import matplotlib.pyplot as plt

## Create polygons and points ##

df = gp.GeoDataFrame([['a',Polygon([(1, 0), (1, 1), (2,2), (1,2)])],
                     ['b',Polygon([(1, 0.25), (2,1.25), (3,0.25)])]],
                     columns = ['name','geometry'])
df = gp.GeoDataFrame(df, geometry = 'geometry')

points = gp.GeoDataFrame( [['box', Point(1.5, 1.115), 4],
                        ['triangle', Point(2.5,1.25), 8]],
                     columns=['name', 'geometry', 'value'], 
                     geometry='geometry')

# generate spatial index
sindex = df.sindex
# define empty list for results
results_list = []
# iterate over the points
for index, row in points.iterrows():
    buffer = row['geometry'].buffer(0.5)  # buffer
    # find approximate matches with r-tree, then precise matches from those approximate ones
    possible_matches_index = list(sindex.intersection(buffer.bounds))
    possible_matches = df.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(buffer)]
    results_list.append(len(precise_matches))
# add list of results as a new column
points['polygons'] = pd.Series(results_list)
...