Obtenez des valeurs raster à partir d'une superposition de polygones dans les solutions SIG Open Source

J'ai deux couches. Une couche en forme de polygone avec de nombreuses tuiles et une couche raster contenant la couverture terrestre CORINE 2006 avec de nombreuses catégories dans une carte des couleurs. Je veux obtenir pour chaque polygone de la couche de formes une somme de chaque catégorie de couverture terrestre de la couche raster.

Par exemple, il y a un polygone avec l'ID '2' et je veux des attributs comme celui-ci pour ce polygone (en pourcentage ou en mètres carrés):

• Terre arable: 15%
• Forêt: 11%
• Rues: 2% (... et donc une)

J'ai essayé de le faire dans l'herbe, qgis (pas de fonction), saga (résume juste chaque à une valeur totale) r (somme totale), mais je n'ai toujours trouvé aucune solution. La plupart des plugins (statistiques zonales dans qgis) ne prennent en charge que les couches raster 0-1. v.rast.stats n'a pas aidé non plus. Je suis ouvert à toute solution bonne et intelligente!. Peut-être que j'ai même utilisé une mauvaise approche ou fait des erreurs.

Dans Arcgis, cette tâche est assez facile, si je me souviens bien, mais il me manque toujours une bonne solution pour votre utilisateur Linux quotidien.

J'utilise un système Linux Debian et c'est pourquoi je ne peux utiliser que des programmes pour ce système d'exploitation.

EDIT: Parce que cette question a encore tant de vues et de visiteurs: j'ai écrit un plugin QGIS, qui est également capable de calculer la couverture terrestre de la couche raster. Je n'ai pas encore codé de superposition de polygones, mais c'est définitivement planifié. Trouvez le plugin ici et installez d'abord la bibliothèque Scipy.

Utilisez «extraire» pour superposer les entités surfaciques d'un SpatialPolygonsDataFrame (qui peut être lu à partir d'un fichier de formes à l'aide de maptools: readShapeSpatial) sur un raster, puis utilisez «table» pour résumer. Exemple:

> c=raster("cumbria.tif") # this is my CORINE land use raster
> summary(spd)
Object of class SpatialPolygonsDataFrame
[...]
> nrow(spd)  # how many polygons:
[1] 2
> ovR = extract(c,spd)
> str(ovR)
List of 2
 $ : num [1:542] 26 26 26 26 26 26 26 26 26 26 ...
 $ : num [1:958] 18 18 18 18 18 18 18 18 18 18 ...

Donc mon premier polygone couvre 542 pixels, et mon deuxième polygone 958. Je peux résumer chacun d'eux:

> lapply(ovR,table)
[[1]]

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792 

Mon premier polygone est donc de 287 pixels de classe 26 et 255 pixels de classe 27. Assez facile à additionner, à diviser et à multiplier par 100 pour obtenir des pourcentages.

Je voulais faire rapport et me voici. La solution de Spacedman a très bien fonctionné et j'ai pu exporter toutes les informations pour chaque polygone de ma forme. Juste au cas où quelqu'un aurait le même problème, voici comment j'ai précédé:

...
tab <- apply(ovR,table)
# Calculate percentage of landcover types for each polygon-field.
# landcover is a datastream with the names of every polygon
for(i in 1:length(tab)){
 s <- sum(tab[[i]])
 mat <- as.matrix(tab[[i]])
 landcover[i,paste("X",row.names(mat),sep="")] <- as.numeric(tab[[i]]/s)
}

si je comprends bien ce que vous voulez, et en supposant que vous avez déjà la couche vectorielle 'mypolygonlayer' et la couche raster 'corina' dans votre base de données GRASS GIS:

Je convertirais d'abord le vecteur en raster. Le chat s'assurera que vous aurez un identifiant unique par polygone. Si vous avez une colonne avec un identifiant numérique unique, vous pouvez utiliser cette colonne à la place. La colonne d'étiquette est facultative:

v.to.rast input = couche mypolygonlayer = 1 output = mypolygons use = cat labelcolumn = NameMappingUnit

Exécutez ensuite r.stats pour obtenir vos statistiques:

r.stats -a -l input = mypolygons, corina separator =; sortie = / home / paulo / corinastats.csv

La dernière étape consiste à ouvrir le corinastats.csv dans par exemple LibreOffice et créer un tableau croisé dynamique ou utiliser R pour créer votre tableau croisé

