¿Lectura y reclasificación eficiente de muchos rasters en R?

Question

Me han encargado que cree un análisis de idoneidad de las condiciones de las olas en el Golfo de México. Tengo unos 2.000 archivos ráster de unos 8 MB cada uno (2.438 columnas, 1.749 filas, tamaño de celda de 1 km). El parámetro de los archivos ráster es el período de las olas y me gustaría reclasificar todos los rásteres de manera que si 4<= wave period <=9 then make cell = 1 A continuación, la suma de todos los raster en un raster final y dividir por el número total de raster para dar un porcentaje total de las observaciones adecuadas y el resultado de la exportación en algún formato compatible con ESRI ... tal vez algo que puede apoyar flotadores si es necesario. No he trabajado mucho con Python o R, pero después de buscar en línea parece tener sentido para llevar a cabo este proceso en una de esas lenguas. He llegado con un poco de código hasta ahora en R, pero estoy confundido acerca de cómo hacer este trabajo.

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 debo también reclasificar y empezar a sumar los datos dentro de este bucle o no. Mi opinión es que sí, ya que de lo contrario creo que posiblemente me quedaría sin memoria en mi ordenador, pero no estoy seguro. Tampoco estoy seguro de cómo reclasificar los datos.

Buscando en Internet creo que utilizaría reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)) pero, ¿se ve bien?

Y para la suma utilizaría sum(stack(my_data)) y de alguna manera la salida que. Además... si esto se puede realizar de forma más eficiente o escribir en Python, también estaría abierto a ello. Realmente soy un principiante cuando se trata de la programación.

Answer 1

Esta es una forma concisa de hacerlo 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')

Answer 2

Como he notado en los comentarios, generalmente se 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í tienes un ejemplo de código 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 conoce cómo obtener todos los archivos del directorio . Tenga en cuenta que los rásteres se representan como arrays, por lo que puede utilizar métodos de arrays para mejorar el rendimiento cuando sea aplicable. Para trabajar con arrays en Python vea Documentación sobre Numpy .

UPD: el código que publiqué inicialmente era una versión truncada de un filtro personalizado que necesitaba un 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)

Answer 3

1. No utilice readGDAL. Lee en un objeto Spatial* lo que podría no ser una buena idea..

2. Utilice el raster paquete. Puede leer cosas GDAL en objetos Raster. Esto es algo bueno. r = raster("/path/to/rasterfile.tif") lo leerá en r .

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

4. La gran pregunta es si hay que imprimirlos en nuevos archivos de trama y luego hacer el paso de resumen en otro bucle. Si tienes capacidad de almacenamiento, yo escribiría en nuevos archivos porque si tu bucle falla (porque uno de esos 2000 archivos es basura) tendrás que empezar de nuevo. Así que usa writeRaster para crear rásteres con umbral si decide hacerlo.

5. Su bucle es entonces algo así como

esto.

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 gestión de la memoria de R podría golpear aquí - cuando se hace count=count+r R bien podría hacer una nueva copia de count . Así que eso es potencialmente 2000 copias de ella - pero la recolección de basura se 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, la operación de umbral tarda aproximadamente 1,5s en una trama de ese tamaño de números aleatorios entre cero y veinte. Por 2000 = haz las cuentas. Como siempre, haz un pequeño conjunto de pruebas para desarrollar el código, luego lánzalo sobre tus datos reales y vete a tomar un café.