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

31

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']})
jburrfischer
la source
Je pense que vous devriez éditer votre deuxième question d'ici pour que celle-ci reste concentrée sur ce que je suppose être la question la plus importante pour vous. L'autre peut être recherché / demandé séparément.
PolyGeo

Réponses:

37

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:

entrez la description de l'image ici

# 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

gène
la source
15

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

claytonrsh
la source
J'avais l'impression que rtreec'était facultatif. Cela ne signifie-t-il pas que vous devez installer rtreeainsi que la libspatialindexbibliothèque C?
kuanb
ça fait un moment mais je pense que lorsque j'ai fait cette installation, les géopandas de github ont été ajoutées automatiquement rtreelors de ma première installation libspatialindex... ils ont fait une version assez importante depuis, donc je suis sûr que les choses ont un peu changé
claytonrsh
9

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.

Mike T
la source
Merci @Mike T ... en utilisant l'objet dict ou PostGIS sont d'excellentes suggestions. Je suis encore un peu confus quant à l'endroit où je peux incorporer Rtree dans mon code (code inclus ci-dessus).
jburrfischer
1

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/

klewis
la source
Merci @klewis ... avez-vous une chance de nous aider avec la deuxième partie? Pour résumer les attributs ponctuels (par exemple, la population) qui se trouvent dans les polygones, j'ai essayé quelque chose de similaire au code ci-dessous, mais cela a généré une erreur. si forme (école ['géométrie']). dans (forme (quartier ['géométrie'])): quartier ['propriétés'] ['population'] + = écoles ['propriétés'] ['population']
jburrfischer
Si vous ouvrez le quartier en mode «r», il peut être en lecture seule. Les deux fichiers de formes ont-ils une population sur le terrain? Quelle ligne # génère l'erreur? Bonne chance.
klewis
Merci encore @klewis ... J'ai ajouté mon code ci-dessus et expliqué l'erreur. De plus, j'ai joué avec rtree et je ne suis toujours pas confus quant à l'ajout de cela dans le code ci-dessus. Désolé d'être un tel ennui.
jburrfischer
Essayez ceci, semble ajouter Aucun à un int est à l'origine de l'erreur. poly_score = poly ['propriétés'] ['score']) point_score = point ['propriétés'] ['score']) si point_score: si poly_score poly ['propriétés'] ['score']) + = point_score sinon: poly ['properties'] ['score']) = point_score
klewis