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

31

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']})
jburrfischer
fonte
Eu acho que você deve editar sua segunda pergunta daqui para que esta permaneça focada no que eu assumo ser a pergunta mais importante para você. O outro pode ser pesquisado / solicitado separadamente.
PolyGeo

Respostas:

37

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:

insira a descrição da imagem aqui

# 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

gene
fonte
15

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

claytonrsh
fonte
Eu estava com a impressão de que rtreeera opcional. Isso não significa que você precisa instalar rtree, assim como a libspatialindexbiblioteca C?
kuanb
já faz um tempo, mas acho que, quando eu fiz isso, a instalação de geopandas do github foi adicionada automaticamente rtreequando eu o instalei libspatialindex... eles fizeram uma versão bastante importante desde então, então eu tenho certeza que as coisas mudaram um pouco
claytonrsh
9

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.

Mike T
fonte
Obrigado @ Mike T ... usando o objeto dict ou PostGIS são ótimas sugestões. No entanto, ainda estou um pouco confuso sobre onde posso incorporar o Rtree no meu código (código incluído acima).
Jburrfischer
1

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/

klewis
fonte
Obrigado @klewis ... tem alguma chance de ajudar com a segunda parte? Para somar os atributos de ponto (por exemplo, população) que se enquadram nos polígonos, tentei algo semelhante ao código abaixo, mas ocorreu um erro. if shape (escola ['geometria']). within (shape (vizinhança ['geometria'])): bairro ['propriedades'] ['população'] + = escolas ['propriedades'] ['população']
jburrfischer
Se você abrir um bairro no modo 'r', ele poderá ser somente leitura. Os dois shapefiles possuem população de campo? Qual linha # está lançando o erro? Boa sorte.
klewis
Mais uma vez obrigado @klewis ... Adicionei meu código acima e expliquei o erro. Além disso, eu tenho brincado com o rtree e ainda estou pouco confuso sobre onde eu adicionaria isso no código acima. Desculpe ser um incômodo.
Jburrfischer
Tente isso, parece que adicionar Nenhum a um int está causando o erro. 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 24/06