gdal_calc amin non riesce quando si superano più di 23 file di input

Question

Ho scritto una funzione R che chiama gdal_calc.py per calcolare il valore minimo in pixel su un RasterStack (serie di file raster di input). L'ho fatto perché è molto più veloce di raster::min per i raster di grandi dimensioni. La funzione funziona bene per un massimo di 23 file, ma genera un avvertimento quando passa 24 o più, restituendo un raster di output pieno di zeri.

Poiché R sta appena preparando una chiamata di sistema a python gdal_calc.py, questa domanda non è specifica per R e incoraggio gli appassionati di python/numpy a leggere.

Ecco la funzione. La struttura dell'eventuale chiamata gdal_calc viene visualizzata nel messaggio di avviso generato dall'utilizzo problematico nella parte inferiore di questo 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) 
} 

E qui è una griglia raster manichino, scritto 'test.tif' nella directory di lavoro attuale. Per semplicità, e per dimostrare che il problema non è dovuto a un particolare raster individuale, passo lo stesso file di input più volte a gdal_calc, piuttosto che passare un numero di file diversi.

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

La funzione funziona bene per i file di input 23:

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

... ma non per ≥ 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 

---

Il traceback dato durante l'esecuzione direttamente quest'ultimo dal terminale:

"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() 

Anche in questo caso, viene eseguito come mi aspettavo durante la rimozione del riferimento a X.

Sto usando Python 3.3.5, gdal 1.11.1 e numpy 1.9.1 su Win 10 x64.

Perché ciò potrebbe accadere?

Answer 1

Il problema si sta eseguendo in è stata che la variabile X dall'ingresso calc è stata la collisione con la variabile creata da un ciclo nella funzione doit:

for X in range(0,nXBlocks): 

Risulta questo è già stato fissato dal GDAL sviluppatori (nel codice non ancora rilasciato). Ora eval s l'input calc in uno spazio dei nomi privato (un dizionario), piuttosto che nel namespace locale della funzione, quindi la variabile X dal ciclo non è più in conflitto con il valore X dall'input.