Je sais que ce message est assez ancien, mais j'ai récemment refusé d'effectuer le même type d'analyse, mais le téléchargement de programmes tels que R est un peu compliqué sur mon ordinateur de travail et doit être approuvé. Après plusieurs heures de recherche sur une méthode que je ne pouvais utiliser qu'avec QGis et Excel, j'ai trouvé que cette méthode fonctionnait le mieux pour moi et je voulais l'offrir aux personnes dans le même genre de situation.

1. Couper le polygone sur la couche raster (Raster → extraction → clipper: fichier d'entrée = couche raster, choisissez votre nom et emplacement de sortie, cliquez sur la couche de masque, choisissez votre polygone → ok)

2. Polygoniser le calque de tondeuse (Raster → Conversion → polygoniser: fichier d'entrée = votre calque de clip, enregistrer la sortie → ok)

3. Calcul du nombre de pixels (Cliquez sur le fichier de forme que vous venez de créer → calculateur de champ ouvert: cochez «créer un nouveau champ» et ajoutez le nom du champ, Fonction = géométrie → zone → ok). Vous devriez maintenant avoir une nouvelle colonne dans votre table attributaire indiquant le nombre de pixels.

4. Enregistrer la couche de polygonisation (clic droit sur la couche de polygonisation, enregistrer sous: format = fichier DBF, enregistrer sous → ok)

5. Résumant le nombre de pixels pour chaque habitat (démarrez Excel, ouvrez le fichier, si vous n'avez pas de titres, ajoutez-en un pour chaque colonne, cliquez sur une cellule vide, accédez à l'onglet DONNÉES, consolidez, assurez-vous que c'est sur la somme, cliquez sur le flèche rouge à côté de "parcourir ..." et sélectionnez deux colonnes (titres inclus), cliquez sur "ajouter" et cochez les cases "Ligne du haut" et "colonne de gauche" → ok)

6. Si, comme moi, vous avez beaucoup de polygones à analyser et devez les comparer dans le même tableau, cette étape vous sera utile. Dans un nouveau classeur Excel, répertoriez vos numéros d'habitats dans la colonne A (pour moi 1 à 48) et placez les deux colonnes que vous venez de consolider dans les colonnes B et C (habitat en B et nombre de pixels en C). Dans la cellule D1, écrivez la formule suivante: = IFNA (INDEX (C: C; MATCH (A2; B: B; 0)); "") et faites glisser (ou double-cliquez dans le coin inférieur droit) jusqu'à votre dernière valeur (donc si vous avez 48 habitats jusqu'à la cellule D48). Le nombre de pixels est maintenant dans les cellules correspondantes de votre habitat et vous pouvez répéter ce processus pour tous vos polygones.

Que diriez-vous de convertir les données CORINE en un jeu de données de polygones vectoriels en utilisant QGIS ( Raster> Conversion> Polygonize ) puis en utilisant la fonction Union ( Vector> Outils de géotraitement> Union ) pour combiner avec les polygones. L'ensemble de données vectorielles résultant contiendrait les zones de chaque classe CORINE dans chaque polygone.

QGIS.

Dans le tronc QGIS, il existe une autre version de ZonalStats disponible, elle s'appelle Zonal Statistics.

Ceci exécute la fonction dont vous avez besoin.

En ce qui concerne le flux de travail, je ne sais pas combien de rasters vous avez ou sont-ils simplement des bandes dans un raster?

Contrairement à la plupart des réponses ci-dessus, je dirais que la meilleure option est de pixelliser vos polygones et de travailler avec deux jeux de données raster au lieu de deux jeux de données polygonaux. Ceci demande beaucoup moins de traitement et est par conséquent la seule solution facile à implémenter pour traiter de grands rasters et de gros fichiers de polygones dans R.

Après avoir pixellisé vos polygones avec exactement la même étendue et résolution des données raster, vous pouvez tabuler des statistiques récapitulatives comme expliqué ici , ce qui est approprié si votre raster tient en mémoire (couches raster petites / moyennes) ou vous pouvez binariser chaque catégorie avec la reclassfonction et que de calculer des zonalstatistiques pour chaque classe. Voici une solution qui incorpore la pixellisation et les statistiques zonales dans une fonction et fonctionne bien avec de très grands ensembles de données.