knitr::opts_chunk$set(echo = TRUE,comment="")
date()
## [1] "Thu Oct 29 10:14:01 2020"

Objectivo del reporte

Explicar como un modelo de ahusamiento puede ser ocupado como un simulador de trozado.

Breve introduccion

Un modelo de ahusamiento es una función matemática que describe el perfil fustal, es decir, el diametro del fuste (generalmente sin corteza) a una altura \(l\) desde la base del árbol (\(d_l\)) en función de características del árbol (e.g., su diámetro a la altura del pecho \(d\) y altura total \(h\)) y la altura fustal \(h_l\).

Al describir el perfil del fuste, es posible mediante integración numérica, calcular el volumen del árbol entre cualquier altura a definir. Por ejemplo, si se integra entre la base y la altura total, se calcularía el volumen total del árbol. No obstante la mayor ventaja de un modelo de ahusamiento es la posibilidad de calcular el volumen comercial de este, que se define de acuerdo a diámetros mínimos de trozos y/o sus largos.

Siendo mas específicos, el modelo de ahusamiento permite:

  1. Predecir el diámetro del fuste a una altura \(l\)

  2. Calcular el volumen entre dos alturas fustales (la diferencia de estas corresponde al largo de la sección) a partir de la predicción del diámetro fustal a dichas alturas.

  3. Estimar la altura fustal para cualquier diámetro fustal de interes o \(diu\) (diámetro a un índice de utilización), mediante la iteración de la función de ahusamiento para un valor fijo de \(d_l\) (lado izquierdo de la función) y los posibles valores de \(h_l\).

A continuación me centrare en explicar fundamentalmente como mediante (3) es posible calcular volúmenes para diferentes tipos de productos.

Un modelo de ejemplo

Uno de los modelos más conocidos en estudios de ahusamiento han sido llevados a cabo por el prof. Kozak de la University of British Columbia (Vancouver, Canada).

Por ejemplo el modelo de Kozak (1988) es el siguiente

\[d_{l}= \alpha_0 d^{\alpha_1}\alpha_2^{d} \left(X \right)^{\beta_1 Z^{2}+\beta_2 \ln(Z+0.001) +\beta_3 \sqrt{Z}+\beta_4 \exp{Z}+\beta_5 \frac{d}{h}} \tag{1}\] donde \(X=\left(1-\sqrt{h_l/h}\right)\left(1-\sqrt{p}\right)\); \(d\) es el diametro a la altura del pecho o dap; \(h\) es la altura total; \(h_l\) es la altura fustal a los \(l\) metros desde la base; \(d_l\) es el diametro fustal (sin corteza) a los \(l\) metros desde la base; \(p=0.25\); y \(Z=\frac{h_l}{h}\)

Este modelo se puede escribir en R como sigue

kozak88 <- function(dap, h, hi) {
 p <- 0.25
   dib <- 1.02453 * dap^0.88809 * 1.00035^dap
   X <- (1.0 - sqrt(hi / h)) / (1.0 - sqrt(p))
 Z=hi/h
 a=0.95086*Z*Z;
   b = -0.18090 * log(Z + 0.001);
   c = 0.61407 * sqrt(Z) +  -0.35105 * exp(Z);
   d = 0.05686 * (dap / h);
 retval<-(dib*X^(a+b+c+d))
 retval[h < hi] <- 0
 retval }

veamos como se ve el ahusamiento de un arbol de ejemplo con un \(d=45\) cm y \(h=27\) m.

dap.arb <- 45
h.arb <- 27
hi.vector <- 0:h.arb
di.hat <- kozak88(dap=dap.arb,h=h.arb,hi=hi.vector)

plot(di.hat~hi.vector, col="blue", type="l",las=1, ylab="Diametro fustal (cm)", xlab="Altura fustal (m)")

Altura comercial

Es clave poder encontrar la altura a la cual se encuentra un diametro de indice de utilizacion determinado, es decir, encontrar \(h_{iu}\) dado un valor de \(d_{iu}\). Para esto es necesario iterar en la ecuacion de ahusamiento, ya que debemos fijar el valor de la variable respuesta primero y luego encontrar los valores de \(h_l\) que logran ese valor.

Para esto escribamos la siguiente funcion

hiu.fx <- function(hi, dap, h, diu) {
   di <- kozak88(dap, h, hi) 
   dif <- di - diu
 dif
}

Por ejemplo, busquemos el \(h_l\) para un valor fijo de \(d_{l}=10\), es decir, busquemos el valor de \(h_{iu10}\).

diu.objetivo <- 10
hiu.hat <- uniroot(hiu.fx,
 c(0, h.arb),
 dap = dap.arb,
 h = h.arb,
 diu = diu.objetivo)

hiu.hat
$root
[1] 23.345

$f.root
[1] 1.5848e-06

$iter
[1] 5

$init.it
[1] NA

$estim.prec
[1] 6.1035e-05

Esto indica que a los 23.345 m de altura del fuste se encuentra el diametro fustal de 10 cm.

Ahora vamos a verificar esto mediante nuestro modelo de ahusamiento ya escrito, como sigue en R

diu.check <- kozak88(dap=dap.arb, h=h.arb, hiu.hat$root)
diu.check
[1] 10

el valor de \(d_{l}\) que predice la funcion de ahusamiento cuando \(h_{l}\) es igual a la altura del indice de utilizacion, es exactamente 10 cm, justamente lo que buscamos. Mediante esto verificamos que el valor que iteramos de \(h_{iu}\) esta correcto.

