Una funzione più veloce per ridurre la risoluzione di un raster R

Question

Sto usando il pacchetto raster di abbassare la risoluzione delle grandi raster, utilizzando la funzione di aggregazione come questo

require(raster)  
x <- matrix(rpois(1000000, 2),1000) 

a <-raster(x) 
plot(a) 

agg.fun <- function(x,...) 
    if(sum(x)==0){ 
     return(NA) 
    } else { 
     which.max(table(x)) 
    } 

a1<-aggregate(a,fact=10,fun=agg.fun) 
plot(a1) 

le immagini raster devo aggregare sono molto più grande 34000x34000 quindi mi piacerebbe sapere se c'è un modo più veloce per implementare la funzione agg.fun.

Answer 1

È possibile utilizzare gdalUtils::gdalwarp per questo. Per me, è meno efficiente di @fasterAgg.Fun di JosephWood per i raster con 1.000.000 di celle, ma per l'esempio più grande di Joseph è molto più veloce. Richiede che il raster esista su disco, quindi calcoli il tempo di scrittura nel sotto se il tuo raster è in memoria.

Di seguito, ho utilizzato la modifica di fasterAgg.Fun che restituisce il valore più frequente di anziché il suo indice nel blocco.

library(raster) 
x <- matrix(rpois(10^8, 2), 10000) 
a <- raster(x) 

fasterAgg.Fun <- function(x,...) { 
    myRle.Alt <- function (x1) { 
     n1 <- length(x1) 
     y1 <- x1[-1L] != x1[-n1] 
     i <- c(which(y1), n1) 
     x1[i][which.max(diff(c(0L, i)))] 
    } 

    if (sum(x)==0) { 
     return(NA) 
    } else { 
     myRle.Alt(sort(x, method="quick")) 
    } 
} 
system.time(a2 <- aggregate(a, fact=10, fun=fasterAgg.Fun)) 

## user system elapsed 
## 67.42 8.82 76.38 

library(gdalUtils) 
writeRaster(a, f <- tempfile(fileext='.tif'), datatype='INT1U') 
system.time(a3 <- gdalwarp(f, f2 <- tempfile(fileext='.tif'), r='mode', 
          multi=TRUE, tr=res(a)*10, output_Raster=TRUE)) 

## user system elapsed 
## 0.00 0.00 2.93 

noti che v'è una leggera differenza nella definizione della modalità quando ci sono vincoli: gdalwarp seleziona il valore più alto, mentre le funzioni passate al aggregate sopra (attraverso comportamenti which.max s') selezionare il più basso (ad esempio , vedi which.max(table(c(1, 1, 2, 2, 3, 4)))).

Inoltre, la memorizzazione dei dati raster come numero intero è importante (se applicabile). Se i dati sono memorizzati come float (il valore predefinito di writeRaster), ad esempio, l'operazione gdalwarp in alto richiede circa 14 secondi sul mio sistema. Vedi ?dataType per i tipi disponibili.

Answer 2

Prova questo:

fasterAgg.Fun <- function(x,...) { 
    myRle.Alt <- function (x1) { 
     n1 <- length(x1) 
     y1 <- x1[-1L] != x1[-n1] 
     i <- c(which(y1), n1) 
     which.max(diff(c(0L, i))) 
    } 

    if (sum(x)==0) { 
     return(NA) 
    } else { 
     myRle.Alt(sort(x, method="quick")) 
    } 
} 

library(rbenchmark) 
benchmark(FasterAgg=aggregate(a, fact=10, fun=fasterAgg.Fun), 
      AggFun=aggregate(a, fact=10, fun=agg.fun), 
      replications=10, 
      columns = c("test", "replications", "elapsed", "relative"), 
      order = "relative") 
     test replications elapsed relative 
1 FasterAgg   10 12.896 1.000 
2 AggFun   10 30.454 2.362 

Per un oggetto di prova più grande, abbiamo:

x <- matrix(rpois(10^8,2),10000) 
a <- raster(x) 
system.time(a2 <- aggregate(a, fact=10, fun=fasterAgg.Fun)) 
    user system elapsed 
111.271 22.225 133.943 

system.time(a1 <- aggregate(a, fact=10, fun=agg.fun)) 
    user system elapsed 
282.170 24.327 308.112 

Se si desidera che i valori attuali come dice @digEmAll nei commenti di cui sopra, è sufficiente modificare il valore di ritorno in myRle.Alt da which.max(diff(c(0L, i))) a x1[i][which.max(diff(c(0L, i)))].

Answer 3

Solo per divertimento ho creato anche una funzione Rcpp (non molto più veloce di @JosephWood):

########### original function 
#(modified to return most frequent value instead of index) 
agg.fun <- function(x,...){ 
    if(sum(x)==0){ 
     return(NA) 
    } else { 
     as.integer(names(which.max(table(x)))) 
    } 
} 
########### @JosephWood function 
fasterAgg.Fun <- function(x,...) { 
    myRle.Alt <- function (x1) { 
     n1 <- length(x1) 
     y1 <- x1[-1L] != x1[-n1] 
     i <- c(which(y1), n1) 
     x1[i][which.max(diff(c(0L, i)))] 
    } 

    if (sum(x)==0) { 
     return(NA) 
    } else { 
     myRle.Alt(sort(x, method="quick")) 
    } 
} 
########### Rcpp function 
library(Rcpp) 
library(inline) 

aggrRcpp <- cxxfunction(signature(values='integer'), ' 
       Rcpp::IntegerVector v(clone(values)); 
       std::sort(v.begin(),v.end()); 
       int n = v.size(); 
       double sum = 0; 
       int currentValue = 0, currentCount = 0, maxValue = 0, maxCount = 0; 
       for(int i=0; i < n; i++) { 
        int value = v[i]; 
        sum += value; 
        if(i==0 || currentValue != value){ 
         if(currentCount > maxCount){ 
         maxCount = currentCount; 
         maxValue = currentValue; 
         } 
         currentValue = value; 
         currentCount = 0; 
        }else{ 
         currentCount++; 
        } 
       } 
       if(sum == 0){ 
        return Rcpp::IntegerVector::create(NA_INTEGER); 
       } 
       if(currentCount > maxCount){ 
        maxCount = currentCount; 
        maxValue = currentValue; 
       } 
       return wrap(maxValue) ; 
     ', plugin="Rcpp", verbose=FALSE, 
     includes='') 
# wrap it to support "..." argument 
aggrRcppW <- function(x,...)aggrRcpp(x); 

Benchmark:

require(raster) 
set.seed(123) 
x <- matrix(rpois(10^8, 2), 10000) 
a <- raster(x) 

system.time(a1<-aggregate(a,fact=100,fun=agg.fun)) 
# user system elapsed 
# 35.13 0.44 35.87 
system.time(a2<-aggregate(a,fact=100,fun=fasterAgg.Fun)) 
# user system elapsed 
# 8.20 0.34 8.59 
system.time(a3<-aggregate(a,fact=100,fun=aggrRcppW)) 
# user system elapsed 
# 5.77 0.39 6.22 

########### all equal ? 
all(TRUE,all.equal(a1,a2),all.equal(a2,a3)) 
# > [1] TRUE 

Answer 4

Se l'obiettivo è l'aggregazione, potrebbe non si desidera che il max funzione?

library(raster)  
x <- matrix(rpois(1000000, 2),1000) 
a <- aggregate(a,fact=10,fun=max)