Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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. Any help or suggestions for a better method would be very helpful.

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. Thanks for any assistance. I truly am a beginner when it comes to programming.

share|improve this question
    
Just use python-gdal. It will be much more efficient in your case. –  SS_Rebelious May 4 '13 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 '13 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. –  SS_Rebelious May 5 '13 at 5:14
    
That'd be super helpful, thanks again. –  Nigel May 5 '13 at 11:29

3 Answers 3

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)
share|improve this answer
2  
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 '13 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. –  SS_Rebelious May 6 '13 at 6:50
    
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 '13 at 7:02
    
@Spacedman, I just updated the answer. –  SS_Rebelious May 6 '13 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 '13 at 21:16
  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.

share|improve this answer
    
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 '13 at 13:52

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')
share|improve this answer
    
+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 '13 at 14:54
    
@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(). –  RobertH May 11 '13 at 17:46
    
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 '13 at 14:25
    
+1 Liking the concision of R vs the provided Python example... –  SlowLearner May 14 '13 at 15:29

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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