## Script: 09_2varEstadoParces.R                  /
## Sobre:  calculo de variables a nivel agregado  /
##  por unidad muestral. Calculo de N, G, Dg y V. /
##
## Profesor: Christian Salas Eljatib
## E-mail: cseljatib AT gmail DOT com
## Web: https://eljatib.com
#====================================================
#Calculo de variables de estado de rodal, en donde  #
#  no se tienen alturas medidas para todos los 
#  arboles dentro de cada unidad muestral.  
#Se considera el ajuste de un modelo de altura-diametro
#  y su aplicacion solo para los arboles no-muestra.
#----------------------------------------------------# 
library(datana)
data(radiatapl)
df <- radiatapl
head(df)
df

unique(df$plot)
tapply(df$dbh, df$plot, length)

df$sup.parcela <- 150 #m2
df$fe <- 10000/df$sup.parcela
head(df)


##@@@@@@@@@@@@@@@@@@@@
## 1. Variables agregadas
##@@@@@@@@@@@@@@@@@@@@
#densidad por hectarea y parcela
nha.ppar<-tapply(df$fe, df$plot, sum)
nha.ppar

#area basal por arbol
df$g <- (pi/40000)*df$dbh^2
head(df)

#area basal del arbol en la hectarea
df$garb.ha <- df$g * df$fe
head(df)

#area basal por hectarea y parcela
gha.ppar<-tapply(df$garb.ha, df$plot, sum)
gha.ppar

#diam.medio.cuadratico por parcela
dmc.ppar<-sqrt((gha.ppar/nha.ppar)*(40000/pi))
dmc.ppar



####ajuste modelo altura
##solo los arboles con medicion de altura
df.h <- df[!is.na(df$height),]
dim(df.h)

summary(df.h$height)
plot(height~dbh, data=df)
plot(height~dbh, data=df, ylim=c(0,23),
xlim=c(0,40))

mod1<-lm( I(log(height))~ I(1/dbh), data=df)
summary(mod1)
b0.hat<-coef(mod1)[1]
b1.hat<-coef(mod1)[2]

df.h$h.final <- df.h$height

##arboles sin medicion de altura
df.noh <- df[is.na(df$height),]
dim(df.noh)

df.noh$h.final <- exp(b0.hat+b1.hat*(1/df.noh$dbh))

df2<- rbind(df.h,df.noh)
head(df2)
tail(df2)



##volumen
## usaremos un factor de forma
## ff=v/v_r
## ff * v_r = v
ff <- 0.3
df2$v  <- df2$g*df2$h.final*ff
head(df)

#volumen del arbol en la hectarea
df2$varb.ha <- df2$v * df$fe

#volumen por hectarea y parcela
vha.ppar<-tapply(df2$varb.ha, df2$plot, sum)
vha.ppar


##variables de estado por parcela
##crearemos nuevos objetos solo por conveniencia
# en la preparacion de la data que contendra las variables por parcela.
nha<-nha.ppar;gha<-gha.ppar;dg<-dmc.ppar;vha<-vha.ppar
varest.ppar<-data.frame(nha,gha,dg,vha)
varest.ppar

write.csv(varest.ppar, file="varestadoPorParcela.csv",row.names = FALSE)

#@@@@@@@@@@@@@@@@@
#eso es todo estimad@ estudiante
#disfRuten!
#saludos
#C
#@@@@@@@@@@@@@@@@@