Suite

Modification des valeurs de calcul raster - raccourcir ce script ?

Modification des valeurs de calcul raster - raccourcir ce script ?


je fais quelque chose d'assez simple concernant l'algèbre raster mais j'ai du mal à trouver la bonne fonction ou à faire fonctionner une fonction correctement;

Fondamentalement, j'ai 2 rasters, représentant des années consécutives, qui représentent une classification (y compris les valeurs NA);

r <- raster(ncol=10,nrow=10) r[] <- sample(c(1,2,4,8),size=100,replace=T) r[runif(10*10) >= 0.50 ] <- NA r1 <- raster(ncol=10,nrow=10) r1[] <- sample(c(1,2,4,8),size=100,replace=T) r1[runif(10*10 ) >= 0,50] <- NA

pour obtenir mon changement d'une année (r) à l'autre (r1), je soustrais simplement r1 de r ;

r2 = r - r1 fréq(r2,chiffres=2)

Chaque nouvel entier code pour un changement spécifique, 0 signifie simplement aucun changement. Cela ne me dérange pas que les cellules deviennent NA si NA n'existe que dans une couche raster, c'est bien. Ce que je veux faire, c'est examiner de plus près les valeurs 0, donc si la valeur d'origine dans r est 1 et que la valeur récente dans r1 est également 1, je veux le savoir dans le r2 résultant - c'est-à-dire quantifier et analyser le " non changer les cellules. Idem pour r=2 & r1=2, r=4 & r1=4 etc - pas seulement tous comme des 0 mais comme une collection de nouveaux codes inutilisés pour représenter quelles cellules sont restées les mêmes (et comment elles sont restées les mêmes) d'un année à l'autre.

Fondamentalement, j'ai extrait la position des cellules de r2 qui sont égales à 0 en tant que trame de données de points spatiaux;

pts <-rasterToPoints(r2,fun=function(x){x==0},spatial=T) plot(pts)

puis j'ai remplacé les cellules de r2 qui sont égales à 0 par des valeurs de r (ou r1, peu importe car elles ont la même valeur) en utilisant les emplacements dérivés ci-dessus ;

r2[r2==0] <- (extract(r,pts))/10 plot(r2) freq(r2,digits=2)

diviser les valeurs de remplacement par 10 garantit qu'elles ne tombent pas dans la même « corbeille » que tout autre changement de code mappé.

Je suis sûr qu'il y a un moyen plus rapide? je pensais pouvoir créer une pile raster à partir de r et r1, créer un r2, puis utiliser une sorte de fonction "où" ou "ifelse" pour effectuer les calculs ci-dessus, mais tout ce que je fais entraîne des erreurs ou une négligence des valeurs NA pour faire fonctionner les fonctions. Ou des erreurs S4/entier.

Je ne suis toujours pas sûr que cela fonctionnera sur un grand échantillon pour le moment, et mes ensembles de données finaux sont très volumineux.


Ce que vous voulez, c'est un calcul conditionnel : retourner la valeur dern'importe quandretr1sont égaux et sinon définissez la sortie sur NA.

Les opérations arithmétiques cellule par cellule semblent être les plus rapides. (Ils sont beaucoup plus rapide que, disons, en utilisantmasquerou les fonctions de reclassification.) Puisqu'elles ne semblent pas offrir un opérateur conditionnel réel, utilisez deux astuces séculaires :

  1. Traitez les logiques comme des nombres.FAUXest 0 etVRAIest 1 dans les opérations arithmétiques.

  2. Créez des valeurs NA (ou, presque aussi efficacement, des valeurs infinies) en utilisant des opérations arithmétiques invalides.

Une solution est

r3 <- r == r1 r3 <- r3 * r * (1/r3)

Cela fonctionne parce que quandretr1sont égaux, les deuxr3et1/r3égale 1 et les multiplications ne changent rien : elles renvoient la valeur der. Lorsqueretr1ne sont pas égaux,1/r3est indéfini, produisant un résultat infini. Par conséquent,fréqne totalise que les valeurs des cellules oùretr1se mettre d'accord.

Sur ma machine, ce calcul prend environ une seconde pour les rasters avec 10 000 000 de cellules. (C'est environ dix fois plus long que la simple comparaisonr - r1.) Il évoluera proportionnellement au nombre de cellules jusqu'à ce que la pagination du disque soit invoquée, auquel cas vous serez à la merci de votre débit de stockage.

(Si vous pouvez mettre toutes les données dans la RAM, c'est encore plus rapide à utiliserRdes opérations intégrées sur le tableau de données, puis le reconvertir en objet raster.)


Depuis que whuber m'a donné le bon morceau de code rapide pour créer un raster avec les mêmes valeurs entre les rasters, j'ai pensé que je finirais tout le travail;

créer des rasters et soustraire l'un de l'autre pour obtenir un raster de « modification » (plein de 0 qui doivent également être examinés);

r <- raster(ncol=10,nrow=10) r[] <- sample(c(1,2,4,8),size=100,replace=T) r[runif(10*10) >= 0.50 ] <- NA r1 <- raster(ncol=10,nrow=10) r1[] <- sample(c(1,2,4,8),size=100,replace=T) r1[runif(10*10 ) >= 0,50] <- NA r2 = r - r1

puis, en utilisant le petit bout de code de whuber, créez un nouveau raster où les valeurs sont les mêmes (en divisant par 10 pour que les valeurs des cellules ne correspondent à rien de calculé précédemment);

r0 <- r == r1 r0 <- (r0 * r * (1/r0))/10

Pour combiner le raster de changement d'origine avec le nouveau raster qui analyse les cellules qui n'ont pas de changement, ajoutez-les ensemble, en changeant d'abord les valeurs NA en 0 ;

r0[est.na(r0)] <- 0 rfinal <- r2 + r0

qui crée un nouveau raster avec toutes les modifications et toutes les valeurs 0 modifiées en un code spécifique pour examiner exactement quel type de « pas de changement » se produit.

exporter au format csv ;

counts <- freq(rfinal,digits=2,useNA='no') write.csv(counts,"counts.csv")