Lendo e reclassificando com eficiência muitos rasters em R?

Fui encarregado de criar uma análise de adequação das condições das ondas no Golfo do México. Eu tenho mais ou menos 2.000 arquivos raster que têm cerca de 8 MB cada (2438 colunas, 1749 linhas, tamanho de célula de 1 km). O parâmetro nos arquivos de varredura é o período de onda e eu gostaria de reclassificar todos os rasters, de modo que 4<= wave period <=9 then make cell = 1, se , outra célula = 0. Em seguida, some todos os rasters em uma varredura final e divida pelo número total de rasters para dar uma porcentagem total de observações adequadas e resultado da exportação para algum formato compatível com ESRI ... talvez algo que possa suportar flutuadores, se necessário. Eu não trabalhei muito com Python ou R, mas depois de pesquisar on-line, parece fazer sentido executar esse processo em uma dessas linguagens. Eu criei algum código até agora no R, mas estou confuso sobre como fazer isso funcionar.

library(rgdal)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
num_files <- length(raster_data)
for (i in raster_data) {                   #read in rasters
   my_data <- readGDAL(raster_data[i])

Nesse ponto, estou confuso quanto a reclassificar e começar a somar os dados dentro desse loop ou não. Meu palpite seria sim, pois caso contrário, acho que ficaria sem memória no meu computador, mas não tenho certeza. Também não tenho certeza sobre como reclassificar os dados.

Ao pesquisar online, acho que usaria reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), mas isso parece certo?

E, para resumir, eu usaria sum(stack(my_data))e, de alguma forma, exibirei isso. Além disso ... se isso puder ser executado ou escrito com mais eficiência em Python, eu também estaria aberto a isso. Eu realmente sou iniciante quando se trata de programação.

Esta é uma maneira concisa de fazer isso no R --- aqui sem arquivos intermediários:

library(raster)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
s <- stack(raster_data)
f <- function(x) { rowSums(x >= 4 & x <= 9) }
x <- calc(s, f, progress='text', filename='output.tif')

Como observei nos comentários, geralmente você deve evitar o uso de R para fins não estatísticos devido a problemas de desempenho em certos aspectos (trabalhar com ciclos é um exemplo). Aqui está um exemplo de código para você em Pyhton (graças a este artigo ) para reclassificação de um único arquivo com uma única banda. Você poderá modificá-lo facilmente para processamento em lote se souber como obter todos os arquivos do diretório . Observe que os rasters são representados como matrizes; portanto, você pode usar métodos de matriz para melhorar o desempenho, quando aplicável. Para trabalhar com matrizes em Python, consulte a documentação do Numpy .

UPD: o código que publiquei inicialmente era uma versão truncada de um filtro personalizado necessário para o processamento de pixels. Mas para esta pergunta, o uso do Numpy aumentará o desempenho (consulte o código).

from osgeo import gdal
import sys
import numpy

gdalData = gdal.Open("path_to_file")
if gdalData is None:
  sys.exit("ERROR: can't open raster")

#print "Driver short name", gdalData.GetDriver().ShortName
#print "Driver long name", gdalData.GetDriver().LongName
#print "Raster size", gdalData.RasterXSize, "x", gdalData.RasterYSize
#print "Number of bands", gdalData.RasterCount
#print "Projection", gdalData.GetProjection()
#print "Geo transform", gdalData.GetGeoTransform()


raster = gdalData.ReadAsArray()
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize

#print xsize, 'x', ysize

## per pixel processing
#for col in range(xsize):
#  for row in range(ysize):
#    # if pixel value is 16 - change it to 7
#    if raster[row, col] == 16:
#      raster[row, col] = 7

#reclassify raster values equal 16 to 7 using Numpy
temp = numpy.equal(raster, 16)
numpy.putmask(raster, temp, 7)

# write results to file (but at first check if we are able to write this format)
format = "GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
if metadata.has_key(gdal.DCAP_CREATE) and metadata[gdal.DCAP_CREATE] == "YES":
  pass
else:
  print "Driver %s does not support Create() method." % format
  sys.exit(1)
if metadata.has_key(gdal.DCAP_CREATECOPY) and metadata[gdal.DCAP_CREATECOPY] == "YES":
  pass
else:
  print "Driver %s does not support CreateCopy() method." % format
  sys.exit(1)

# we already have the raster with exact parameters that we need
# so we use CreateCopy() method instead of Create() to save our time
outData = driver.CreateCopy("path_to_new_file", gdalData, 0)
outData.GetRasterBand(1).WriteArray(raster)

1. Não use readGDAL. Ele lê um objeto Spatial * que pode não ser uma boa ideia.

2. Use o rasterpacote. Ele pode ler coisas GDAL em objetos Raster. Isso é uma coisa boa. r = raster("/path/to/rasterfile.tif")irá lê-lo r.

3. Sua classificação é então t = r > 4 & r <= 9

4. A grande questão é se esses arquivos serão exibidos em novos arquivos de varredura e, em seguida, execute a etapa de resumo em outro loop. Se você tiver o armazenamento, eu os escreveria em novos arquivos apenas porque se o seu loop falhar (porque um desses arquivos de 2000 é lixo), você terá que iniciar novamente. Portanto, use writeRasterpara criar rasters limiares se você decidir fazer isso.

5. Seu loop é então apenas algo como

isto.

count = raster(0,size of your raster)
for(i in 1:number of rasters){
  r = raster(path to binary raster file 'i')
  count = count + r
}
return(count)

O gerenciamento de memória de R pode ser útil aqui - quando você faz count=count+rR, pode fazer uma nova cópia de count. Então, são potencialmente 2000 cópias dele - mas a coleta de lixo é ativada e limpa, e é aqui que o loop de R costumava ser muito ruim.

Em termos de tempo, no meu PC de 4 anos, o limiar de operação demora cerca de 1,5s em uma varredura desse tamanho de números aleatórios entre zero e vinte. Vezes 2000 = você faz as contas. Como sempre, faça um pequeno conjunto de testes para desenvolver o código, inicie-o nos dados reais e tome um café.