#07.06.2023 fchianucci@gmail.com

rm(list=ls())

if (!requireNamespace(c("tidyverse",'readxl','devtools'), quietly = TRUE)){
  install.packages(c("tidyverse",'readxl','devtools'))
}

library(tidyverse)
library(readxl)
library(devtools)

#upload the dataset
path <- 'https://zenodo.org/record/7773762/files/dataset_carpaneta_2023-03-27.xlsx'

p1f <- tempfile()
download.file(path, p1f, mode="wb")
carpaneta<-read_excel(path = p1f, sheet = 1)

#inspect the dataset:
head(carpaneta)
str(carpaneta)

#plot:
carpaneta |> 
  ggplot(aes(x=BX, y=BY, size=DBH, col=SP))+
  geom_point()+
  theme_bw()+
  facet_grid(~inv_year)



#tree structure analysis
carpaneta |> 
  filter(VIT==1)|> #only alive trees 
  group_by(CA, inv_year)|> 
  summarise(N=n(),
            QMD=sqrt(sum(DBH^2,na.rm=T)/N),#quadratic mean diameter
            MD=mean(DBH,na.rm=TRUE),#mean diameter
            G=pi/4*sum((DBH/100)^2,na.rm=T),#basal area
            MT=mean(TH,na.rm=TRUE), #mean tree height
              N_ha=mean(N*(70*140)/10000,na.rm=T),#tree density
            G_ha=mean(G*(70*140)/10000,na.rm=T))#basal area per ha


#install treespat for spatial analysis
devtools::install_gitlab('fchianucci/treespat')

library(treespat)

#calculate diameter dominance:
DDOM(carpaneta,
     .x = BX,
     .y = BY,
     .mark = DBH,
     xmax = 70, ymax = 140, max.k = 4, shape='square',
     .groups=c('CA','inv_year'))

#calculate diameter differentiation:
DIFF(carpaneta,
     .x = BX,
     .y = BY,
     .mark = DBH,
     xmax = 70, ymax = 140, max.k = 4, shape='square',
     .groups=c('CA','inv_year'))


#calculate mingling:
MING(carpaneta,
            .x = BX,
            .y = BY,
            .species = SP,
            xmax = 70, ymax = 140, max.k = 4,
            .groups=c('CA','inv_year'))

#calculate Stand Complexity Index (SCI):
SCI(carpaneta,
    .x = BX,
    .y = BY,
    .mark = DBH,
    .groups=c('CA','inv_year'))