Более эффективное пространственное соединение в Python без QGIS, ArcGIS, PostGIS и т. Д.

31

Я пытаюсь выполнить пространственное соединение, очень похожее на приведенный здесь пример: есть ли опция Python для «объединения атрибутов по местоположению»? , Однако такой подход кажется действительно неэффективным / медленным. Даже выполнение этого со скромными 250 точками занимает почти 2 минуты, и оно полностью не работает на шейп-файлах с> 1000 точками. Есть ли лучший подход? Я хотел бы сделать это полностью в Python без использования ArcGIS, QGIS и т. Д.

Мне также было бы интересно узнать, возможно ли суммировать атрибуты (то есть совокупность) всех точек, которые попадают в полигон, и присоединять это количество к шейп-файлу полигона.

Вот код, который я пытаюсь преобразовать. Я получаю сообщение об ошибке в строке 9:

poly['properties']['score'] += point['properties']['score']

который говорит:

Ошибка типа: неподдерживаемые типы операндов для + =: 'NoneType' и 'float'.

Если я заменю «+ =» на «=», он будет работать нормально, но это не суммирует поля. Я также пытался сделать их как целые числа, но это тоже не получается.

with fiona.open(poly_shp, 'r') as n: 
  with fiona.open(point_shp,'r') as s:
    outSchema = {'geometry': 'Polygon','properties':{'region':'str','score':'float'}}
    with fiona.open (out_shp, 'w', 'ESRI Shapefile', outSchema, crs) as output:
        for point in s:
            for poly in n:
                if shape(point['geometry']).within(shape(poly['geometry'])):  
                    poly['properties']['score']) += point['properties']['score'])
                    output.write({
                        'properties':{
                            'region':poly['properties']['NAME'],
                            'score':poly['properties']['score']},
                        'geometry':poly['geometry']})
jburrfischer
источник
Я думаю, что вы должны отредактировать свой второй вопрос отсюда, чтобы этот вопрос был сосредоточен на том, что, как я полагаю, является для вас более важным вопросом. Другой может быть исследован / задан отдельно.
PolyGeo

Ответы:

37

Фиона возвращает словари Python, и вы не можете использовать poly['properties']['score']) += point['properties']['score'])со словарем.

Пример суммирования атрибутов с использованием ссылок, данных Майком Т:

введите описание изображения здесь

# read the shapefiles 
import fiona
from shapely.geometry import shape
polygons = [pol for pol in fiona.open('poly.shp')]
points = [pt for pt in fiona.open('point.shp')]
# attributes of the polygons
for poly in polygons:
   print poly['properties'] 
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])

# attributes of the points
for pt in points:
    print i['properties']
 OrderedDict([(u'score', 1)]) 
 .... # (same for the 8 points)

Теперь мы можем использовать два метода, с или без пространственного индекса:

1) без

# iterate through points 
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     #iterate through polygons
     for j, poly in enumerate(polygons):
        if point.within(shape(poly['geometry'])):
             # sum of attributes values
             polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

2) с индексом R-дерева (вы можете использовать pyrtree или rtree )

# Create the R-tree index and store the features in it (bounding box)
 from rtree import index
 idx = index.Index()
 for pos, poly in enumerate(polygons):
       idx.insert(pos, shape(poly['geometry']).bounds)

#iterate through points
for i,pt in enumerate(points):
  point = shape(pt['geometry'])
  # iterate through spatial index
  for j in idx.intersection(point.coords[0]):
      if point.within(shape(multi[j]['geometry'])):
            polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

Результат с двумя решениями:

for poly in polygons:
   print poly['properties']    
 OrderedDict([(u'score', 2)]) # 2 points in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon
 OrderedDict([(u'score', 1)]) # 1 point in the polygon

