Eine schnellere Funktion die Auflösung eines Rasters R zu senken

Question

ich das Raster Paket bin mit der Auflösung der großen Rastern zu senken, indem die Funktion Aggregat wie diese

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) 

die Rasterbilder Ich habe zu aggregieren viel sind größer 34000x34000, also würde ich gerne wissen, ob es einen schnelleren Weg gibt, die Funktion agg.fun zu implementieren.

Answer 1

Sie können hierfür gdalUtils::gdalwarp verwenden. Für mich ist es weniger effizient als @ JosephWoods fasterAgg.Fun für Raster mit 1.000.000 Zellen, aber für Josephs größeres Beispiel ist es viel schneller. Voraussetzung dafür ist, dass das Raster auf der Festplatte vorhanden ist. Schreiben Sie die Zeit in den unteren Bereich, wenn sich Ihr Raster im Speicher befindet.

Unten habe ich die Änderung von fasterAgg.Fun verwendet, die den häufigsten Wert zurückgibt, anstatt seinen Index im Block.

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 

Man beachte, dass es ein leichter Unterschied Bindungen in der Definition des Modus, wenn es: gdalwarp den höchsten Wert auswählt, während die zu aggregate oben geleitet Funktionen (über which.max ‚Verhalten) Wählen Sie die niedrigste (zB , siehe which.max(table(c(1, 1, 2, 2, 3, 4)))).

Auch das Speichern der Rasterdaten als Ganzzahl ist wichtig (falls zutreffend). Wenn die Daten als Float gespeichert sind (z. B. writeRaster), dauert die obige Operation gdalwarp ~ 14 Sekunden auf meinem System. Verfügbare Typen finden Sie unter ?dataType.

Answer 2

Try this:

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 

Für eine größere Testobjekt, die wir haben:

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 

Wenn Sie die aktuellen Werte wollen als @digEmAll oben in den Kommentaren sagt, einfach den Wert Rückkehr ändern in myRle.Alt von which.max(diff(c(0L, i))) bis x1[i][which.max(diff(c(0L, i)))].

Answer 3

Just for fun habe ich auch eine RCPP Funktion (nicht viel schneller als @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

Wenn Ihr Ziel Aggregation ist, würden Sie nicht wollen, dass die max Funktion?

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