# Why does aggregation of MODIS pixels not sum to 1?

I want to find the fractional cover of different tree species from MODIS (30 m resolution). I use R to stack, aggregate (1020 m resolution) and thereby find the fractional tree species cover for a given region. The script is here:

``````r = raster("path/tree.tif")
s = do.call(stack, lapply(unique(r[]), function(v){(r==v)*v})) #Split raster into unique layers in a stack.
``````

The unique values are 1,2,3, so I get 4 layers (an additional layer of nonsense (everything is 0)).

``````#Aggregate MODIS pixel (30m) to a ~1km resolution. Sum all values of 1 (x == 1) and divide by the total number of pixels being aggregated (34*34)
vegetation1 = aggregate(x = s[[3]], fact = 34, fun = function(x, ...){(sum(x == 1, ...)/1156)})
vegetation2 = aggregate(x = s[[4]], fact = 34, fun = function(x, ...){(sum(x == 2, ...)/1156)})
vegetation3 = aggregate(x = s[[2]], fact = 34, fun = function(x, ...){(sum(x == 3, ...)/1156)})
``````

I then stack them, and the end result does not add up to 1 per pixel when I open it in Qgis.

``````stacked <- stack(veg1,veg2,veg3)
writeRaster(stacked, "path/stacked_tree.tif")
``````

Not all pixels sum up to 1 when the fractions are added together. I have looked at the nonsense layer, but it only consists of 0's. Why can I not get every pixel's fractional tree species cover, to sum up to 1 in total? This is the data I use (plotted):

EDIT: Here is the data: tree

The problem occurs when I sum up the pixels and divide them with 1156, because some of the pixels that are included have `NA` values. Therefore I will never make it to a sum of 1.

Can I somehow account for this?

• Should that be `lapply` in your first code block instead of `apply`? Why are you aggregating the cells into 34x34 cells groups and where does 1156 come from? Is your data 34x34 cells (I see 1156 is 34*34)... Why not do `table(r[])` to get the number of cells of each value? Sep 29, 2021 at 15:44
• Does the updates help? Yes it is `lapply` and the numbers come from aggregating 30 m resolution to a 1020 m resolution. Thus, multiply 30 m with a factor 34, and then divide with the total number of cells that have been aggregated (34x34 = 1156). Sep 29, 2021 at 16:33
• can you supply the source tree.tif file? or something similar that shows the problem? or code to generate a sample raster that shows the problem? currently we don't know anything about the data. Sep 29, 2021 at 18:14
• Here is the data: drive.google.com/file/d/1aUG0fx4BJ9SrJvUhUcXuLj5_Vpt4mBYA/… Sep 30, 2021 at 6:21
• If you run the script, you will see that especially the layer with 0's gives some problems. I end up having band 1 summing up to 1 and the other bands still having some fractional tree cover too, thus exceeding 1. Sep 30, 2021 at 6:23

## 1 Answer

First I'm working on what I hope is a representative subset in order to fail faster (and hence get the solution faster). You should always do this. Don't work with a full dataset until you have a solution on a small dataset:

``````e = extent(c(xmin=-550000, xmax=-400000, ymin=2020000, ymax=2170000))
r = raster("tree.tif")
r = crop(r, e)
``````

Then looping over the values from `unique(r)`, which, unlike `unique(r[])` doesn't include the `NA` data, aggregate by a factor of 34 on the raster being equal to the value:

``````aas = do.call(stack, lapply(unique(r), function(x){aggregate(r==x, 34)}))
``````

By default `aggregate` uses `mean`, so this will return the fraction of each of the aggregated 34x34 cells equal to each value.

Now, do they add up to 1? Let's see:

``````suma = stackApply(aas, 1, sum)
plot(suma == 1)
``````

Mostly, except where there's missing data, and those odd speckles. Those are just where arithmetic rounding is very close to 1...

``````plot(suma > 0.999999999999999)
``````

• I accept your answer @Spacedman, however, I have simply forgotten to reclassify some of the data and that caused the (in hindinsight) obvious problem. Oct 1, 2021 at 10:08