gdal_calc amin 在传递超过 23 个输入文件时失败

Question

我编写了一个调用 gdal_calc.py 的 R 函数来计算 RasterStack（一系列输入栅格文件）中的像素最小值。我这样做是因为它比大型栅格的 raster::min 快得多。该函数适用于最多 23 个文件，但在传递 24 个或更多文件时会引发警告，并返回一个用零填充的输出栅格。

由于 R 只是准备对 python gdal_calc.py 的系统调用，所以这个问题并不特定于 R，我鼓励 python/numpy 爱好者继续阅读。

这是函数。最终 gdal_calc 调用的结构显示在此 post.

底部有问题的用法引发的警告消息中

gdal_min <- function(infile, outfile, return_raster=FALSE, overwrite=FALSE) {
  require(rgdal)
  if(return_raster) require(raster)
  # infile:        The multiband raster file (or a vector of paths to multiple 
  #                 raster files) for which to calculate cell minima.
  # outfile:       Path to raster output file.
  # return_raster: (logical) Should the output raster be read back into R?
  # overwrite:     (logical) Should outfile be overwritten if it exists?
  gdal_calc <- Sys.which('gdal_calc.py')
  if(gdal_calc=='') stop('gdal_calc.py not found on system.')
  if(file.exists(outfile) & !overwrite) 
    stop("'outfile' already exists. Use overwrite=TRUE to overwrite.", 
         call.=FALSE)
  nbands <- sapply(infile, function(x) nrow(attr(GDALinfo(x), 'df')))
  if(length(infile) > 26 || nbands > 26) stop('Maximum number of inputs is 26.')
  if(length(nbands) > 1 & any(nbands > 1)) 
    warning('One or more rasters have multiple bands. First band used.')

  if(length(infile)==1) {
    inputs <- paste0('-', LETTERS[seq_len(nbands)], ' ', infile, ' --', 
                     LETTERS[seq_len(nbands)], '_band ', seq_len(nbands), collapse=' ')
    n <- nbands
  } else {
    inputs <- paste0('-', LETTERS[seq_along(nbands)], ' ', infile, ' --', 
                     LETTERS[seq_along(nbands)], '_band 1', collapse=' ')
    n <- length(infile)
  }

  message('Calculating minimum and writing to ', basename(outfile))
  system2('python',
          args=c(gdal_calc, inputs,
                 sprintf("--outfile=%s", outfile),
                 sprintf('--calc="amin([%s], axis=0)"', 
                         paste0(LETTERS[seq_len(n)], collapse=',')),
                 '--co="COMPRESS=LZW"',
                 if(overwrite) '--overwrite'),
          stdout=FALSE
  )
  if(return_raster) raster(outfile) else invisible(NULL)
}

这是一个虚拟栅格网格，写入当前工作目录中的 'test.tif'。为简单起见，并为了证明问题不是由特殊的单个栅格引起的，我将同一个输入文件多次传递给 gdal_calc，而不是传递多个不同的文件。

library(raster)
writeRaster(raster(matrix(runif(100), 10)), 'test.tif')

该函数适用于 23 个输入文件：

r <- gdal_min(rep('test.tif', 23), f_out <- tempfile(fileext='.tif'),
              return_raster=TRUE)
plot(r)

...但不适用于 ≥ 24:

r <- gdal_min(rep('test.tif', 24), f_out <- tempfile(fileext='.tif'),
              return_raster=TRUE)
## Calculating minimum and writing to file2310254bcc.tif
## Warning message:
## running command '"python" C:\Python33\Scripts\GDAL_C~1.PY -A test.tif --A_band 1 -B test.tif --B_band 1 -C test.tif --C_band 1 -D test.tif --D_band 1 -E test.tif --E_band 1 -F test.tif --F_band 1 -G test.tif --G_band 1 -H test.tif --H_band 1 -I test.tif --I_band 1 -J test.tif --J_band 1 -K test.tif --K_band 1 -L test.tif --L_band 1 -M test.tif --M_band 1 -N test.tif --N_band 1 -O test.tif --O_band 1 -P test.tif --P_band 1 -Q test.tif --Q_band 1 -R test.tif --R_band 1 -S test.tif --S_band 1 -T test.tif --T_band 1 -U test.tif --U_band 1 -V test.tif --V_band 1 -W test.tif --W_band 1 -X test.tif --X_band 1 --outfile=C:\Users\John\AppData\Local\Temp\RtmpK4heMM\file2310254bcc.tif --calc="amin([A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X], axis=0)" --co="COMPRESS=LZW"' had status 1

直接从终端执行后者时给出的回溯是：

"python" C:\Python33\Scripts\GDAL_C~1.PY -A test.tif --A_band 1 -B test.tif --B_band 1 -C test.tif --C_band 1 -D test.tif --D_band 1 -E test.tif --E_band 1 -F test.tif --F_band 1 -G test.tif --G_band 1 -H test.tif --H_band 1 -I test.tif --I_band 1 -J test.tif --J_band 1 -K test.tif --K_band 1 -L test.tif --L_band 1 -M test.tif --M_band 1 -N test.tif --N_band 1 -O test.tif --O_band 1 -P test.tif --P_band 1 -Q test.tif --Q_band 1 -R test.tif --R_band 1 -S test.tif --S_band 1 -T test.tif --T_band 1 -U test.tif --U_band 1 -V test.tif --V_band 1 -W test.tif --W_band 1 -X test.tif --X_band 1 --outfile=C:\Users\John\AppData\Local\Temp\RtmpK4heMM\file2310254bcc.tif --calc="amin([A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X], axis=0)" --co="COMPRESS=LZW"

0.. evaluation of calculation amin([A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X], axis=0) failed
Traceback (most recent call last):
  File "c:\Python33\lib\site-packages\numpy\core\fromnumeric.py", line 2208, in amin
    amin = a.min
AttributeError: 'list' object has no attribute 'min'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Python33\Scripts\GDAL_C~1.PY", line 329, in <module>
    main()
  File "C:\Python33\Scripts\GDAL_C~1.PY", line 326, in main
    doit(opts, args)
  File "C:\Python33\Scripts\GDAL_C~1.PY", line 275, in doit
    myResult = eval(opts.calc)
  File "<string>", line 1, in <module>
  File "c:\Python33\lib\site-packages\numpy\core\fromnumeric.py", line 2211, in amin
    out=out, keepdims=keepdims)
  File "c:\Python33\lib\site-packages\numpy\core\_methods.py", line 29, in _amin
    return umr_minimum(a, axis, None, out, keepdims)
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

同样，在删除对 X.

的引用时，它按我预期的方式运行

我在 Win 10 x64 上使用 Python 3.3.5、gdal 1.11.1 和 numpy 1.9.1。

为什么会这样？

Answer 1

您运行遇到的问题是来自 calc 输入的 X 变量与 doit 函数中循环创建的变量发生冲突：

for X in range(0,nXBlocks):

gdal 开发人员似乎已经修复了这个问题（在尚未发布的代码中）。它现在 eval 是私有名称空间（字典）中的 calc 输入，而不是函数的本地名称空间中的输入，因此循环中的 X 变量不会与 X 来自输入的值。

gdal_calc amin 在传递超过 23 个输入文件时失败

gdal_calc amin fails when passing more than 23 input files

python

numpy

r

raster

gdal