##########################################################################################
##########
########## Curvas-de-Sitio-script-R
##########
##########                 Ajuste de Curvas de Sítio pelos métodos
##########
##########           * Método da Curva Guia
##########           * Método da Equação da Diferença
##########
########## Batista: 2025/08
##########
##########################################################################################
# Leitura dos Dados
euca <- read.table("inventario-continuo.csv",as.is=TRUE,header=TRUE)

  # Separação dos dados por Rotacao
euca1 <- subset(euca, subset = rotacao == 1)
euca2 <- subset(euca, subset = rotacao == 2)


# Métdo da Curva Guia
## Gráficos Exploratórios

  # Altura pela Idade na ESCALA ORIGINAL
scatter.smooth(euca1$idade, euca1$mhdom, lpars=list(lwd=2, col="red"))

  # Altura pela Idade na ESCALA DE AJUSTE
scatter.smooth(1/euca1$idade, log(euca1$mhdom), lpars=list(lwd=2, col="red"))

## Ajuste da Curva Guia
schum.curvag <- lm( log(mhdom) ~ I(1/idade) , euca1)
plot(schum.curvag)
summary(schum.curvag)
coef.curvag <- coef(schum.curvag)
coef.curvag 

## Gráficos das Curvas de Sítio
  # Escala de Ajuste
scatter.smooth(1/euca1$idade, log(euca1$mhdom), lpars=list(lwd=2, col="red"),
	main="Método da Curva Guia", cex.main=2)
curve( coef.curvag[1] + coef.curvag[2]*x, 0.1, 0.35, add=TRUE, lwd=2, col="blue")
abline(v=1/5, lty=9, lwd=2, col="blue")

  # Escala Original
scatter.smooth(euca1$idade, euca1$mhdom, lpars=list(lwd=2, col="red"),
	main="Método da Curva Guia", cex.main=2)
curve( exp(coef.curvag[1] + coef.curvag[2]*(1/x)), 3, 9, add=TRUE, lwd=2, col="blue")
abline(v=5, lty=9, lwd=2, col="blue")

  # Curvas de Sítio
graf.curvag <- function()
 {
	plot(mhdom ~ idade, euca1, 
     	xlim=c(3,12.5),
     	xlab="Idade (anos)", ylab="Alt. Média das Dominantes (m)",
     	main="Curvas de Sítio - Método da Curva Guia", cex.main=1.8)
	abline(v=5, lwd=2, lty=2, col="blue")
	S <- seq(18,24,by=2)
	ypos <- seq(22,30,by=2.3)
	for(i in 1:length(S))
 	{
   		curve( S[i]*exp(coef.curvag[2]*(1/x - 1/5)), 3, 10, add=TRUE, lwd=2, col="red")
   		text(10.2, ypos[i], paste("Índice de Sítio", S[i], "m"), adj=0, col="blue", cex=1.3)
 }
 }

graf.curvag() 

## Predição do Índice de Sítio

head(euca1)
coef.curvag
euca1$S.cg <- round(euca1$mhdom/exp(coef.curvag[2]*(1/euca1$idade - 1/5)),0)

plot(euca1$idade, euca1$mhdom )
points(rep(5,126), euca1$S.cg, pch=16, col="red")


subset(euca1, subset = parcela == 1101 )

# Métdo da Equação da Diferença
## Gerando os Dados para Equação da Diferença

dados.eqdif <- function(fdata=euca1)
 {
   	PARCELA <- unique(fdata$parcela)
	out <- data.frame(parcela=NA, idade_1=NA, mhdom_1=NA, idade_2=NA, mhdom_2=NA)
	for(k in 1:length(PARCELA))
	{
	  	tmp <- subset( fdata, subset = parcela == PARCELA[k],
	       			select = c("parcela","idade","mhdom")	)
		tmp <- tmp[ order(tmp$idade), ]
		n <- dim(tmp)[1]
		tmp.par <- cbind(tmp[1,] , tmp[2,-1])
		for(i in 1:(n-1))
		{
		     for(j in (i+1):n)
		     {
				tmp2 <- cbind(tmp[i,] , tmp[j,-1])
				tmp.par <- rbind(tmp.par, tmp2)
		     }
		}
		names(tmp.par) <- c("parcela","idade_1","mhdom_1","idade_2","mhdom_2")
		out <- rbind(out, tmp.par[-1,])	
	}
	out <- out[-1,]
	dimnames(out)[[1]] <- 1:dim(out)[1]
	return(out)
 }

euca1.par <- dados.eqdif()

## Ajuste do Modelo pela Equação da Diferença

schum.eqdif <- lm( I(log(mhdom_2) - log(mhdom_1) ) ~ -1 + I(1/idade_2  - 1/idade_1), euca1.par)
plot(schum.eqdif)
summary(schum.eqdif)
coef.eqdif <- coef(schum.eqdif)
coef.eqdif 


## Gráficos das Curvas de Sítio

  # Índice de Sítio de Duas Parcelas

PARCELA.exemplo <- euca1.par[ c(1, 210), ]
PARCELA.exemplo 

scatter.smooth(euca1$idade, euca1$mhdom, lpars=list(lwd=2, col="red"),
		main="Método da Equação da Diferença", cex.main=2)