В чем разница ?

  • Без индекса вы должны пройти через все геометрии (многоугольники и точки).
  • С помощью ограничивающего пространственного индекса ( Spatial Index RTree ) вы перебираете только те геометрии, которые могут пересекаться с вашей текущей геометрией («фильтр», который может сэкономить значительное количество вычислений и времени ...)
  • но Пространственный Индекс не волшебная палочка. Когда требуется извлечь очень большую часть набора данных, пространственный индекс не может дать никакого преимущества в скорости.

После:

schema = fiona.open('poly.shp').schema
with fiona.open ('output.shp', 'w', 'ESRI Shapefile', schema) as output:
    for poly in polygons:
        output.write(poly)

Чтобы продолжить, посмотрите на Использование пространственной индексации Rtree с OGR, Shapely, Fiona

ген
источник
15

Кроме того - геопанда теперь необязательно включает rtreeв качестве зависимости, см. Репозиторий github

Таким образом, вместо того, чтобы следовать всему (очень хорошему) коду выше, вы можете просто сделать что-то вроде

import geopandas
from geopandas.tools import sjoin
point = geopandas.GeoDataFrame.from_file('point.shp') # or geojson etc
poly = geopandas.GeoDataFrame.from_file('poly.shp')
pointInPolys = sjoin(point, poly, how='left')
pointSumByPoly = pointInPolys.groupby('PolyGroupByField')['fields', 'in', 'grouped', 'output'].agg(['sum'])

Чтобы получить эту привлекательную функциональность, сначала установите C-библиотеку libspatialindex.

РЕДАКТИРОВАТЬ: исправлен импорт пакетов

claytonrsh
источник
Я был под впечатлением rtreeнеобязательно. Разве это не означает, что вам нужно установить rtreeтак же, как libspatialindexC-библиотеку?
17
Прошло много времени, но я думаю, что когда я делал это, установка геопанд из github автоматически добавлялась, rtreeкогда я только что установил libspatialindex... с тех пор они сделали довольно серьезную версию, так что я уверен, что все немного изменилось
claytonrsh
9

Используйте Rtree в качестве индекса для выполнения гораздо более быстрых объединений, затем Shapely для выполнения пространственных предикатов, чтобы определить, действительно ли точка находится внутри многоугольника. Если все сделано правильно, это может быть быстрее, чем большинство других ГИС.

Смотрите примеры здесь или здесь .

Вторая часть вашего вопроса о «СУММЕ» - используйте dictобъект для накопления населения, используя идентификатор многоугольника в качестве ключа. Хотя с PostGIS этот тип вещей гораздо удобнее.

Майк Т
источник
Спасибо @Mike T ... отличные предложения - использование объекта dict или PostGIS. Однако я все еще немного сбит с толку, когда я могу включить Rtree в свой код (включая код выше).
jburrfischer
1

На этой веб-странице показано, как использовать поиск точки-полигона в ограничивающей рамке перед более дорогим пространственным запросом в Shapely.

http://rexdouglass.com/fast-spatial-joins-in-python-with-a-spatial-index/

klewis
источник
Спасибо @klewis ... есть ли шанс, что вы можете помочь со второй частью? Чтобы суммировать точечные атрибуты (например, совокупность), которые попадают в полигоны, я попытался сделать что-то похожее на код ниже, но он выдал ошибку. if shape (школа ['геометрия']). внутри (форма (окрестность ['геометрия'])): окрестность ['свойства'] ['население'] + = школы ['свойства'] ['население']
jburrfischer
Если вы открываете соседство в режиме 'r', оно может быть доступно только для чтения. Имеются ли оба шейп-файла в поле? Какая строка # выдает ошибку? Удачи.
Klewis
Еще раз спасибо @klewis ... Я добавил свой код выше и объяснил ошибку. Кроме того, я поигрался с rtree, и я все еще немного растерялся, когда бы добавил это в код выше. Извините за беспокойство.
jburrfischer
Попробуйте это, кажется, добавление None к int вызывает ошибку. poly_score = poly ['properties'] ['score']) point_score = point ['properties'] ['score']) if point_score: if poly_score poly ['properties'] ['score']) + = point_score else: poly ['properties'] ['score']) = point_score
klewis