Efficiently reading and reclassifying many rasters in R?

Question

I've been tasked to create a suitability analysis of wave conditions in the Gulf of Mexico. I have 2 thousand or so raster files that are about 8 MB each (2438 columns, 1749 rows, 1km cell size). The parameter in the raster files is wave period and I'd like to reclassify all of the rasters such that if 4<= wave period <=9 then make cell = 1, else cell = 0. Then sum up all the rasters into a final raster and divide by the total number of rasters to give a total percentage of suitable observations and export result into some ESRI-compatible format...maybe something that can support floats if need be. I've not worked much at all with either Python or R, but after searching online it seems to make sense to perform this process in one of those languages. I've come up with some code so far in R, but am confused about how to make this work.

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

At this point I'm confused as to whether I should also reclassify and start summing the data within this loop or not. My guess would be yes since otherwise I think I would possibly run out of memory on my computer, but just not sure. I'm also not sure about how to reclassify the data.

In researching online I think I would use reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), but does that look right?

And for summing would I use sum(stack(my_data)) and somehow output that. Also...if this might be more efficiently performed or written in Python I'd be open to that as well. I truly am a beginner when it comes to programming.

Answer 1

This is a concise way to do that in R --- here without intermediate files:

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

As I noticed in comments, generally you should avoid using R for non-statistical purposes due to performance issues in certain aspects (working with cycles is an example). Here is code example for you in Pyhton (thanks to this article) for reclassification of a single file with a single band. You will be able to modify it easily for batch processing if you know how to get all files from the directory. Notice that rasters are represented as arrays, so you may use array methods to improve performance when applicable. For working with arrays in Python see Numpy documentation.

UPD: the code I posted initially was a truncated version of a custom filter that needed per pixel processing. But for this question Numpy usage will boost performance (see code).

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. Don't use readGDAL. It reads into a Spatial* object which might not be a good idea..

2. Use the raster package. It can read GDAL things into Raster objects. These are a good thing. r = raster("/path/to/rasterfile.tif") will read it into r.

3. Your classification is then t = r > 4 & r <= 9

4. The big question is whether to output these to new raster files and then do the summary step in another loop. If you've got the storage, I'd write them to new files just because if your loop fails (because one of those 2000 files is junk) you'll have to start again. So use writeRaster to create thresholded rasters if you decide to do that.

5. Your loop is then just something like

this.

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)

R's memory management might hit you here - when you do count=count+r R might well make a new copy of count. So that's potentially 2000 copies of it - but garbage collection will kick in and clean up, and its here that R's looping used to be very bad.

In terms of timing, on my 4yo PC, the threshold op takes about 1.5s on a raster of that size of random numbers between zero and twenty. Times 2000 = you do the maths. As always, make a small test set for developing the code, then launch it on your real data and go have a coffee.