abline(v=5, lty=2, lwd=2, col="blue")
points(mhdom_1 ~ idade_1, PARCELA.exemplo, pch=16, cex=2, col="blue")
curve(PARCELA.exemplo$mhdom_1[1] * exp(coef.eqdif*(1/x - 1/PARCELA.exemplo$idade_1[1])), 
      PARCELA.exemplo$idade_1[1], 5, add=TRUE, lwd=2, col="blue")
curve(PARCELA.exemplo$mhdom_1[2] * exp(coef.eqdif*(1/x - 1/PARCELA.exemplo$idade_1[2])), 
      PARCELA.exemplo$idade_1[2], 5, add=TRUE, lwd=2, col="blue")

  # Curvas de Sítio pela Equação da Diferença

graf.eqdif <- function()
 {
	plot(mhdom ~ idade, euca1, 
     	xlim=c(3,12.5),
     	xlab="Idade (anos)", ylab="Alt. Média das Dominantes (m)",
     	main="Curvas de Sítio - Método da Equação da Diferença", cex.main=1.8)
	abline(v=5, lwd=2, lty=2, col="blue")
	S <- seq(18,24,by=2)
	ypos <- seq(22.4,30,by=2.5)
	for(i in 1:length(S))
 	{
   		curve( S[i]*exp(coef.eqdif*(1/x - 1/5)), 3, 10, add=TRUE, lwd=2, col="red")
   		text(10.2, ypos[i], paste("Índice de Sítio", S[i], "m"), adj=0, col="blue", cex=1.3)
       }
 }

graf.eqdif()

## Predição do Índice de Sítio

head(euca1)
coef.eqdif
euca1$S.ed <- round(euca1$mhdom/exp(coef.eqdif*(1/euca1$idade - 1/5)),0)

subset(euca1, subset = parcela == 1101 )


# Comparação dos Métodos
## Comparação das Curvas de Sítio

graf.compar <- function()
{
	plot(mhdom ~ idade, euca1, 
     		xlim=c(3,10),
     		xlab="Idade (anos)", ylab="Alt. Média das Dominantes (m)",
     		main="Curvas de Sítio", 
		cex.main=1.8, cex.lab=1.5)
	abline(v=5, lwd=2, lty=2, col="blue")
	S <- seq(18,24,by=2)
	ypos <- seq(22.4,30,by=2.5)
	for(i in 1:length(S))
 	{
   		curve( S[i]*exp(coef.eqdif*(1/x - 1/5)), 3, 10, 
			 add=TRUE, lwd=2, col="red")
   		curve( S[i]*exp(coef.curvag[2]*(1/x - 1/5)), 3, 10, 
			 add=TRUE, lwd=2, col="blue", lty=9)
        }
	legend(7, 18, legend=c("Eq. Diferença","Curva Guia"),
                      lwd=c(2,2), lty=c(1,9), col=c("red","blue"), 
		      title = "Métodos de Construção",
		      cex = 1.5	, title.cex = 1.5	)
}

graf.compar()

## Comparações das Predições do Índice de Sítio
euca1

tmp.bas <- aggregate( cbind(idade,mhdom) ~ parcela, euca1, function(x) x[1] )
tmp.min <- aggregate( cbind(S.cg,S.ed) ~ parcela, euca1, min )
names(tmp.min)[2:3] <- c("S.cg.min","S.ed.min")
tmp.max <- aggregate( cbind(S.cg,S.ed) ~ parcela, euca1, max )
names(tmp.max)[2:3] <- c("S.cg.max","S.ed.max")

euca1.sitio <- merge(tmp.bas, tmp.min)
euca1.sitio <- merge(euca1.sitio, tmp.max)

euca1.sitio <- euca1.sitio[ order(euca1.sitio$S.cg.min), ]
range(euca1.sitio$S.cg.max)
range(euca1.sitio$S.ed.max)

graf.S.pred <- function()
 {
	plot(1:21-0.125, euca1.sitio$S.cg.min,xaxt="n", ylim=c(13,30) , 
	     col="blue", pch=16,
	     xlab="Parcelas", ylab="", cex.lab=1.5,
	     main="Comparação das Predições do Índice de Sítio", cex.main=2,
	     cex.axis=1.3)
	mtext( "Índice de Sítio", 2, 2, cex=1.5)
	axis(1, at=1:21, labels=euca1.sitio$parcela, cex=1.3)
	abline(v=1:21, lty=9, col="darkgreen")
	points(1:21-0.125, euca1.sitio$S.cg.max,xaxt="n", pch=16, col="blue")
	segments(1:21-0.125, euca1.sitio$S.cg.min, 1:21-0.125, euca1.sitio$S.cg.max, 
		 col="blue")
	points(1:21+0.125, euca1.sitio$S.ed.min,xaxt="n", pch=16, col="red")
	points(1:21+0.125, euca1.sitio$S.ed.max,xaxt="n", pch=16, col="red")
	segments(1:21+0.125, euca1.sitio$S.ed.min, 1:21+0.125, euca1.sitio$S.ed.max, 
		 col="red")
	legend(2, 29, legend=c("Método da Curva Guia","Método da Equação da Diferença"),
	       col=c("blue","red"), lwd=c(2,2), pch=16, cex=1.5, bg="white")
 }

graf.S.pred()


# EXERCÍCIO

  # Repita o código acima com os dados da tabela 'euca2'
  # Antes porém faça a correção da idade:
euca2$idade <- euca2$idade - 7

scatter.smooth(euca2$idade, euca2$mhdom, lpars=list(col="red",lwd=2))

##########################################################################################
##########################################################################################
##########################################################################################
save.image()