I'm working on global climate projections raster maps. My starting data are global raster maps of daily future climate, where every map represent one single day (e.g. I have 365 maps for 1 year). Furthermore, for every cell in the single map I have a temperature value.

I want to create maps of indices of extreme climatic events (ECEs) based on my starting data. An example of ECE indices is "number of days over 25°C in a year" (SD, Summer Days).

  1. I was able to do it for single maps. Use the "IF" operation in the Raster Calculator. The output is a binary global map for a single day, where 1 are the "t > 20°C" values of the map and 0 are the "t < 20°C" values.
  2. I wasn't able to do it, not even with just two maps. "SUM" the binary maps all together (365 maps for 1 year, but possibly I would like to have a 20 year mean), so that I can get an indices as "number of days over 25°C in a year".

My real questions are:

  • Why the "SUM" formula does not work with just two raster (if the "IF" formula works)? (Less important question, I think I can solve it somehow)
  • Is it possible to compute together the "IF" and the "SUM" operations in the Raster Calculator? (Important question)
  • Is it possible to summarize together all the formulas related to one year (or more years) in just one calculation (for one year, it would mean: 365 "IF" plus 365 "SUM")? (Much Important question)
  • 3
    Dont ask multiple questions in one post, your question might get closed. Focus on the most important one – BERA Sep 27 at 16:28
  • 1
    Welcome to Geographic Information Systems! Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works. – Ian Turton Sep 27 at 16:43

If your files are named t_<yyyy>_<doy>, then you can do this in 2 steps: (1) Compute your ECE index and (2) sum.

(1), using GNU parallel to do it quickly and save writing for loops:

parallel "r.mapcalc \"ece_{1}_{2} = if(t_{1}_{2} > 20)\"" ::: $(seq 2020,2099) ::: $(seq 1 365)

This r.mapcalc should match what it sounds like you already did successfully. Perhaps the statement is a bit different (Using >=?, or if(test, 1, 0) or if(test, 1, null()), etc... Also, the parallel loops through 365 days, I'm not sure how to best handle leap years if they are in your model outputs. Perhaps 366 and ignore errors?

(2) should use r.series:

r.series input=$(g.list type=raster pattern=ece_* separator=comma), output=output, method=sum

The above is inefficient - it computes and stores a whole bunch of new rasters, but they should be small, because they're just 1s and 0s. More efficient is:

r.mapcalc "ece = 0"
for y in $(seq 2020 2099); do
  for doy in $(seq 0 365); do
    r.mapcalc "ece = ece + if(t_${y}_${doy} > 20)" --o

Note that r.mapcalc says you should not use the same variable on the LHS and RHS of an equation. I've actually never run into a bug with this (as far as I know). So perhaps the inner code above should be:

r.mapcalc "tmp = ece + if(t_${y}_${doy} > 20)" --o
g.rename raster=tmp,ece --o

Finally, I don't think this 2nd method can be parallelized, but if you want to use parallel to save yourself writing for loops, you can with the -j 1 flag.

| improve this answer | |
  • Hi mankoff! For step 1. The file are named like this: "tasmax_day_CMCC_CMS_rcp45_r1i1p1_20300101_20391231.1" where the LAST number of name represent the day (1=1/1/2030; 3650 is ten years after, it is 31/12/3039). Formula that works in GRASS calculator is: if( tasmax_day_CMCC_CMS_rcp45_r1i1p1_20300101_20391231.1@PERMANENT >20,1,). How do I put the iteration? Maybe what you are saying is that it is possible to write something as: if( tasmax_day_CMCC_CMS_rcp45_r1i1p1_20300101_20391231.i@PERMANENT >20,1,), for i=1:3650 Can the calculator understand the file I want it to read? Thank you – matteo s. Sep 28 at 18:55
  • I don't use GRASS GUI, only command line interface. I'm not sure how to loop with the GUI. – mankoff Sep 29 at 12:56

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