Nesse curso são apresentados os procedimentos básicos para o processamento dos dados de um Inventário Florestal. Esses procedimentos envolvem fundamentalmente as operações de processamento de dados já vistas (nos cursos anteriores) e a aplicação dos modelos dendrométricos para predição de alguns atributos das árvores individualmente.

Ao final desse curso você deverá ser capaz de:

1. organizar os dados de forma que sua estrutura seja compatível com o processamento no R;
2. realizar as operações de agregação dos dados do nível de árvore para o nível de parcela;
3. obter estimativas a nível de povoamento (talhao, propriedade, etc.) num delineamento amostral simples.

Ilustraremos a estrutura dos dados e processamento de inventário com o arquivo exemplo-inventario.csv. Exmploremos um pouco esses dados:

exemplo = read.csv("exemplo-inventario.csv",header=T,as.is=T)
dim(exemplo)
head(exemplo)
colnames(exemplo)
table(exemplo$cod_parcela, exemplo$cod_talhao)

Abaixo mostramos as primeiras linhas desse arquivo na forma de uma tabela: A primeira linha, que está numerada, não existe no arquivo. Ela é apenas uma referência às 17 colunas (variáveis) que compõem o arquivo. <csv> 1,2,3,4,5,6,7,8 “cod_unidade”,“cod_projeto”,“desc_projeto”,“cod_talhao”,“area”,“espacamento”,“clone_semente”,“data_plantio” “CBO”,“F531”,“Ligiana”,“F53101-1996-01”,9.3,300250,“S”,08/02/96 “CBO”,“F531”,“Ligiana”,“F53101-1996-01”,9.3,300250,“S”,08/02/96 “CBO”,“F531”,“Ligiana”,“F53101-1996-01”,9.3,300250,“S”,08/02/96 </csv>

<csv> 9,10,11,12,13,14,15,16,17 “data_medicao”,“clid”,“idade”,“cod_parcela”,“area_parc”,“cod_arvore”,“categoria”,“dap”,“ht” “15/04/2001”,5,“5.182.751”,1,466.8,1,“Normal”,18.46,29 “15/04/2001”,5,“5.182.751”,1,466.8,2,“Normal”,18.08,26.5 “15/04/2001”,5,“5.182.751”,1,466.8,3,“Normal”,16.3,24.3 </csv> Note que da coluna 1 (cod_unidade) até a coluna 11 (idade), todas as variáveis se referem a informação de talhão ou parcela. As coluna 12 (cod_parcela) e 13 (area_parc) trazem informação da parcela e as demais colunas (14 - cod_arvore a 17 - ht) são relativas às árvores individualmente.

A redundância nesse arquivo é muito grande, pois os dados de talhão e parcela são repetidos uma grande quantidade de vezes para acompanhar os dados das árvores. Além disso, para o processamento no nível de árvore, os dados de talhão são raramente necessários (desde que a variável cod_parcela identifique unicamente cada parcela).

Esses dados podem ser melhor organizados na forma de duas tabelas de dados (data frames):

• cadastro - data frame com os dados de talhao e parcela, e
• inventario - data frame com os dados das árvores.

Façamos a re-organização desses dados. Primeiro façamos o cadastro:

cadastro = aggregate(exemplo[ , 1:12], list(exemplo$cod_parcela), function(x) x[1] )
head(cadastro)
cadastro = cadastro[,-1]
head(cadastro)
dim(cadastro)

Façamos agora o data frame inv (inventário):

inv = exemplo[ , c(4,12:17)]
head(inv)
dim(inv)

Note que as variáveis cod_talhao e cod_parcela são redundantes pois estão nos dois data frames. Esse nível de redundância é necessário para podemos ligar novamente os dados de cadastro e o inventário quando necessário.

Será que houve ganho em termos de armazenamento?

dim(exemplo)
prod(dim(exemplo))

dim(cadastro)
prod(dim(cadastro))
dim(inv)
prod(dim(inv))

prod(dim(cadastro)) + prod(dim(inv))

O número de registros caiu de 35.496 (data frame exemplo) para 14.748 (data frame cadastro e inv). Como a maior parte do processamento envolverá apenas o data frame inv, estaremos trabalhando com 14.616 registros na maior parte do tempo.

