Your code isn't computing monthly averages. Its doing something very very different.
The numerical argument to the stackApply function defines the grouping that is passed into the function. This numerical argument is repeated to the length of the number of layers. If you have 36 layers, then the numeric argument is 1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5... etc. It then computes the mean of all the "1" layers, the mean of all the "2" layers and so on up to the mean of all the "12" layers. In your output "January" is actually the mean of layer 1, layer 13 and layer 25, "February" is the mean of layers 2, 14 and 26.
Here's a little test. I'll make two rasters, one with values 1 to 12 and one with all zeroes:
> r = raster(matrix(1:12, 3,4))
> r0 = raster(matrix(0, 3,4))
and stack them, first 12 are 1:12 and second 12 are zeroes:
> s = stack(r,r,r,r,r,r,r,r,r,r,r,r,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0,r0)
Then stackApply(s, 1:12, mean) returns:
> stackApply(s, 1:12, mean)
class : RasterBrick
dimensions : 3, 4, 12, 12 (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.3333333 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : index_1, index_2, index_3, index_4, index_5, index_6, index_7, index_8, index_9, index_10, index_11, index_12
min values : 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
max values : 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6
and you see all 12 output layers are identical because they are the mean of one r layer and one r0 layer. The first layer is mean of layer 1 and layer 13, and so on up to layer 12 and 24
What you seem to need to do is group the first 3 layers, then the next three, then the next three and so on. To do this you need to construct a vector of three ones, three twos, three threes and so on.
My example above has 24 layers, which in your scheme with three layers per month would be 8 months. I can construct a grouping vector eg using rep:
> grp = rep(1:8, rep(3,8))
> grp
[1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8
Then I can stackApply with that grouping vector to get my 8 monthly means:
> stackApply(s, grp, mean)
class : RasterBrick
dimensions : 3, 4, 12, 8 (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.3333333 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : index_1, index_2, index_3, index_4, index_5, index_6, index_7, index_8
min values : 1, 1, 1, 1, 0, 0, 0, 0
max values : 12, 12, 12, 12, 0, 0, 0, 0
and you can see the first four months are 1:12 and the last four months are all zero, as expected.
A safer way to do this which would be robust to any missing monthly data would be to derive the grouping from the data itself, for example from the file name. If you can get a categorical value you can group on that. For example, if you have 2006 monthly data and can compute this:
> grp
[1] "200601" "200601" "200601" "200602" "200602" "200602" "200603" "200603"
[9] "200603" "200604" "200604" "200604" "200605" "200605" "200605" "200606"
[17] "200606" "200606" "200607" "200607" "200607" "200608" "200608" "200608"
Then stackApply will still work and as a bonus you get more meaningful layer names:
> stackApply(s, grp, mean)
class : RasterBrick
dimensions : 3, 4, 12, 8 (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.3333333 (x, y)
extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : index_200601, index_200602, index_200603, index_200604, index_200605, index_200606, index_200607, index_200608
min values : 1, 1, 1, 1, 0, 0, 0, 0
max values : 12, 12, 12, 12, 0, 0, 0, 0
And as a final point you should probably look at using the terra package for raster handling instead of raster. The terra package is newer, under development, but does all this differently (the relevant function is tapp, I think).
Best Answer
I've made a folder with 162 tiffs in it that should roughly replicate your data. Here's the first 10:
When I make the stack, I get the year-month-date from the filename in the names:
with an X to make it non-numeric. Now my grouping is going to be based on the year-month, which is the first 7 characters of the name (including the X, we'll keep that):
Now apply means over the groupings in the stack. The resulting stack has useful names:
So we can loop over the names, extract the layer, and write it using the name:
and that gives me 54 individual single-layer TIFF files output:
Use whatever string processing you want to create according to whatever naming convention you want to use.
Also, try and use the
terrapackage.