# Processing large amounts of Rasters in R

I have 6500 processed and clipped rasters that indicate snow levels on a daily basis in the northern hemisphere. I need to create a raster that represents a sum of the snow for the season. The winter season would then go from September 1999 to march 2000. I was just wondering how I would go about recalling the daily rasters into R using the loop function/ how i would go about calculating the seasonal sums. Keep in mind that the rasters span two different years and that I don't want to include all 1999 snow data because that would include the previous winter(98/99) as well.

• What is the naming system of your files? Do you have to deal with missing files? – Mikkel Lydholm Rasmussen Jun 12 '15 at 18:39

Index each day with its year and julian day. September 1 is 244, March 31 is 90. Not knowing the format of your files, you can probably figure out a date format to use to parse and turn it into a POSIX date (from https://stackoverflow.com/questions/21414847/convert-a-date-vector-into-julian-day-in-r):

tmp <- as.POSIXlt("16Jun10", format = "%d%b%y")
tmp$yday  Assuming you have extracted your dates from filenames: tmp <- as.POSIXlt(myFormattedDates, format = "%d%b%y") julianDay <- tmp$yday


Next, to determine the winter year for each file, take the year and julian day as a 7-digit number (e.g., September 1, 1998 is 1998244 and March 31, 1999 is 1999090), subtract 91, divide by 1000, and take the floor:

winterYear <- floor(((year*1000+julianDay)-91)/1000)


this assumes you have three vectors with the same length, corresponding to the same rasters:

1. a vector year containing the year of all the files
2. a vector julianDay containing the julian day of all the files
3. a vector rasterNames with the paths to all the rasters.

Next, you can select the files based on the given year:

for (year in unique(winterYear)) {
winterSumRaster <- NULL
thisWinterRasterNames <- rasterNames[which(winterYear == year)]
for (winterRasterName in thisWinterRasterNames) {
dailyRaster <- raster(winterRasterName)
if (is.null(winterSumRaster)) {
winterSumRaster <- raster(dailyRaster)
setValues(winterSumRaster,getValues(dailyRaster))
} else {
setValues(winterSumRaster,getValues(dailyRaster)+getValues(winterSumRaster))
}
}
writeRaster(winterSumRaster,paste0("winterSum",year,".grd"))
}


You can save the raster in other formats of course, but .grd does not require anything special like gdal.

• Nice answer but why the getValues, setValues approach? Why not just read all the data as a stack, use which on the names to create a subset index and then use calc with sum to derive seasonal sums? – Jeffrey Evans Jun 12 '15 at 19:20
• @JeffreyEvans, I would like to see what you mean. It sounds more elegant than my solution. – aaryno Jun 12 '15 at 19:59
• setValues(winterSumRaster,getValues(dailyRaster)) - this is not actually persisted anywhere since there's no assigment – mdsumner Jun 12 '15 at 20:49
• @mdsumner what's the appropriate method call there? – aaryno Jun 12 '15 at 22:17
• I think the call is ok, you just want to assign it somewhere, .i.e. winterSumRaster <- . – mdsumner Jun 13 '15 at 0:37