Jointure spatiale plus efficace en Python sans QGIS, ArcGIS, PostGIS, etc.

J'essaie de faire une jointure spatiale un peu comme l'exemple ici: Existe - t-il une option python pour "joindre les attributs par emplacement"? . Cependant, cette approche semble vraiment inefficace / lente. Même le faire avec un modeste 250 points prend presque 2 minutes et il échoue entièrement sur les fichiers de formes avec> 1000 points. Est-ce qu'il y a une meilleure approche? Je voudrais le faire entièrement en Python sans utiliser ArcGIS, QGIS, etc.

Je serais également intéressé de savoir s'il est possible de sommer les attributs (c'est-à-dire la population) de tous les points qui se trouvent dans un polygone et de joindre cette quantité au fichier de formes du polygone.

Voici le code que j'essaie de convertir. Je reçois une erreur sur la ligne 9:

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

qui dit:

TypeError: type (s) d'opérande non pris en charge pour + =: 'NoneType' et 'float'.

Si je remplace le "+ =" par "=" cela fonctionne bien mais cela ne résume pas les champs. J'ai également essayé de les faire sous forme d'entiers, mais cela échoue également.

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 renvoie des dictionnaires Python et vous ne pouvez pas les utiliser poly['properties']['score']) += point['properties']['score'])avec un dictionnaire.

Exemple de sommation d'attributs à l'aide des références données par 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)

Maintenant, nous pouvons utiliser deux méthodes, avec ou sans index spatial:

1) sans

# 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) avec un index R-tree (vous pouvez utiliser pyrtree ou 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']

Résultat avec les deux solutions:

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

Quelle est la différence ?

• Sans l'index, vous devez parcourir toutes les géométries (polygones et points).
• Avec un index spatial englobant ( Spatial Index RTree ), vous itérez uniquement à travers les géométries qui ont une chance de se croiser avec votre géométrie actuelle (`` filtre '' qui peut économiser une quantité considérable de calculs et de temps ...)
• mais un indice spatial n'est pas une baguette magique. Lorsqu'une très grande partie de l'ensemble de données doit être récupérée, un index spatial ne peut pas apporter d'avantage de vitesse.

Aprè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)

Pour aller plus loin, consultez Utilisation de l'indexation spatiale Rtree avec OGR, Shapely, Fiona

De plus - les géopandas incluent désormais en option rtreecomme dépendance, voir le dépôt github

Ainsi, au lieu de suivre tout le code (très agréable) ci-dessus, vous pouvez simplement faire quelque chose comme:

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

Pour obtenir cette fonctionnalité élégante, assurez-vous d'installer la bibliothèque C libspatialindex premier

EDIT: importations de paquets corrigées

Utilisation Rtree comme index pour effectuer les jointures beaucoup plus rapides, puis Shapely pour effectuer les prédicats spatiaux afin de déterminer si un point se trouve réellement dans un polygone. Si cela est fait correctement, cela peut être plus rapide que la plupart des autres SIG.

Voir des exemples ici ou ici .

La deuxième partie de votre question concernant 'SUM', utilisez un dictobjet pour accumuler des populations en utilisant un identifiant de polygone comme clé. Cependant, ce type de chose se fait beaucoup mieux avec PostGIS.

Cette page Web montre comment utiliser une recherche ponctuelle dans un cadre de sélection avant la requête spatiale Within plus coûteuse de Shapely.

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