| [[http://cmq.esalq.usp.br|{{:publico:tree-of-life.gif?100|CMQ: Centro de Métodos Quantitativos}}]] ^ **__Centro de Métodos Quantitativos__** \\ Departamento de Ciências Florestais \\ Escola Superior de Agricultura "Luiz de Queiroz" \\ **UNIVERSIDADE DE SÃO PAULO** | [[publico:tutoriais:uso-r:start|{{:publico:tutoriais:r-relampago:r-logo.png?90}}]] |
^ |
Curso Básico de R
para Mensuração Florestal:
Inventário Florestal
^ |
| ^ Inventário com Delineamento Simples | |
\\
====== Objetivo ======
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:
- organizar os dados de forma que sua estrutura seja compatível com o processamento no R;
- realizar as operações de agregação dos dados do nível de árvore para o nível de parcela;
- obter estimativas a nível de povoamento (talhao, propriedade, etc.) num delineamento amostral simples.
====== Estrutura dos Dados de Inventário ======
Ilustraremos a estrutura dos dados e processamento de inventário com o arquivo {{.:exemplo-inventario.csv|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.
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
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
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).
===== Estrutura em Data Frames Relacionados =====
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.
===== Observações Inexistentes =====
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 ) )
=== Exercício ===
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?
\\
====== Passos para o Processamento do Inventário ======
===== (1) Predição das Alturas Não Medidas =====
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)
===== (2) Predição do Volume das Árvores =====
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)
===== (3) Agregação no Nível de Parcela =====
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)
===== (4) Obtenção das Estimativas =====
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
\\
\\
====== Autor ======
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