Na linguagem de banco de dados, os dados agora formam um banco de dados relacional com duas tabelas (data frames) onde as variáveis cod_talhao e cod_parcela estabelecem a relação. Embora seja uma estrutura relacional bastante simples, ela é provavelmente o formato mais eficiente e prático para um inventário florestal simples.

Estudemos um pouco as duas variáveis quantitativas que utilizaremos no processamento: DAP e altura total das árvores:

head(inv)
par(mfrow=c(2,1))
hist(inv$dap, main="dap")
hist(inv$ht, main="ht")

plot(inv$dap, main="dap")
abline(h=0, col="red")
plot(inv$ht, main="ht")
par(mfrow=c(1,1))

Note que nos dados de inventário há uma série de observações ZERO, o que não acontece com a altura. Quantas observações são?

sum( inv$dap == 0 )

Como são elas?

head( inv[ inv$dap == 0, ] )
table( inv$categoria, inv$dap == 0 )

Conclusão: toda Falha (árvore inexistente) recebeu DAP zero. Mas DAP igual a zero é o mesmo que dap inexistente? NÃO. É necessário utilizar a expressão de valor inexistente: NA (not avaiable).

Note que no caso da altura isso aconteceu naturalmente, pois as árvores que não tiveram suas alturas medidas ficaram com o espaço reservado para a variável altura EM BRANCO. Ao ler o arquivo no R, ele automaticamente coloca o rótulo inexistente (NA) nessas observações.

Nota muito importante: a expressão NA não é um valor é a inexistência do valor e, portanto, se aplica a qualquer classe de variável. Note a diferença entre o valor zero e o NA:

x = seq(1, 10, by=0.3)
x
a = 0
a
b = NA
b

a * x
b * x

sum( x )
sum( c(x,a) )
sum( c(x,b) )

prod( x )
prod( c(x,a) )
prod( c(x,b) )

O valor inexistente é INDETERMINADO. Qualquer operação matemática com um valor indeterminado resultará num valor indeterminado. Inclusive as operações lógicas:

x == NA
a == NA
b == NA
sum( inv$ht == NA )
table( inv$categoria, inv$ht == NA )

Mas o R traz a possibilidade de identificar essas observações inexistentes. É a função is.na:

is.na( x )
is.na( a )
is.na( b )

sum( is.na(inv$ht) )
table( inv$categoria, is.na( inv$ht ) )

É necessário converter os dados de dap zero em dados inexistentes para que não tenhamos problemas no processamento de dados (por exemplo no cálculo das médias e desvios padrão).

sum( inv$dap == 0 )
inv$dap [ inv$categoria == "Falha" ] = NA
sum( is.na(inv$dap) )
table( inv$categoria, is.na( inv$dap ) )

Nesse conjunto de dados, existe mais alguma observação que deve ser estudada ou excluída para não gerar resultados equivocados no processamento? Quais? Por que?

Para predição das alturas, utilizaremos um relação hipsométrica por parcela:

require(lattice)
xyplot( ht ~ dap | cod_parcela, data=inv )

relhipso = lm( ht ~ (dap + I(dap^2)) * factor(cod_parcela), data=inv)
plot(relhipso)
summary(relhipso)

Obter a altura predita é simples no R, quando se possui um modelo linear (relhipso):

inv$h.pred = predict(relhipso, inv)
head(inv)

A predição do volume pode partir de um modelo linear ajustado, como no caso da altura, ou a partir de uma equação de volume disponível:

inv$vol = (0.7340 + 0.0350 * tmp2$dap^2 * tmp2$h.pred) * (1/1000)
head(inv)

Nessa etapa, todas as informações no nível de árvores são agregadas para o nível de parcela. Uma grande quantidade de informações podem ser obtidas. Nesse exemplo teremos as mais frequentemente utilizadas.

Primeiramente o fator de expansão da parcela:

inv.par = aggregate(inv$area_parc, list(cod_parcela = inv$cod_parcela), function(x) 10000/x[1] )
colnames(inv.par)[2] = c("fd")
head(inv.par)

Informações de totalização:

