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

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

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:

# 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

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

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.

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/