Suite

Accélérer le temps d'extraction de la proportion de types de couverture terrestre dans la zone tampon à partir du raster à l'aide de R?


Je voudrais extraire des données d'occupation du sol dans une zone tampon de 10 km autour des objets de la classe SpatialLines et calculer la proportion de chaque type d'occupation du sol. J'ai donc utilisé la fonctionextrait(paquetraster). j'ai aussi utilisé la fonctionrecadrerpour recadrer mon raster. Dans mon raster, j'ai 10 types d'occupation du sol. Voici mon code :

buf_line <- gBuffer(seg_line, width=10000) # seg_line = Lines objects ha <-extract(x=data_raster,y=buf_line) # La proportion de chaque type d'occupation du sol doit être en colonnes (une colonne = un type d'occupation du sol) ha_1 <-length(ha[[1]][ha[[1]]==1])/length(ha[[1]]) ha_2 <-length(ha[[1]][ha[[1] ]==2])/longueur(ha[[1]]) ha_3 <-longueur(ha[[1]][ha[[1]]==3])/longueur(ha[[1]]) ha_4 <-length(ha[[1]][ha[[1]]==4])/length(ha[[1]]) ha_5 <-length(ha[[1]][ha[[1]] ==5])/longueur(ha[[1]]) ha_6 <-longueur(ha[[1]][ha[[1]]==6])/longueur(ha[[1]]) ha_7 < -length(ha[[1]][ha[[1]]==7])/length(ha[[1]]) ha_8 <-length(ha[[1]][ha[[1]]= =8])/longueur(ha[[1]]) ha_9 <-longueur(ha[[1]][ha[[1]]==9])/longueur(ha[[1]]) ha_10 <- longueur(ha[[1]][ha[[1]]==10])/longueur(ha[[1]])

j'ai utilisé ce code en boucles'appliquerpour extraire des données d'occupation du sol dans une zone tampon de 10 km autour de 30 000 lignes. Le problème est que l'extraction autour de 30 000 lignes est très lente.

Comment accélérer le temps de traitement ?


Je trouve que, pour des problèmes comme celui-ci, on peut accélérer un peu les choses si vous recadrez le raster dans le tampon, forcez le vecteur/matrice, puis effectuez des calculs.

Jetez également un œil à table et prop.table pour calculer vos proportions d'occupation du sol. Voici un exemple de polygone, ce que vous recherchez étant donné le tampon de ligne. Le résultat ici sera une liste "prop" qui est ordonnée de la même manière que les données de polygone source. Vous pouvez utiliser lapply ou do.call sur la liste pour résumer ou réduire les résultats dans un data.frame où les colonnes sont les classes.

Tout d'abord, créez des exemples de données et tracez-les

bibliothèque(raster) bibliothèque(sp) p <- raster(nrow=10, ncol=10) p[] <- runif(ncell(p)) * 10 p <- rasterToPolygons(p, fun=function(x){x > 9}) r <- raster(nrow=100, ncol=100) r[] <- round(runif(ncell(r), 1,5),0) plot(r) plot(p, add=TRUE, lwd=4)

Maintenant, nous créons une liste vide "prop" pour stocker les résultats et écrivons une boucle for qui parcourt les polygones, recadre le raster à chaque polygone de sous-ensemble, utilise prop.table pour obtenir les proportions des valeurs nominales (par exemple, classe de couverture terrestre) et écrit les résultats dans l'objet liste.

prop <- list() for(i in 1:nrow(p)){ lsub <- p[i, ] cr <- raster::crop(r, raster::extent(lsub), snap = "out") fr <- raster::rasterize(lsub, cr) r.sub <- raster::mask(x = cr, mask = fr) prop[[i]] <- prop.table(table(getValues(r.sub) )) } as.data.frame( do.call("rbind", prop) )

Ajouter que vous auriez besoin de modifier pour votre analyse serait de passer le code d'un objet de ligne et d'utiliser gBuffer pour créer le polygone pour chaque ligne.

Si vous disposez de la RAM disponible pour traiter les rasters, une alternative serait d'utiliser la classe sp "SpatialPixlesDataFrame" pour les objets raster et la fonction over. Vous pouvez importer un raster, dans un raster de classe sp, à partir du disque en utilisant readGDAL dans le package rgdal. En commençant par notre exemple précédent, voici à quoi ressemblerait le code.

r <- as(r, "SpatialPixelsDataFrame") # forcer r dans un objet sp prop <- list() for(i in 1:nrow(p)){ lsub <- p[i, ] r.sub <- r [!is.na(sp::over(r, sp::geometry(lsub))), ] prop[[i]] <- prop.table(table([email protected][,1])) }