##########################################################################################
##########
########## 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: 2026/08
##########
##########################################################################################
# Leitura dos Dados
egr1 <- read.csv("dados/Inventario-E_grandis-Rotacao_1.csv")

 #Inventario-E_grandis-Rotacao_2.csv
 #Inventario-E_saligna-Rotacao_1.csv
 #Inventario-E_saligna-Rotacao_2.csv

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

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

  # Altura pela Idade na ESCALA DE AJUSTE
scatter.smooth(1/egr1$idade, log(egr1$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/egr1$idade, log(egr1$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(egr1$idade, egr1$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(fdata=egr1)
 {
	plot(mhdom ~ idade, fdata, 
     	xlim=c(3,12.5),
     	ylim=c(14,30),
     	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(egr1)
coef.curvag
egr1$S.cg <- round(egr1$mhdom/exp(coef.curvag[2]*(1/egr1$idade - 1/5)),0)

graf.pred.cg <- function(fdata=egr1)
 {
	plot(fdata$idade, fdata$mhdom, 
	     xlab="Idade (anos)", ylab="",
	     cex.axis=1.2, cex.lab=1.4)
	mtext("Altura Med. Dominantes (m)", 2, 2.5, cex=1.4 )
	pred <- fdata$S.cg
	n <- length(pred)
	idade.base <- 5
	points(rep(idade.base,n), pred, pch=16, col="red", cex=1.8)
	abline(v=idade.base, lty=2, col="blue")
 }

graf.pred.cg() 

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

dados.eqdif <- function(fdata=egr1)
 {
   	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)
 }

egr1.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), 
		  egr1.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 <- egr1.par[ c(1, 150), ]
PARCELA.exemplo 

graf.projeta <- function(fdata=egr1)
 {
  scatter.smooth(fdata$idade, fdata$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")
 }

graf.projeta()


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

graf.eqdif <- function(fdata=egr1)
 {
	plot(mhdom ~ idade, egr1, 
     	   xlim=c(3,12.5), ylim=c(13,30),
     	   xlab="Idade (anos)", ylab="",
     	   main="Curvas de Sítio - Método da Equação da Diferença", 
	   cex.lab=1.4, cex.main=1.8)
     	mtext("Alt. Média das Dominantes (m)", 2, 2.5, cex.lab=1.4)
	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(egr1)
coef.eqdif
egr1$S.ed <- round(egr1$mhdom/exp(coef.eqdif*(1/egr1$idade - 1/5)),0)



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

graf.compar <- function(fdata=egr1)
{
	plot(mhdom ~ idade, fdata, 
     		xlim=c(3,10), ylim=c(13,30),
     		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
egr1

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

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

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

graf.S.pred <- function(fdata=egr1.sitio)
 {
        n <- dim(fdata)[1]
	plot((1:n)-0.125, fdata$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:n, labels=fdata$parcela, cex=1.3)
	abline(v=1:n, lty=9, col="darkgreen")
	points(1:n-0.125, fdata$S.cg.max,xaxt="n", pch=16, col="blue")
	segments(1:n-0.125, fdata$S.cg.min, 1:n-0.125, fdata$S.cg.max, 
		 col="blue")
	points(1:n+0.125, fdata$S.ed.min,xaxt="n", pch=16, col="red")
	points(1:n+0.125, fdata$S.ed.max,xaxt="n", pch=16, col="red")
	segments(1:n+0.125, fdata$S.ed.min, 1:n+0.125, fdata$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 de 
  #      * E. grandis, 2a. rotação 
  #      * E. saligna, 1a. rotação 
  #      * E. saligna, 2a. rotação 

  # Não se esqueça de subtrair "7" da idade dos dados de 2a.rotação


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