9

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.

3
  • Just use python-gdal. It will be much more efficient in your case. May 4, 2013 at 18:27
  • Thanks, Rebelious. Just curious to know why python-gdal is more efficient in this situation? Also would it be possible to see some of the steps in code that would be necessary to do this? Trying to figure out how best to do this while utilizing as little memory, and cpu, as possible...it's confusing to figure out how to write the code such that it will read in the data, process, take it out of memory and then move to the next raster.
    – Nigel
    May 4, 2013 at 22:27
  • I can't tell you exactly why, but the general cause is that R was designed for other purposes and is known to perform slow with cycles. As for code example, if no one will provide it, I will share one with you in about 10 hours when I will get access to the machine where corresponding script is stored. May 5, 2013 at 5:14

3 Answers 3

10

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')
4
  • 1
    +1 This is good for small problems, but let's do the math for this one: 2438 columns times 1749 rows times 8 bytes/value times 2 thousand grids = 63.6 GB, all of which R must keep in RAM to create s. (Likely twice as much RAM is needed because s probably won't replace raster_data.) I hope you have loads of RAM! Your solution can be made practicable by breaking the 2000 grids into smaller groups, performing the calculation for each group, and then combining those calculations.
    – whuber
    May 11, 2013 at 14:54
  • 3
    @whuber: 's' is a small object, just a bunch of pointers to the files. The calc function, like other functions in the raster package, will not load all the values into memory; it will process them in chunks. That is, the breaking up into groups as you suggest is done automatically, behind the scenes. The chunk size can be optimized for the amount of RAM available via rasterOptions(). May 11, 2013 at 17:46
  • 1
    Thank you for explaining that! I had assumed, without checking, that stack and calc worked like most other R functions by first loading all data into RAM.
    – whuber
    May 12, 2013 at 14:25
  • +1 Liking the concision of R vs the provided Python example... May 14, 2013 at 15:29
5

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)
7
  • 5
    The "R is slow when doing cycles (loops)" is often misused as a reason to avoid R. Yes, if you looped over the cells of a raster in R it would be slow, but the raster package works on entire rasters at once, and has a lot of C code and so runs at C speed. For a raster that size, most of the work would be at C speeds, the looping overhead would be insignificant.
    – Spacedman
    May 5, 2013 at 18:28
  • @Spacedman, yes 'raster' is a useful package (and I do like it), but I never was satisfied with its performance even when loops were not involved. May 6, 2013 at 6:50
  • 2
    Okay, well compare the time it takes in R with the time it takes in Python. Can you not operate on a whole numpy array rather than looping?
    – Spacedman
    May 6, 2013 at 7:02
  • @Spacedman, I just updated the answer. May 6, 2013 at 7:10
  • Many thanks to you both. I'm going to try tinkering with both the Python code you provided and some R and see what I can accomplish. I'll update with results or issues.
    – Nigel
    May 6, 2013 at 21:16
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.

1
  • I'm curious about what you mean by "the threshold op." On my system (running R 2.15.0), reading a 1.6 megapixel grid (in native ESRI floating point format) takes 0.19 seconds and performing the five loop operations--two comparisons, a conjunction, an addition, and a store--takes another 0.09 seconds, for 0.28 seconds per iteration. When instead the loop is performed by caching 100 grids at a time into an array and using rowSums to do the totals, RAM usage of course goes up: but, equally obviously, there is no garbage collection. Timing drops only to 0.26 seconds per iteration.
    – whuber
    May 6, 2013 at 13:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.