Effizienterer räumlicher Join in Python ohne QGIS, ArcGIS, PostGIS usw

31

Ich versuche, eine räumliche Verknüpfung zu erstellen, ähnlich wie im folgenden Beispiel: Gibt es eine Python-Option zum "Verknüpfen von Attributen nach Speicherort"? . Dieser Ansatz scheint jedoch wirklich ineffizient / langsam zu sein. Selbst das Ausführen mit bescheidenen 250 Punkten dauert fast 2 Minuten und schlägt bei Shapefiles mit> 1.000 Punkten vollständig fehl. Gibt es einen besseren Ansatz? Ich möchte dies vollständig in Python tun, ohne ArcGIS, QGIS usw. zu verwenden.

Es würde mich auch interessieren, ob es möglich ist, Attribute (dh Grundgesamtheit) aller Punkte, die in ein Polygon fallen, zu summieren und diese Menge mit dem Polygon-Shapefile zu verbinden.

Hier ist der Code, den ich konvertieren möchte. Ich erhalte eine Fehlermeldung in Zeile 9:

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

was sagt:

TypeError: Nicht unterstützte Operandentypen für + =: 'NoneType' und 'float'.

Wenn ich das "+ =" durch "=" ersetze, läuft es gut, aber das summiert nicht die Felder. Ich habe auch versucht, diese als ganze Zahlen zu machen, aber das schlägt auch fehl.

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
quelle
Ich denke, Sie sollten Ihre zweite Frage hier herausarbeiten, damit sich diese weiterhin auf die Frage konzentriert, von der ich annehme, dass sie für Sie wichtiger ist. Der andere kann separat recherchiert / angefragt werden.
PolyGeo

Antworten:

37

Fiona gibt Python-Wörterbücher zurück und kann nicht poly['properties']['score']) += point['properties']['score'])mit einem Wörterbuch verwendet werden.

Beispiel für das Summieren von Attributen unter Verwendung der von Mike T angegebenen Referenzen:

Bildbeschreibung hier eingeben

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

Jetzt können wir zwei Methoden mit oder ohne räumlichen Index verwenden:

1) ohne

# 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) mit einem R-Tree-Index (Sie können pyrtree oder rtree verwenden )

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

Ergebnis mit den beiden Lösungen:

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

Was ist der Unterschied ?

  • Ohne den Index müssen Sie alle Geometrien (Polygone und Punkte) durchlaufen.
  • Mit einem begrenzenden räumlichen Index ( Spatial Index RTree ) iterieren Sie nur durch die Geometrien, die sich mit Ihrer aktuellen Geometrie überschneiden können ('Filter', wodurch eine erhebliche Menge an Berechnungen und Zeit gespart werden kann ...).
  • Ein Raumindex ist jedoch kein Zauberstab. Wenn ein sehr großer Teil des Datensatzes abgerufen werden muss, kann ein Raumindex keinen Geschwindigkeitsvorteil bieten.

Nach:

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

Weitere Informationen finden Sie unter Verwenden der räumlichen Indexierung von Rtree mit OGR, Shapely, Fiona

Gen
quelle
15

Zusätzlich - Geopandas enthält nun optional rtreeals Abhängigkeit das Github-Repo

Anstatt dem obigen (sehr netten) Code zu folgen, können Sie einfach Folgendes tun:

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

Um diese schicke Funktionalität zu erhalten, müssen Sie zuerst die C-Bibliothek libspatialindex installieren

EDIT: Korrigierte Paketimporte

claytonrsh
quelle
Ich hatte den Eindruck rtreeoptional zu sein. Bedeutet das nicht, dass Sie rtreeneben der libspatialindexC-Bibliothek auch die Installation durchführen müssen ?
Kuanb
es ist ein gewesen , während , aber ich denke , wenn ich diese Installation geopandas von GitHub habe automatisch hinzugefügt , rtreewenn ich zum ersten Mal installiert hatte libspatialindex... sie eine ziemlich große Version getan hat , da so ich habe sicher die Dinge ein wenig verändert
claytonrsh
9

Verwenden Sie Rtree als Index, um die viel schnelleren Verknüpfungen durchzuführen, und dann Shapely, um die räumlichen Prädikate zu erstellen , um festzustellen, ob sich ein Punkt tatsächlich innerhalb eines Polygons befindet. Bei richtiger Ausführung kann dies schneller sein als bei den meisten anderen GIS.

Beispiele finden Sie hier oder hier .

Im zweiten Teil Ihrer Frage zu 'SUM' verwenden Sie ein dictObjekt, um Populationen unter Verwendung einer Polygon-ID als Schlüssel zu akkumulieren. Mit PostGIS ist dies jedoch wesentlich einfacher.

Mike T
quelle
Danke @Mike T ... mit dem Diktierobjekt oder PostGIS sind tolle Vorschläge. Ich bin immer noch ein wenig verwirrt, wo ich Rtree in meinen Code einbauen kann (Code oben enthalten).
jburrfischer
1

Diese Webseite zeigt, wie eine Begrenzungsrahmen-Point-in-Polygon-Suche vor der teureren räumlichen Abfrage von Shapely verwendet wird.

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

klewis
quelle
Danke @klewis ... hast du eine Chance, beim zweiten Teil zu helfen? Um die Punktattribute (z. B. Grundgesamtheit) zu summieren, die in die Polygone fallen, habe ich etwas Ähnliches wie den folgenden Code versucht, aber es wurde ein Fehler ausgegeben. if shape (schule ['geometrie']). within (form (nachbarschaft ['geometrie']): nachbarschaft ['eigenschaften'] ['bevölkerung'] + = schulen ['
eigenschaften
Wenn Sie die Nachbarschaft im "r" -Modus öffnen, ist sie möglicherweise schreibgeschützt. Haben beide Shapefiles eine Feldpopulation? Welche Zeile # löst den Fehler aus? Viel Glück.
Klewis
Nochmals vielen Dank @klewis ... Ich habe meinen Code oben hinzugefügt und den Fehler erklärt. Außerdem habe ich mit rtree rumgespielt und bin immer noch ein bisschen verwirrt, wo ich das in den obigen Code einfügen würde. Tut mir leid, dass ich so nervig bin.
jburrfischer
Versuchen Sie dies, scheint das Hinzufügen von None zu einem Int den Fehler zu verursachen. 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