diff --git a/01_Extract_data_Project_31.Rmd b/01_Extract_data_Project_31.Rmd
index de59e4fb2735f66450fd26b0cc9ca84514054ac2..ebfd3858d74135c5ca41c3c12be7f03074e31a1e 100755
--- a/01_Extract_data_Project_31.Rmd
+++ b/01_Extract_data_Project_31.Rmd
@@ -87,7 +87,7 @@ DT.xylem <- DT2 %>%
 
 ```
 
-Out of the `r length(species_list)` species in the sRoot list, `r sum(unique(DT2$species) %in% species_list)` species are present in sPlot, for a total of `r nrow(DT.xylem %>% filter(species %in% species_list))` plots., across `r length(plot.sel)` plots.
+Out of the `r length(species_list)` species in the sRoot list, `r sum(unique(DT2$species) %in% species_list)` species are present in sPlot, for a total of `r nrow(DT.xylem %>% filter(species %in% species_list))` records, across `r length(plot.sel)` plots.
 
 # 2 Extract woody species
 This is partial selection, as we don't have information on the growth form of all species in sPlot
@@ -95,8 +95,7 @@ This is partial selection, as we don't have information on the growth form of al
 #load list of woody species, as provided to me by Alexander Zizka, within sPlot project #21
 load("../Project_21/_input/evowood_species_list.rda")
 
-#Select all woody species, as well as those without GrowthForm information
-#extract relevant traits from TRY
+#Select all woody species and extract relevant traits from TRY
 woody_species_traits <- sPlot.traits %>%
   dplyr::select(species, GrowthForm, is.tree.or.tall.shrub, n,
                 starts_with("StemDens"),
@@ -111,17 +110,20 @@ woody_species_traits <- sPlot.traits %>%
           (species %in% synonyms$name_binomial) |
           grepl(pattern = "tree|shrub", x = GrowthForm) |
           is.tree.or.tall.shrub==T
-              )
+              ) %>% 
+  #counter proof - exclude species NOT herb
+  filter(GrowthForm != "herb" | is.na(GrowthForm))
 
 table(woody_species_traits$GrowthForm, exclude=NULL)
 # MEMO: some standardization needed in sPlot 3.0
-# Including the data from A.Zizka in the selection 
-# increases the number of selected species only marginally (from ~21k to ~22k)
+#
+# Using data from A.Zizka whhen selecting species
+# improves the selection only marginally (from ~21k to ~22k)
 ```
 
 ```{r, echo=F}
 knitr::kable(woody_species_traits %>% 
-               sample_n(20), caption="Example of gap-filled trait data from TRY (20 randomly selected species") %>%
+               sample_n(20), caption="Example of gap-filled trait data from TRY (20 randomly selected species)") %>%
     kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, position = "center")
 ```
 
@@ -146,13 +148,39 @@ DT.xylem <- DT.xylem %>%
 nrow(DT.xylem)
 ```
 
+Merge relative cover across vegetation layers, if needed, and normalize to 1 (=100%)
+```{r, cache=T}
+###combine cover values across layers
+combine.cover <- function(x){
+    while (length(x)>1){
+      x[2] <- x[1]+(100-x[1])*x[2]/100
+      x <- x[-1]
+    }
+    return(x)
+}
+
+DT.xylem <- DT.xylem %>%
+  dplyr::select(PlotObservationID, species,Layer, Relative.cover) %>%
+  # normalize relative cover to 100
+  left_join({.} %>% 
+              group_by(PlotObservationID, Layer) %>% 
+              summarize(Tot.Cover=sum(Relative.cover), .groups="drop"), 
+            by=c("PlotObservationID", "Layer")) %>% 
+  mutate(Relative.cover=Relative.cover/Tot.Cover) %>% 
+  group_by(PlotObservationID, species) %>%
+  summarize(Relative.cover=combine.cover(Relative.cover), .groups="drop") %>%
+  ungroup()
+
+nrow(DT.xylem)
+```
+
 
 # 3 Calculate CWMs and trait coverage 
 Calculate CWM and trait coverage for each trait and each plot. Select plots having more than 80% coverage for at least one trait.
 ```{r, cache=T,  cache.lazy=F}
 # Merge species data table with traits
 CWM.xylem0 <- DT.xylem %>%
-  as.tbl() %>%
+  as_tibble() %>%
   dplyr::select(PlotObservationID, species, Relative.cover) %>%
   left_join(xylem_data %>%
               dplyr::rename(species=Species) %>%
@@ -163,7 +191,8 @@ CWM.xylem0 <- DT.xylem %>%
 CWM.xylem1 <- CWM.xylem0 %>%
   group_by(PlotObservationID) %>%
   summarize_at(.vars= vars(P50:Ks),
-               .funs = list(~weighted.mean(., Relative.cover, na.rm=T))) %>%
+               .funs = list(~weighted.mean(., Relative.cover, na.rm=T)), 
+               .groups="drop") %>%
   dplyr::select(PlotObservationID, order(colnames(.))) %>%
   pivot_longer(-PlotObservationID, names_to="trait", values_to="trait.value")
 
@@ -239,8 +268,8 @@ CWM.xylem08 <- CWM.xylem %>%
 
 ```{r, echo=F}
 knitr::kable(CWM.xylem08 %>%
-               group_by(variable) %>%
-               summarize("num.plots"=n()), caption="Number of plots with >=.8 coverage per trait") %>%
+               group_by(trait) %>%
+               summarize("num.plots"=n(),.group="drop"), caption="Number of plots with >=.8 coverage per trait") %>%
     kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, position = "center")
 ```
 
@@ -282,11 +311,11 @@ ggworld
 
 Summarize data across data sets in sPlot, and create list of data custodians
 ```{r, warning=F, message=F}
-db.out <- fread("/data/sPlot/users/Francesco/_sPlot_Management/Consortium/Databases.out.csv") %>%
+db.out <- read_csv("/data/sPlot/users/Francesco/_sPlot_Management/Consortium/Databases.out.csv") %>%
   dplyr::select(`GIVD ID`, Custodian)
 data.origin <- header.xylem %>% 
    group_by(`GIVD ID`) %>% 
-   summarize(Num.plot=n()) %>%
+   summarize(Num.plot=n(), .groups="drop") %>%
    left_join(db.out)
 ```
 
@@ -307,7 +336,7 @@ soilclim.xylem <- soilclim %>%
 
 ```{r, echo=F}
 knitr::kable(soilclim.xylem %>% 
-               sample_n(20), caption="Example of climatic and soil variables for 20 randomly selected plots. All values represent the median in a circle centered on the plot coordinates, having a radius equal to the plot's location uncertainty (capped to 50 km for computing reasons)") %>%
+               sample_n(20), caption="Example of climatic and soil variables for 20 randomly selected plots. All values represent the mean in a circle centered on the plot coordinates, having a radius equal to the plot's location uncertainty (capped to 50 km for computing reasons). Sd is also reported.") %>%
     kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, position = "center")
 ```
 
@@ -354,7 +383,7 @@ P.ret.cat                  Phosphorous Retention - Categorical value, see [ISRIC
 \newline \newline
 
 
-# 4 Export & SessionInfo
+# 5 Export & SessionInfo
 
 ```{r, eval=T}
 save( woody_species_traits, DT.xylem, CWM.xylem, header.xylem,