¿Por qué la agregación de píxeles MODIS no suma 1?

Question

Quiero encontrar la cobertura fraccional de diferentes especies de árboles a partir de MODIS (30 m de resolución). Utilizo R para apilar, agregar (1020 m de resolución) y así encontrar la cobertura fraccional de especies arbóreas para una región determinada. El script está aquí:

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.

Los valores únicos son 1,2,3, por lo que obtengo 4 capas (una capa adicional sin sentido (todo es 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)})

Luego los apilo, y el resultado final no suma 1 por píxel cuando lo abro en Qgis.

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

No todos los píxeles suman 1 cuando se suman las fracciones. He mirado la capa sin sentido, pero sólo está formada por 0's. ¿Por qué no puedo conseguir que la cobertura fraccional de especies arbóreas de cada píxel sume 1 en total? Estos son los datos que utilizo (trazados):

EDITAR: Aquí están los datos: árbol

El problema se produce cuando sumo los píxeles y los divido con 1156, porque algunos de los píxeles que se incluyen tienen NA valores. Por lo tanto, nunca llegaré a una suma de 1.

¿Puedo tenerlo en cuenta de alguna manera?

Answer 1

Primero estoy trabajando en lo que espero que sea un subconjunto representativo para fallar más rápido (y, por tanto, obtener la solución más rápido). Deberías hacer esto siempre. No trabajes con un conjunto de datos completo hasta que tengas una solución en un conjunto de datos pequeño:

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

A continuación, se realiza un bucle sobre los valores de unique(r) que, a diferencia de unique(r[]) no incluye el NA datos, agregar por un factor de 34 en la trama que es igual al valor:

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

Por defecto aggregate utiliza mean por lo que devolverá la fracción de cada una de las celdas agregadas 34x34 igual a cada valor.

Ahora bien, ¿suman 1? Veamos:

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

Casi siempre, excepto cuando faltan datos, y esas motas de impar. Esos son sólo donde el redondeo aritmético está muy cerca de 1...

plot(suma > 0.999999999999999)