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

31

Estoy intentando hacer una unión espacial como en el ejemplo aquí: ¿Hay una opción de python para "unir atributos por ubicación"? . Sin embargo, ese enfoque parece realmente ineficiente / lento. Incluso ejecutar esto con unos modestos 250 puntos lleva casi 2 minutos y falla completamente en archivos de forma con> 1,000 puntos. ¿Hay un mejor enfoque? Me gustaría hacer esto completamente en Python sin usar ArcGIS, QGIS, etc.

También me interesaría saber si es posible SUMAR atributos (es decir, población) de todos los puntos que se encuentran dentro de un polígono y unir esa cantidad al archivo de forma del polígono.

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

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

que dice:

TypeError: tipos de operando no admitidos para + =: 'NoneType' y 'float'.

Si reemplazo el "+ =" con "=" funciona bien pero eso no suma los campos. También intenté hacerlos como enteros, pero eso 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']})
jburrfischer
fuente
Creo que deberías editar tu segunda pregunta de aquí para que ésta permanezca enfocada en lo que supongo que es la pregunta más importante para ti. El otro puede investigarse / preguntarse por separado.
PolyGeo

Respuestas:

37

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

Ejemplo de suma de atributos usando las referencias dadas por Mike T:

ingrese la descripción de la imagen aquí

# 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 usar dos métodos, con o sin un í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 R-tree (puede usar 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(multi[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 recorrer en iteración todas las geometrías (polígonos y puntos).
  • Con un índice espacial delimitador ( Spatial Index RTree ), itera solo a través de las geometrías que tienen la oportunidad de cruzarse con su 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 se debe recuperar una gran parte del conjunto de datos, un índice espacial no puede proporcionar ningún beneficio 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 Rtree con OGR, Shapely, Fiona

gene
fuente
15

Además, las geopandas ahora se incluyen opcionalmente rtreecomo una dependencia, consulte el repositorio de Github

Entonces, en lugar de seguir todo el código (muy bueno) anterior, simplemente podría 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 funcionalidad elegante, asegúrese de instalar primero la biblioteca C libspatialindex

EDITAR: importaciones de paquetes corregidos

Claytonrsh
fuente
Tenía la impresión de que rtreeera opcional. ¿No significa eso que necesita instalar rtreetan bien como la libspatialindexbiblioteca C?
kuanb
que ha sido un tiempo, pero creo que cuando hice esto geopandas Instalación desde github añaden automáticamente rtreecuando había instalado primero libspatialindex... que han hecho una versión bastante importante ya que por lo que estoy seguro que las cosas han cambiado un poco
claytonrsh
9

Use Rtree como índice para realizar las uniones mucho más rápidas, 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 sobre 'SUMA', use un dictobjeto para acumular poblaciones usando una identificación de polígono como clave. Aunque, este tipo de cosas se hace mucho mejor con PostGIS.

Mike T
fuente
Gracias @Mike T ... usar el objeto dict o PostGIS son excelentes sugerencias. Sin embargo, todavía estoy un poco confundido sobre dónde puedo incorporar Rtree en mi código (incluido el código anterior).
jburrfischer
1

Esta página web muestra cómo usar una búsqueda de punto en el polígono de Bounding Box antes de la consulta espacial Dentro más costosa de Shapely.

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

klewis
fuente
Gracias @klewis ... ¿alguna posibilidad de que puedas ayudar con la segunda parte? Para resumir los atributos de punto (por ejemplo, población) que se encuentran dentro de los polígonos, intenté algo similar al código a continuación, pero arrojó un error. if shape (school ['geometry']). inside (shape (vecindario ['geometry'])): vecindario ['propiedades'] ['población'] + = escuelas ['propiedades'] ['población']
jburrfischer
Si abre vecindario en modo 'r', podría ser de solo lectura. ¿Ambos shapefiles tienen población de campo? ¿Qué línea # arroja el error? Buena suerte.
klewis
Gracias nuevamente @klewis ... Agregué mi código arriba y expliqué el error. Además, he estado jugando con rtree y todavía estoy un poco confundido sobre dónde agregaría esto en el código anterior. Lamento ser una molestia.
jburrfischer
Pruebe esto, parece que agregar None a un int está causando el error. 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