Uniones espaciales más eficientes en Python sin QGIS, ArcGIS, PostGIS, etc.

27

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

También me gustarí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 archivo de forma de polígono.

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

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

que dice:

  

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

Si sustituyo el "+=" por "=" funciona bien, pero eso no suma los campos. También he intentado hacer estos 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']})
    
pregunta jburrfischer 23.06.2014 - 20:37

4 respuestas

35

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

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

#readtheshapefilesimportfionafromshapely.geometryimportshapepolygons=[polforpolinfiona.open('poly.shp')]points=[ptforptinfiona.open('point.shp')]#attributesofthepolygonsforpolyinpolygons:printpoly['properties']OrderedDict([(u'score',0)])OrderedDict([(u'score',0)])OrderedDict([(u'score',0)])#attributesofthepointsforptinpoints:printi['properties']OrderedDict([(u'score',1)])....#(sameforthe8points)

Ahora,podemosusardosmétodos,conosinuníndiceespacial:

1)sin

#iteratethroughpointsfori,ptinenumerate(points):point=shape(pt['geometry'])#iteratethroughpolygonsforj,polyinenumerate(polygons):ifpoint.within(shape(poly['geometry'])):#sumofattributesvaluespolygons[j]['properties']['score']=polygons[j]['properties']['score']+points[i]['properties']['score']

2)conun índice de árbol R (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 todas las geometrías (polígonos y puntos).
  • Con un índice espacial delimitador ( RTatial Index RTree ), solo se itera las geometrías que tienen la posibilidad de intersectarse 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, consulte Uso de la indexación espacial de Rtree con OGR , Shapely, Fiona

    
respondido por el gene 24.06.2014 - 21:30
10

Además, las geopandas ahora incluyen opcionalmente rtree como una dependencia, consulte el github repo

Entonces, en lugar de seguir todo el código (muy bueno) de arriba, simplemente puedes 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 la biblioteca C libspatialindex primero

EDITAR: importaciones de paquetes corregidos

    
respondido por el claytonrsh 06.10.2015 - 05:01
9

Utilice Rtree como un índice para realizar uniones mucho más rápidas, luego Shapely para hacer los predicados espaciales para determinar si un punto es 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.

Consulte los ejemplos aquí o aquí .

En la segunda parte de su pregunta relativa a 'SUM', use un objeto dict para acumular poblaciones usando una identificación de polígono como clave. Aunque, este tipo de cosas se hace mucho más bien con PostGIS.

    
respondido por el Mike T 23.06.2014 - 23:50
1

Esta página web muestra cómo utilizar una búsqueda de punto en polígono de Cuadro delimitador antes de la más costosa dentro de la consulta espacial de Shapely.

enlace

    
respondido por el klewis 23.06.2014 - 20:53

Lea otras preguntas en las etiquetas