Una unión espacial más eficiente en Python sin QGIS, ArcGIS, PostGIS, etc.

Question

Estoy intentando hacer una unión espacial muy parecida al ejemplo de aquí: ¿Existe una opción en python para "unir atributos por ubicación"? . Sin embargo, ese enfoque parece realmente ineficiente / lento. Incluso ejecutando esto con unos modestos 250 puntos se tarda casi 2 minutos y falla por completo en shapefiles con > 1.000 puntos. ¿Existe un enfoque mejor? Me gustaría hacer esto completamente en Python sin usar ArcGIS, QGIS, etc.

También me interesaría saber si es posible sumar los atributos (es decir, la población) de todos los puntos que caen dentro de un polígono y unir esa cantidad al shapefile del polígono.

Aquí está el código que estoy tratando de convertir. Me da un error en la línea 9:

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

que dice:

TypeError: tipo(s) de operando no soportado(s) para +=: 'NoneType' y 'float'.

Si sustituyo el "+=" por "=" funciona bien pero no suma los campos. También he probado a hacerlos como enteros pero también falla.

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

Answer 1

Fiona devuelve diccionarios de Python y no se puede utilizar poly['properties']['score']) += point['properties']['score']) con un diccionario.

Ejemplo de suma de atributos utilizando las referencias dadas 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)

Ahora, podemos utilizar dos métodos, con o sin índice espacial:

1: sin

# 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: con un Índice del árbol R (puede utilizar pyrtree o 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(polygons[j]['geometry'])):
            polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']

Resultado con las dos soluciones:

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

¿Cuál es la diferencia?

• Sin el índice, debe iterar por todas las geometrías (polígonos y puntos).
• Con un índice espacial delimitador ( Índice espacial RTree ), se itera sólo a través de las geometrías que tienen la posibilidad de intersecarse con la geometría actual ("filtro" que puede ahorrar una cantidad considerable de cálculos y tiempo...)
• pero un Índice Espacial no es una varita mágica. Cuando hay que recuperar una parte muy grande del conjunto de datos, un índice espacial no puede aportar ninguna ventaja de velocidad.

Después:

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 más lejos, mira Uso de la indexación espacial de Rtree con OGR, Shapely, Fiona

Answer 2

Además - geopandas ahora incluye opcionalmente rtree como una dependencia, vea el repo de github

Así que, en lugar de seguir todo el código (muy bonito) de arriba, podrías simplemente hacer 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 obtener esta elegante funcionalidad, asegúrese de instalar la biblioteca C libspatialindex primero

EDITADO: se ha corregido la importación de paquetes

Answer 3

Utilice Rtree como índice para realizar las uniones mucho más rápidas, y luego Shapely para hacer los predicados espaciales para determinar si un punto está realmente dentro de un polígono. Si se hace correctamente, esto puede ser más rápido que la mayoría de los otros SIG.

Ver ejemplos aquí o aquí .

La segunda parte de su pregunta relativa a la "SUMA", utiliza un dict para acumular poblaciones utilizando un id de polígono como clave. Aunque, este tipo de cosas se hace mucho mejor con PostGIS.

Answer 4

Esta página web muestra cómo utilizar una búsqueda de puntos en polígonos de Bounding Box antes de la más costosa consulta espacial Within de Shapely.

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