Как реализовать дерево для расчета площади пересечения? - PullRequest
Новая
       44

Как реализовать дерево для расчета площади пересечения?

2 голосов
/ 07 мая 2019

Вот некоторые примеры данных: 

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

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

df = gp.GeoDataFrame(df, geometry = 'geometry')
df['area'] = df.area

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

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

И это выглядит так:

enter image description here

В основном я пытаюсьнайти область пересечения между этими объектами.

До сих пор я делал это с помощью этого кода: 

data = []
for index, geo in df.iterrows():
    for index2, poin in points.iterrows():
        if geo['geometry'].intersects(poin['geometry']):
          data.append({'geometry':geo['geometry'].intersection(poin['geometry']), 'area': geo['geometry'].intersection(poin['geometry']).area})

df2 = gp.GeoDataFrame(data, columns = ['geometry', 'area'])

Однако реальные данные, которые я буду использовать для этого, имеют 100 000 с.полигоны, поэтому этот код будет невероятно трудоемким.Я знаю, что мог бы ускорить это с помощью r-деревьев.Однако я не могу реализовать это должным образом.

Я пробовал что-то вроде этого: 

spatial_index = df.sindex
results_list = []
for index, row in points.iterrows():
    buffer = row['geometry']
    possible_matches_index = list(spatial_index.intersection(buffer.bounds)) 
    possible_matches = df.iloc[possible_matches_index]
    results_list.append({'geometry':possible_matches['geometry'].intersection(row['geometry']), 'area': possible_matches['geometry'].intersection(row['geometry']).area})

df = gp.GeoDataFrame(results_list, columns = ['geometry', 'area'])

Но это помещает все пересечения для каждого квадрата в одну линию.

    geometry                                        area
0   name a POLYGON ((1.615 1.615, 1 1, 0.7 1, 0...  name a 0.000037 b 0.000003 dtype: float64
1                                    name a ...     name a 0.000000 b 0.000013 dtype: float64

Как я могу получить это для создания кадра данных с одной линией для каждого пересечения со столбцами для его геометрии и области?

1 Ответ

1 голос
/ 07 мая 2019

Чтобы рассчитать пересечения между двумя GeoDataFrames, вы можете использовать функцию geopandas.overlay: 

geopandas.overlay(df, points, how='intersection')

Это использует пространственный индекс rtree под капотом (поэтому он должен быть более эффективным, чем грубая сила, удваивающая цикл for), и будет возвращать пересечения всех комбинаций геометрий обоих наборов данных в качестве нового GeoDataFrame ( для которого вы можете затем рассчитать площади).

См. https://geopandas.readthedocs.io/en/latest/set_operations.html для документации по этому вопросу.

...