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): plot

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?

6
  • 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?
    – Spacedman
    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).
    – Thomas
    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.
    – Spacedman
    Sep 29, 2021 at 18:14
  • Here is the data: drive.google.com/file/d/1aUG0fx4BJ9SrJvUhUcXuLj5_Vpt4mBYA/…
    – Thomas
    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.
    – Thomas
    Sep 30, 2021 at 6:23

1 Answer 1

1

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.

enter image description here

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

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

enter image description here

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)

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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