大型光栅频率表/计数

Question

我尝试使用 freq() 计算 R 中栅格像素值的频率/计数。 创建两个示例栅格进行比较：

library(raster)
RastSmall <- raster(nrow=70, ncol=70) 
RastBig   <- raster(nrow=7000, ncol=7000) 
set.seed(0)
RastSmall[] <- round(runif(1:ncell(r_hr), 1, 5))
RastBig[] <- round(runif(1:ncell(r_hr), 1, 5))

使用freq()获取像素数

freq(RastSmall)
value    count
[1,]     1  6540000
[2,]     2 12150000
[3,]     3 12140000
[4,]     4 11720000
[5,]     5  6450000  

但是，它是一个相当大的文件，需要很长时间，即长达数小时。 R中有更快的方法吗？ 这里是小型和大型栅格的速度差异：

system.time(freq(RastSmall))
  user  system elapsed 
  0.008   0.000   0.004 
system.time(freq(RastBig))
  user  system elapsed 
  40.484   0.964  41.445 

有没有办法加快速度？或者，这可以在命令行中使用 gdal 工具之类的工具来完成吗？

Answer 1

上周我确实做到了，但是我找不到其他更快的方法在 R 中做到这一点。我尝试通过调用 GRASS 的 r.report 来使用 rqgis 包来做到这一点。它可以工作，但比 R 本机函数慢。也许你会有更好的运气。如果您想尝试一下，这是我的草代码：

library(RQGIS)
monqgis <- set_env("C:\\Mrnmicro\\Applic\\OSGeo4W")
find_algorithms(search_term = "report", qgis_env = monqgis)
get_usage(alg = "grass7:r.report", qgis_env = monqgis)
params <- get_args_man(alg = "grass7:r.report", qgis_env = monqgis)

get_usage(alg = "grass7:r.report", qgis_env = monqgis)

params$map <- classif
params$units <- 5
params$rawoutput <- "C:\\temp\\outputRQGIS_raw"
params$html <- "C:\\temp\\outputRQGIS"

system.time(asas <- run_qgis(alg = "grass7:r.report", params=params,load_output = params$OUTPUT, qgis_env = monqgis))

Answer 2

不是一个惊人的节省，但如果您从栅格中获取值，然后运行 ​​base::table 函数，它可以节省大约 20%。我的栅格是 c.500m 单元格。

# read in raster to obtain frequency table
r <- raster("./path/myraster.tif")

# perform tests; traditional freq() vs. getValues() & table()
require(microbenchmark)
  mbm <- microbenchmark(
    Freq = {freqf <- freq(r,useNA="no"); 
            freq.df <- data.frame(CODE=freqf[,1], N=freqf[,2]},
    GetVals = {v <- getValues(r);
               vt <- table(v); 
               getval.df <- data.frame(CODE=as.numeric(names(vt)),N=as.numeric(as.matrix(vt)))},
    times=5
  )
  mbm

Unit: seconds
    expr      min       lq     mean   median       uq      max neval
    Freq 191.1649 191.8001 198.8567 192.5256 193.0986 225.6942     5
 GetVals 153.5552 154.8776 156.9173 157.0539 159.0400 160.0598     5

# check the routines have identical results
identical(freq.df,getval.df)
[1] TRUE

我想有点节省 （注意，我制作数据框的原因是我继续处理来自频率分析的数据）

Answer 3

我认为最有效的计算方法是使用来自 GDAL 的GetHistogram( )。不幸的是，我找不到从 R 中使用它的方法。最接近的方法是使用 R 中的 gdalUtilities::gdalinfo，并使用标志 -hist 或 hist = TRUE，但计算范围限制在 0 - 255 之间。 另一个选项是使用rasterDT::freqDT，它比常规选项更快。举个例子：

library(gdalUtilities)
library(raster)
library(rasterDT)
library(microbenchmark)

RastBig   <- raster(nrow=7000, ncol=7000) 
set.seed(0)
RastBig[] <- round(runif(1:ncell(RastBig), 1, 5))
writeRaster(RastBig, filename = 'C:/temp/RastBig.tif')


mbm <- microbenchmark(times = 50,
  freq1 = freq(RastBig),
  freq2 = table(RastBig[]),
  freq3 = freqDT(RastBig),
  freq4 = ({
    
    gdalLog <- capture.output(gdalUtilities::gdalinfo(datasetname = 'C:/temp/RastBig.tif', hist = TRUE));
    (bucxml <- as.numeric(sub('buckets.+', '', grep('buckets ', gdalLog, value = TRUE))));
    (minxml <- as.numeric(gsub('.+from | to.+', '', grep('buckets ', gdalLog, value = TRUE)) ));
    (maxxml <- as.numeric(gsub('.+to |:', '', grep('buckets ', gdalLog, value = TRUE))));
    (histxml <- as.numeric(strsplit(split = '[[:space:]]', gsub("^ |^  ", "", gdalLog[grep('buckets', gdalLog)+1]))[[1]]));
    
    labs <- seq(from = minxml, to = maxxml, length.out = bucxml);
    df <- data.frame(labs, nwlab = c(ceiling(labs[1]),
                                      round(labs[2:(bucxml-1)]),
                                      floor(labs[bucxml])), 
                      val = histxml);
    hist <- aggregate(df$val, by = list(df$nwlab), sum)})
)

结果：

> freq1
 value    count
[1,]     1  6127755
[2,]     2 12251324
[3,]     3 12249376
[4,]     4 12248938
[5,]     5  6122607

> freq2

       1        2        3        4        5 
 6127755 12251324 12249376 12248938  6122607 

> freq3
   ID     freq
1:  1  6127755
2:  2 12251324
3:  3 12249376
4:  4 12248938
5:  5  6122607

> freq4
  Group.1        x
1       1  6127755
2       2 12251324
3       3 12249376
4       4 12248938
5       5  6122607


Unit: milliseconds
  expr          min           lq         mean       median        uq        max neval
 freq1 58628.486301 59100.539302 59400.304887 59383.913701 59650.412 60841.3975    50
 freq2 55912.170401 56663.025202 56954.032395 56919.905051 57202.001 58307.9500    50
 freq3  3785.767301  4006.858102  4288.699531  4292.447250  4536.382  4996.0598    50
 freq4     7.892201     8.883102     9.255641     9.154001     9.483    15.6072    50

编辑：使用它比选项 3 快得多：

rB <- raster('C:/temp/RastBig.tif')
freq3B <- freqDT(rB)