%\VignetteIndexEntry{TapeR: Flexible stem taper modelling} \documentclass[english,11pt]{article} %\documentclass[a4paper,11pt]{article} \usepackage[a4paper]{geometry} \geometry{verbose,tmargin=2cm,bmargin=2cm,lmargin=2cm,rmargin=2cm,footskip=1cm} \usepackage{Sweave} \usepackage{color} \usepackage{graphicx} \usepackage{booktabs} \usepackage[authoryear]{natbib} \usepackage{url} \usepackage[utf8]{inputenc} % \setlength{\oddsidemargin}{1cm} % \setlength{\textwidth}{15cm} % \setlength{\topmargin}{-2.2cm} % \setlength{\textheight}{22cm} % \usepackage{amssymb} % \usepackage{amsmath} %\usepackage{hyperref} \usepackage[linkcolor=blue, bookmarks=true, citecolor=blue, colorlinks=true, linktocpage, a4paper]{hyperref} \title{TapeR: Flexible stem taper modelling} \author{Edgar Kublin\thanks{Forest Research Institute of Baden-W{\"u}rttemberg, Freiburg, Germany}, Johannes Breidenbach\thanks{Norwegian Institute of Bioeconomy Research, 1431 {\AA}s, Norway, job@nibio.no}} %%\date{} %\address{Norwegian Institute for Forest and Landscape} \begin{document} \SweaveOpts{concordance=TRUE} \maketitle %\tableofcontents \section{Introduction} The \texttt{TapeR} package implements methods described in \citep{kublin2013}. If \texttt{R} is running, the \texttt{TapeR} package can be installed by typing <>= install.packages("TapeR") @ into the console\footnote{The character "\texttt{>}" is not part of the command. A working Internet connection is required.}. The command <<>>= library(TapeR) @ loads the package into the current workspace. We can get an overview of the packages' contents by typing <>= ?TapeR @ \section{Using the TapeR functions} \subsection{Fitting a taper model} A taper model can be fit to measured pairs of diameter and corresponding heights. Typically, those data are recorded for regular stem sections but the data can also be availabe at irregular distances. Important is for the method implemented in this packages is that the stem height is available too. We will load some example data and fit a taper model. Note that in real applications, the dataset must be considerably bigger in order to provide reliable taper models that are valid over larger spatial extents. For example, the data set \texttt{SK.par.lme} provides the taper model parameters based on 338 Norway spruce trees in Germany. <>= #load example data data(DxHx.df) #prepare the data (could be defined in the function directly) Id = DxHx.df[,"Id"] x = DxHx.df[,"Hx"]/DxHx.df[,"Ht"]#calculate relative heights y = DxHx.df[,"Dx"] #plot the example data plot(x,y,pch=as.character(Id),xlab="Relative height (m)", ylab="Diameter (cm)") #define the relative knot positions and order of splines knt_x = c(0.0, 0.1, 0.75, 1.0) # B-Spline knots: fix effects ord_x = 4 # ord = order (4 = cubic) knt_z = c(0.0, 0.1, 1.0); ord_z = 4 # B-Spline knots: rnd effects #fit the model taper.model <- TapeR_FIT_LME.f(Id, x, y, knt_x, ord_x, knt_z, ord_z, IdKOVb = "pdSymm") #save model parameters for documentation or dissimination #parameters can be load()-ed and used to predict the taper #or volume using one or several measured dbh spruce.taper.pars <- taper.model$par.lme #save(spruce.taper.pars, file="spruce.taper.pars.rdata")##uncomment to save @ \section{Example} This is from the R-help of the exported TapeR-functions: Let's take a look on the taper curve of the first tree in the example data set if it is only calibrated using the diameter measurement in 2m height. <>= #example data data(DxHx.df) #taper curve parameters based on all measured trees data(SK.par.lme) #select data of first tree Idi <- (DxHx.df[,"Id"] == unique(DxHx.df$Id)[1]) tree1 <- DxHx.df[Idi,] ## Predict the taper curve based on the diameter measurement in 2 m ## height and known height tc.tree1 <- E_DHx_HmDm_HT.f(Hx=1:tree1$Ht[1], Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1], sHt = 0, par.lme = SK.par.lme) #plot the predicted taper curve plot(tc.tree1$Hx, tc.tree1$DHx, type="l") #lower CI lines(tc.tree1$Hx, tc.tree1$CI_Mean[,1], lty=2) #upper CI lines(tc.tree1$Hx, tc.tree1$CI_Mean[,3], lty=2) #lower prediction interval lines(tc.tree1$Hx, tc.tree1$CI_Pred[,1], lty=3) #upper prediction interval lines(tc.tree1$Hx, tc.tree1$CI_Pred[,3], lty=3) #add measured diameter points(tree1$Hx[3], tree1$Dx[3], pch=3, col=2) #add the observations points(tree1$Hx, tree1$Dx) ## Add the population average taper curve (without calibration) to the ## plot (not of high practical interest but good to know how to get it). Ht = tree1$Ht[1] Hx = tree1$Hx #get fixed-effects design matrix for the Splines X = TapeR:::XZ_BSPLINE.f(x=Hx/Ht, knt = SK.par.lme$knt_x, ord = SK.par.lme$ord_x) #Calculate population average taper curve DHx_PA = X %*% SK.par.lme$b_fix #add to plot lines(tree1$Hx, DHx_PA, lwd=2, lty=4) @ Let's see how intervals change, if some uncertainty in height is assumed. <>= ## Let's see how CI's change if there's some uncertainty in the height ## measurement tc.tree1 <- E_DHx_HmDm_HT.f(Hx=1:tree1$Ht[1], Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1], sHt = 1, par.lme = SK.par.lme) #plot the predicted taper curve plot(tc.tree1$Hx, tc.tree1$DHx, type="l", xlab="Height (m)", ylab="Diameter (cm)") #lower CI lines(tc.tree1$Hx, tc.tree1$CI_Mean[,1], lty=2) #upper CI lines(tc.tree1$Hx, tc.tree1$CI_Mean[,3], lty=2) #lower prediction interval lines(tc.tree1$Hx, tc.tree1$CI_Pred[,1], lty=3) #upper prediction interval lines(tc.tree1$Hx, tc.tree1$CI_Pred[,3], lty=3) #add measured diameter points(tree1$Hx[3], tree1$Dx[3], pch=3, col=2) #add the observations points(tree1$Hx, tree1$Dx) @ Additionally, exact confidence intervals can be calculated: <>= plot(tc.tree1$Hx, tc.tree1$DHx, type="l", xlab="Height (m)", ylab="Diameter (cm)") ## Calculate "exact" CIs. Careful: This takes a while! #library(pracma)# for numerical integration with gaussLegendre() tc.tree1.exact <- E_DHx_HmDm_HT_CIdHt.f(Hx=1:tree1$Ht[1], Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1], sHt = 1, par.lme = SK.par.lme) #add exact confidence intervals to approximate intervals above - fits #quite well lines(tc.tree1.exact[,1], tc.tree1.exact[,2], lty=2,col=2) lines(tc.tree1.exact[,1], tc.tree1.exact[,4], lty=2,col=2) ## Calculate the height given a certain diameter threshold, say 8.5 cm ht.tree1.d8.5 <- E_HDx_HmDm_HT.f (Dx=8.5, Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1], sHt = 1, par.lme = SK.par.lme) #add to above created plot abline(v=ht.tree1.d8.5) ## Calculate the timber volume #for the whole stem E_VOL_AB_HmDm_HT.f(Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1] , sHt = 1, par.lme = SK.par.lme) #Calculate the timber volume for a selected section given a height #threshold (0.3 - 5 m) E_VOL_AB_HmDm_HT.f(Hm=tree1$Hx[3], Dm=tree1$Dx[3], mHt = tree1$Ht[1] , sHt = 1, par.lme = SK.par.lme, A=0.3, B=5, iDH = "H") #Calculate the timber volume for a selected section given a diameter #threshold (30 cm - 15 cm) (negative value if A