¿Leyendo y reclasificando eficientemente muchos rásteres en R?

Me encargaron crear un análisis de idoneidad de las condiciones de las olas en el Golfo de México. Tengo aproximadamente 2 mil archivos ráster que son de aproximadamente 8 MB cada uno (2438 columnas, 1749 filas, 1 km de tamaño de celda). El parámetro en los archivos ráster es el período de onda y me gustaría reclasificar todos los rásteres de modo que si 4<= wave period <=9 then make cell = 1, de lo contrario, la celda = 0. Luego, sume todos los rásteres en un ráster final y divida por el número total de rásteres para obtener un porcentaje total de observaciones adecuadas y resultados de exportación en algún formato compatible con ESRI ... tal vez algo que pueda admitir flotantes si es necesario. No he trabajado mucho con Python o R, pero después de buscar en línea parece tener sentido realizar este proceso en uno de esos idiomas. Se me ocurrió un código hasta ahora en R, pero estoy confundido sobre cómo hacer que esto funcione.

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])

En este punto, estoy confundido sobre si también debería reclasificar y comenzar a sumar los datos dentro de este ciclo o no. Supongo que sí, ya que de lo contrario creo que posiblemente me quedaría sin memoria en mi computadora, pero no estoy seguro. Tampoco estoy seguro de cómo reclasificar los datos.

Al investigar en línea, creo que usaría reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), pero ¿se ve bien?

Y para resumir, usaría sum(stack(my_data))y de alguna manera generaría eso. Además ... si esto pudiera realizarse o escribirse de manera más eficiente en Python, también estaría abierto a eso. Realmente soy un principiante cuando se trata de programación.

Esta es una manera concisa de hacer eso en R --- aquí sin archivos intermedios:

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 noté en los comentarios, generalmente debe evitar el uso de R para fines no estadísticos debido a problemas de rendimiento en ciertos aspectos (trabajar con ciclos es un ejemplo). Aquí hay un ejemplo de código para usted en Pyhton (gracias a este artículo ) para la reclasificación de un solo archivo con una sola banda. Podrá modificarlo fácilmente para el procesamiento por lotes si sabe cómo obtener todos los archivos del directorio . Tenga en cuenta que los rásteres se representan como matrices, por lo que puede usar métodos de matriz para mejorar el rendimiento cuando corresponda. Para trabajar con matrices en Python, consulte la documentación de Numpy .

UPD: el código que publiqué inicialmente era una versión truncada de un filtro personalizado que necesitaba procesamiento por píxel. Pero para esta pregunta, el uso de Numpy aumentará el rendimiento (ver 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. No uses readGDAL. Se lee en un objeto espacial * que podría no ser una buena idea.

2. Usa el rasterpaquete. Puede leer cosas GDAL en objetos Raster. Estas son una buena cosa. r = raster("/path/to/rasterfile.tif")lo leerá en r.

3. Tu clasificación es entonces t = r > 4 & r <= 9

4. La gran pregunta es si enviarlos a nuevos archivos ráster y luego hacer el paso de resumen en otro bucle. Si tiene el almacenamiento, los escribiría en archivos nuevos solo porque si su bucle falla (porque uno de esos 2000 archivos es basura) tendrá que comenzar de nuevo. Por lo tanto, use writeRasterpara crear rásteres trillados si decide hacerlo.

5. Tu ciclo es entonces algo así como

esta.

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)

La administración de la memoria de R podría afectarlo aquí; cuando lo haga, count=count+rR bien podría hacer una nueva copia de count. Por lo tanto, eso es potencialmente 2000 copias, pero la recolección de basura se activará y limpiará, y es aquí donde el bucle de R solía ser muy malo.

En términos de tiempo, en mi PC de 4 años, el umbral de operación toma alrededor de 1.5 segundos en una trama de ese tamaño de números aleatorios entre cero y veinte. Veces 2000 = haces los cálculos. Como siempre, haga un pequeño conjunto de prueba para desarrollar el código, luego ejecútelo en sus datos reales y tome un café.