Suite

Statistiques zonales pondérées en Python et/ou R


J'ai des fichiers raster avec des données de précipitations (4x4 km), mensuellement, depuis environ 100 ans. J'ai un fichier de formes pour un bassin versant (qui comprend plusieurs polygones en tant que sous-bassins). J'ai besoin des précipitations moyennes pour le bassin versant pour chaque mois. Je ne veux pas le faire à la main dans ArcMap car cela devra probablement être répété plusieurs fois (pour d'autres bassins versants et lorsque davantage de données seront disponibles), donc je le code.

J'ai maintenant écrit du code en R, en utilisant des packagesrasteretrgdal. Cela fonctionne très bien dans le sens où il a une fonctionnalité pour calculer la moyenne pondérée par sous-bassin ; il vérifie quelle fraction d'une cellule raster appartient au polygone et l'utilise pour calculer une moyenne pondérée pour le polygone. Étant donné que certains des sous-bassins sont petits par rapport au raster, il s'agit d'une fonctionnalité intéressante. Et en plus, parce que je peux accéder aux informations du fichier de formes telles que la superficie des sous-bassins, je peux calculer une moyenne pondérée pour l'ensemble du bassin versant.

Cependant - c'est lent. Il le dit aussi dans la description de la fonctionextrait(). Le calcul de la valeur pondérée prend 4 à 5 secondes, donc lorsque j'ai exécuté ce script pendant 30 ans de données mensuelles, cela a pris 37 minutes (il y a quelques calculs supplémentaires en plus de la moyenne pondérée).

J'ai le sentiment qu'il pourrait être plus rapide d'utiliser Python, mais je ne sais pas comment aborder cela. J'ai commencé avec cela avant le code R (voir ci-dessous) mais j'ai déplacé mon attention vers R. Cela me semblait un peu intimidant car cela faisait un moment que je code en Python, et je n'ai jamais utilisé Python pour des analyses géospatiales.

De plus, un test initial semble montrer qu'en Python, il n'y a pas de moyenne pondérée et que les petits polygones qui ne contiennent pas le centre d'une cellule raster sont attribués0ou alorsRien.

Alors mes questions sont :

  1. Si c'est possible en Python, la moyenne pondérée de cette façon sera-t-elle plus rapide ?
  2. Comment accéder aux autres données associées au fichier de formes (telles que la zone) ?

Vous trouverez ci-dessous le code très très basique que j'ai jusqu'à présent.

from rasterstats import zonal_stats import gdal bassin="Steinhatchee_Project.shp" rain="PRISM_ppt_stable_4kmM3_201504_bil.bil" # REMARQUE De la documentation Python : "Il n'y a pas de bonne ou de mauvaise façon de rastériser un vecteur. # La stratégie par défaut est d'inclure tous les pixels le long du chemin de rendu de ligne (pour les lignes), # ou les cellules où le point central est dans le polygone (pour les polygones) stats=zonal_stats(basin,rain) # obtenir la moyenne pour chaque polygone (il y en a 19) data=([f['mean '] pour f dans les statistiques]) # supprimez les valeurs "Aucun"… data2=[x pour x dans les données si x n'est pas Aucun] #… afin que je puisse calculer la moyenne pour l'ensemble du bassin versant (données2)

Les donnéesressemble à ce qui suit. Noter laRien

Out [76]: [115,64999389648438, 83,07333374023438, 92,33000183105469, 73,35416666666667, 97,36499786376953, 92,03500366210938, 80,30999755859375, 102,8699951171875, 94,3499984741211, 110,99250030517578, 80,18000030517578, 67,6749979654948, 107,825, 82,63999938964844, 88,375, 59,769999186197914, 63,48249816894531, None, None]

Mes spécifications système : Windows 7, Python 2.7.9 (en utilisant l'IDE Enthought Canopy) 32 bits, GDAL-1.11.3, rasterio-0.29.0, Shapely-1.5.13


« Si c'est possible en Python, est-ce que la moyenne pondérée de cette façon sera plus rapide ? »

Voir la réponse à la question suivante. La collecte des statistiques agrégées ainsi que la rastérisation sous le capot est ce qui prend du temps. Vous pouvez obtenir votre moyenne pondérée en utilisant la réponse suivante. En ce qui concerne ce qui est plus rapide, je ne peux pas le commenter… si vous êtes prêt à recoder les rasterstats, cela peut être rendu un peu plus rapide mais vous perdriez son utilité.

« Comment accéder aux autres données associées au fichier de formes (telles que la zone) ?"

Pour les données renvoyées par statistiques_zonales, il y a une valeur "count" qui, je crois, est le nombre de cellules rastérisées qui contribuent. Si vous connaissez la résolution de votre raster (et donc la superficie par pixel), vous avez la superficie (area = area_per_pixel * 'count').