Junção espacial mais eficiente em Python sem QGIS, ArcGIS, PostGIS, etc

Estou tentando fazer uma junção espacial muito parecida com o exemplo aqui: Existe uma opção python para "juntar atributos por local"? . No entanto, essa abordagem parece realmente ineficiente / lenta. Mesmo executando isso com 250 pontos modestos leva quase 2 minutos e falha totalmente em arquivos de forma com mais de 1.000 pontos. Existe uma abordagem melhor? Eu gostaria de fazer isso inteiramente em Python sem usar ArcGIS, QGIS, etc.

Também gostaria de saber se é possível adicionar atributos SUM (ou seja, população) de todos os pontos que se enquadram em um polígono e associar essa quantidade ao shapefile do polígono.

Aqui está o código que estou tentando converter. Eu recebo um erro na linha 9:

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

que diz:

TypeError: tipo (s) de operando não suportado para + =: 'NoneType' e 'float'.

Se eu substituir o "+ =" por "=", ele funciona bem, mas isso não soma os campos. Eu também tentei fazer isso como números inteiros, mas isso também falha.

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']})

Fiona retorna dicionários Python e você não pode usar poly['properties']['score']) += point['properties']['score'])com um dicionário.

Exemplo de soma de atributos usando as referências fornecidas por Mike T:

# 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)

Agora, podemos usar dois métodos, com ou sem um índice espacial:

1) sem

# 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) com um índice R-tree (você pode usar pyrtree ou 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']

Resultado com as duas soluções:

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

Qual é a diferença ?

• Sem o índice, você deve percorrer todas as geometrias (polígonos e pontos).
• Com um índice espacial delimitador ( Spatial Index RTree ), você itera apenas pelas geometrias que têm a chance de se cruzar com a geometria atual ('filtro' que pode economizar uma quantidade considerável de cálculos e tempo ...)
• mas um índice espacial não é uma varinha mágica. Quando uma parte muito grande do conjunto de dados precisa ser recuperada, um Índice Espacial não pode oferecer nenhum benefício de velocidade.

Depois de:

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

Para ir além, consulte Usando a indexação espacial da Rtree com OGR, Shapely, Fiona

Além disso - as geopandas agora incluem opcionalmente rtreecomo dependência, consulte o repositório do github

Então, em vez de seguir todo o código (muito bom) acima, você pode simplesmente fazer algo como:

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'])

Para obter essa funcionalidade sofisticada, instale primeiro a biblioteca C libspatialindex

EDIT: importação de pacotes corrigidos

Use Rtree como um índice para executar as junções muito mais rápidas e, em seguida, Shapely para executar os predicados espaciais para determinar se um ponto está realmente dentro de um polígono. Se feito corretamente, isso pode ser mais rápido do que a maioria dos outros GIS.

Veja exemplos aqui ou aqui .

A segunda parte da sua pergunta sobre 'SUM', use um dictobjeto para acumular populações usando um ID de polígono como chave. Embora, esse tipo de coisa seja feito muito mais bem com o PostGIS.

Esta página da web mostra como usar uma pesquisa de ponto no polígono da caixa delimitadora antes da consulta espacial Within mais cara da Shapely.

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