4

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
5

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.

  • 2
    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
  • 1
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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