Productos comerciales

Existen diferentes productos madereros potenciales a obtener de un árbol. Las características de éstos dependenden de una serie de factores, dentro de los cuales los más importantes son el mercado y la especie. A continuación se utiliza a modo de ejemplo un cuadro que describe las dimensiones para distintos tipos de productos.

Producto/destino Dimensiones (diámetro de sección)
Pulpable \(2 < d_l< 5\)
Aserrable delgado \(5 < d_l< 12\)
Aserrable mediano \(12 < d_l< 18\)
Aserrable grande \(18 < d_l< 32\)
Debobinable \(d_l>32\)
trozo.breaks=c(2,5,12,18,32,999)

trozo.produ=c("pulpa","ase3","ase2","ase1","debo")

produ.names <- c("Pulpable","Aserrable #3","Aserrable #2",
                  "Aserrable #1", "Debobinable")

Ahora debemos almacenar los hiu y asi saber bien los quiebres entre productos

hiu.bks <- rep(0, length(trozo.breaks))

 for(i in 1:length(trozo.breaks)) {
    dap <- dap.arb; h <- h.arb; 
    if(trozo.breaks[i] <= kozak88(dap, h, 0.0)) {
      hiu <- uniroot(hiu.fx,
                    c(0, h),
                    dap = dap,
                    h = h,
                    diu = trozo.breaks[i])
   hiu.bks[i] <- hiu$root
} else {   
      hiu.bks[i] <- 0.0
}
}

los resultados estan

trozo.breaks
[1]   2   5  12  18  32 999
hiu.bks
[1] 26.4974 25.4056 22.4512 19.4072  4.5883  0.0000

El vector resultante hiu.bks tiene las alturas fustales a las cuales el fuste no puede ser comercializado segun las dimensiones definidas en trozo.break.

Por ejemplo, el trozo más largo que se puede obtener con un diámetro mínimo de 2 cm es de 26.49745. De la misma manera, el trozo más largo a obtener con un diámetro mínimo de 5 cm es 25.4056 m. Por lo tanto, el largo del trozo con diámetro mínimo de 2 cm y máximo de 5 cm es de \(26.49745-25.4056=1.0989m\).

Calculando el volumen de los trozos

El proceso es entonces: (a) determinar las alturas comerciales dado el set de diametros limites (trozo.breaks ); (b) calcular el largo del fuste para cada trozo.produ; y (c) usar un loop para calcular el volumen de cada trozo dentro del fuste, con el largo de cada trozo siendo la diferencia entre los hiu.

Primero creare una función simple que cubica trozos empleando la fórmula analítica de Smalian Prodan et al. (1997)

volsecc.fx <- function(d1, d2, l) {
 volsecc <- (pi/40000)/4*(d1^2 + d2^2) * l 
 volsecc
}
#la cual revisamos mediante
volsecc.fx(d1=c(10,5),d2=c(5,2),l=c(2,1))
[1] 0.00490874 0.00056941

Entonces ahora, solo nos falta calcular el largo de cada sección, pero tambien estimar (en base a la función de ahusamiento) el diámetro de tocón (sin corteza) para lo cual asumiremos una altura de tocón de 30 cm (0.3 m).

#el diametro de cada seccion
trozo.breaks
[1]   2   5  12  18  32 999
#la altura de cada seccion
hiu.bks
[1] 26.4974 25.4056 22.4512 19.4072  4.5883  0.0000
#el largo de cada seccion
largo.trozos1 <- diff(-hiu.bks)
largo.trozos1
[1]  1.0919  2.9543  3.0440 14.8189  4.5883
diam.min <- (trozo.breaks[1:(length(trozo.breaks)-1)])
diam.min
[1]  2  5 12 18 32
diam.max <- (trozo.breaks[2:(length(trozo.breaks)-1)])
diam.max
[1]  5 12 18 32
largo.trozo <-largo.trozos1
largo.trozo
[1]  1.0919  2.9543  3.0440 14.8189  4.5883
###calculo del diam. de tocon
hst <- 0.3
diu.tocon <- kozak88(45, 27, hst)
diu.tocon
[1] 43.444
diam.max <- c(diam.max,diu.tocon)
diam.max
[1]  5.000 12.000 18.000 32.000 43.444
#el largo del trozo mas grueso es corregido
largo.trozo[length(largo.trozo)]<- largo.trozo[length(largo.trozo)]-hst

#calculo de volumen de cada trozo, segun Smalian
vol.m3<- volsecc.fx(d1=diam.min,d2=diam.max,l=largo.trozo)

Entonces

df.trozos<-data.frame(diam.min,diam.max,largo.trozo, vol.m3)
df.trozos
  diam.min diam.max largo.trozo     vol.m3
1        2    5.000      1.0919 0.00062174
2        5   12.000      2.9543 0.00980341
3       12   18.000      3.0440 0.02797194
4       18   32.000     14.8189 0.39222650
5       32   43.444      4.2883 0.24513983

Bibliography

Kozak, A. 1988. “A Variable-Exponent Taper Equation.” Can. J. For. Res. 18: 1363–8.

Prodan, M., R. Peters, F. Cox, and P. Real. 1997. Mensura Forestal. San José, Costa Rica: Serie investigación y educación de desarrollo sostenible. Instituto interamericano de cooperación para la agricultura (IICA)/BMZ/GTZ.