Load required libraries

library(plyr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)

## -- Attaching packages --------------------- tidyverse 1.2.1 --

## v ggplot2 3.0.0     v readr   1.1.1
## v tibble  1.4.2     v purrr   0.2.5
## v tidyr   0.8.1     v stringr 1.3.1
## v ggplot2 3.0.0     v forcats 0.3.0

## -- Conflicts ------------------------ tidyverse_conflicts() --
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()

library(TapeR)

## Warning: package 'TapeR' was built under R version 3.5.2

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: splines

## Loading required package: pracma

## Warning: package 'pracma' was built under R version 3.5.2

## 
## Attaching package: 'pracma'

## The following object is masked from 'package:purrr':
## 
##     cross

First test using TapeR pachage

getwd()

## [1] "C:/Users/ro/Documents/WORK/RICERCA/GitLab/Forest-Management-Tools_pinaster-and-radiata-pine_Arzana/TaperFunctions"

load("TaperFunctions.RData")
aaa_test <- ProfiliFustiPrincipali %>%
  right_join(select(FustiCampioneEtAdF, KeyID, h_ipso, ID_fustoCampione))

## Joining, by = "KeyID"

aaa_test$ID_fustoCampione <- as.numeric(aaa_test$ID_fustoCampione)
aaa_test$d_sez <- as.numeric(aaa_test$d_sez)

aaa_test <- aaa_test %>%
  select(ID_fustoCampione, d_sez, distSuolo, h_ipso)

plot(x = aaa_test$distSuolo/aaa_test$h_ipso, y = aaa_test$d_sez, pch = as.character(aaa_test$ID_fustoCampione), 
       xlab="Relative height (m)", ylab="Diameter (cm)")

Kublin <- TapeR_FIT_LME.f(Id = aaa_test$ID_fustoCampione, 
                x = aaa_test$distSuolo/aaa_test$h_ipso,
                y = aaa_test$d_sez,
                knt_x = c(0.0, 0.1, 0.75, 1.0),
                ord_x = 4,
                knt_z = c(0.0, 0.1, 1.0),
                ord_z = 4)

Testing taper function from paper ID10

# Testing functions
parameters <- tribble(
  ~param,      ~fnct,    ~value,
  "a0",        "FMOLS",  0.989, 
  "a1",        "FMOLS",  0.963, 
  "a2",        "FMOLS",  0.0458,
  "b1",        "FMOLS",  0.367, 
  "b2",        "FMOLS", -0.335, 
  "b3",        "FMOLS",  0.519, 
  "b4",        "FMOLS",  0.847, 
  "b5",        "FMOLS",  0.0178,
  "b6",        "FMOLS", -0.0265,
  "a0",        "MM3",    1.05,  
  "a1",        "MM3",    0.943, 
  "a2",        "MM3",    0.0473,
  "b1",        "MM3",    0.362, 
  "b2",        "MM3",   -0.691, 
  "b3",        "MM3",    0.585, 
  "b4",        "MM3",    1.13,  
  "b5",        "MM3",    0.0227,
  "b6",        "MM3",    0.0581
)

d_sez <- function(h_sez, d_130, h_ipso, FN = "FMOLS") {
  FN <- unique(FN)
  stopifnot(length(FN)==1, FN %in% parameters$fnct)
  fpar <- parameters %>% filter(fnct == FN)
  for(i in 1:nrow(fpar)) assign(fpar$param[i], fpar$value[i])
  q <- h_sez/h_ipso
  w <- 1-q^(1/3)
  x <- w/(1-(1.3/h_ipso)^(1/3))
  d <- a0*d_130^a1*h_ipso^a2*x^
      (b1*q^4+
       b2*(1/exp(d_130/h_ipso))+
       b3*x^0.1+
       b4*(1/d_130)+
       b5*h_ipso^w+
       b6*x)
}

npp <- 6 # Number of Profiles per Plot

aaa_test_ID10 <- ProfiliFustiPrincipali %>%
  right_join(select(FustiCampioneEtAdF, KeyID, h_ipso, d_130, ID_fustoCampione)) %>%
  filter(ID_fustoCampione <= npp*2) %>% 
  arrange(ID_fustoCampione)

## Joining, by = "KeyID"

output <- aaa_test_ID10 %>%
  mutate(FN = "obs") %>% 
  bind_rows(mutate(aaa_test_ID10, FN = "FMOLS", 
                   d_sez = d_sez(distSuolo, d_130,
                                 h_ipso, FN))) %>%
  bind_rows(mutate(aaa_test_ID10, FN = "MM3", 
                   d_sez = d_sez(distSuolo, d_130,
                                 h_ipso, FN)))
for(p in seq(1,max(output$ID_fustoCampione), npp)){
  plt <- output %>%
    filter(between(ID_fustoCampione, p, p+npp-1)) %>% 
      ggplot(aes(distSuolo, d_sez)) +
        geom_point(aes(1.3, d_130)) +
        geom_line(aes(linetype = FN)) +
        scale_linetype_manual(values=
              c("twodash", "dotted", "solid")) +
        facet_wrap(vars(ID_fustoCampione), 
                   ncol = 3, scales = "free")
  print(plt)
}