# Área Basal
tmp = aggregate(inv$dap, list(cod_parcela = inv$cod_parcela), function(x) sum( (pi/4)*(x[!is.na(x)]/100)^2 ) )
colnames(tmp)[2] = c("g")
inv.par = merge(inv.par, tmp)
inv.par$g = inv.par$g * inv.par$fd
head(inv.par)

# Área Basal
tmp = aggregate(inv$vol, list(cod_parcela = inv$cod_parcela), function(x) sum( x[!is.na(x)] ) )
colnames(tmp)[2] = c("v")
inv.par = merge(inv.par, tmp)
inv.par$v = inv.par$v * inv.par$fd
head(inv.par)

Informações de tamanho médio das árvores:

# DAP
tmp1 = aggregate(inv$dap, list(cod_parcela = inv$cod_parcela), mean, na.rm = TRUE )
colnames(tmp1)[2] = c("mdap")
tmp2 = aggregate(inv$dap, list(cod_parcela = inv$cod_parcela), sd, na.rm = TRUE )
colnames(tmp2)[2] = c("devdap")
tmp = merge(tmp1, tmp2)
inv.par = merge(inv.par, tmp)
head(inv.par)

# Altura
tmp1 = aggregate(inv$ht, list(cod_parcela = inv$cod_parcela), mean, na.rm = TRUE )
colnames(tmp1)[2] = c("mht")
tmp2 = aggregate(inv$ht, list(cod_parcela = inv$cod_parcela), sd, na.rm = TRUE )
colnames(tmp2)[2] = c("devht")
tmp = merge(tmp1, tmp2)
inv.par = merge(inv.par, tmp)
head(inv.par)

# Altura Médias das Árvores Dominantes (aproximadamente 5 por parcela)
tmp = aggregate(inv$ht, list(cod_parcela= inv$cod_parcela), function(x) mean(sort(x[!is.na(x)], decreasing=T )[1:5])  )
colnames(tmp)[2] = c("mhdom")
inv.par = merge(inv.par, tmp)
head(inv.par)

Para obtenção das estimativas necessitamos escolher o nível de estimação. Nesse exemplo utilizaremos o talhão:

head(cadastro)
inv.cad = cadastro[ , c("cod_talhao", "cod_parcela")]
head(inv.cad)
inv.par = merge(inv.cad, inv.par)
head(inv.par)

Obtendo o número de parcelas por talhão

inv.n = aggregate(inv.par$cod_parcela, list(cod_talhao = inv.par$cod_talhao), function(x) length(unique(x)) )
colnames(inv.n)[2] = c("n")
head(inv.n)

Obtendo o valor médio dos atributos das parcelas por talhão

inv.med = aggregate( inv.par[ , c(4:6,8,10)], list(cod_talhao = inv.par$cod_talhao), mean )
varnames = colnames(inv.med)[-1] 
varnames
colnames(inv.med)[-1] = paste(colnames(inv.med)[-1], "media", sep=".")
inv.med

Obtendo o valor desvio padrão dos atributos das parcelas por talhão

inv.dev = aggregate( inv.par[ , c(4:6,8,10)], list(cod_talhao = inv.par$cod_talhao), sd )
colnames(inv.dev)[-1] = paste(colnames(inv.dev)[-1], "devpad", sep=".")
inv.dev

Obtendo o valor intervalo de confiança de 95% dos atributos das parcelas por talhão

tmp1 = inv.dev[,1] 
tmp2 = inv.dev[,-1] * qt(0.975, inv.n[,-1] -1)
inv.ic = cbind(tmp1,tmp2)
colnames(inv.ic)[1] = "cod_talhao"
colnames(inv.ic)[-1] = paste(varnames, "ic95", sep=".")
inv.ic

Estimativas por Talhao:

out = merge(inv.n, inv.med)
out = merge(out, inv.dev)
out = merge(out, inv.ic)
out

João Luís Ferreira Batista

Laboratório de Biometria Ecológica
Centro de Métodos Quantitativos
Departamento de Ciências Florestais
Escola Superior de Agricultura "Luiz de Queiroz"
UNIVERSIDADE DE SÃO PAULO