# Retake the NMDS ----

# Load the data
download.file(
  url = "https://kamapu.gitlab.io/multivar/data/course-data.zip",
  destfile = "course-data.zip", method = "curl")
unzip("course-data.zip", overwrite = TRUE)
unlink("course-data.zip")

# Load Packages
library(vegan)
library(vegtable)

# Load data
wetlands <- readRDS("wetlands.rds")
names(wetlands)

# NMDS
nmds_ord <- metaMDS(wetlands$cross)
nmds_ord

plot(nmds_ord)

hist(wetlands$obs$cover_perc)
hist(wetlands$cross)

hist(log(1 + wetlands$cross))

plot(nmds_ord)

# Customizing the plot
plot(nmds_ord, type = "n")
abline(h = 0, lty = "dashed")
abline(v = 0, lty = "dashed")
points(scores(nmds_ord, display = "sites"), pch = 16,
       col = "grey")

head(wetlands$spp)
subset(wetlands$spp, grepl("papyrus", taxon_name))

sel_plots <- subset(wetlands$obs, taxon_code == "Cypepapy")
sel_plots <- unique(sel_plots$plot_id)
sel_plots <- paste(sel_plots)
sel_plots
class(sel_plots)

head(scores(nmds_ord, display = "sites"))
class (rownames(scores(nmds_ord, display = "sites")))

points(nmds_ord, pch = 16, col = "darkgreen", select = sel_plots)
points(nmds_ord, display = "species", pch = 3, col = "red",
       select = "Cypepapy")

head(wetlands$env)

scores_nmds <- scores(nmds_ord, display = "sites")
scores_nmds <- data.frame(
  plot_id = as.integer(rownames(scores_nmds)),
  scores_nmds)
head(scores_nmds)

wetlands$env <- merge(wetlands$env, scores_nmds)
head(wetlands$env)

pairs(wetlands$env[c("groundwater", "tmean", "prsum",
                     "NMDS1", "NMDS2")])

plot(nmds_ord)
ordisurf(nmds_ord ~ tmean, data = wetlands$env, add = TRUE)
ordisurf(nmds_ord ~ prsum, data = wetlands$env, add = TRUE,
         col = "blue")

nmds_ord

# Introduction to vegtable ----
rm(list = ls())

library(vegtable)

wetlands <- readRDS("wetlands-vegtable.rds")
wetlands

slotNames(wetlands)

wetlands@species
summary(wetlands@species, "papyrus")
indented_list(wetlands@species, "papyrus")
indented_list(wetlands@species, "Cyperaceae", synonyms = TRUE)

print_name(wetlands@species, 206)

## Functional traits / life forms -----
head(wetlands@species@taxonTraits)
head(wetlands@header)
head(wetlands@samples)

wetlands <- trait_proportion(
  trait = "life_form",
  object = wetlands,
  weigth = "cover_percentage"
)
head(wetlands@header)

boxplot(reed_plant_prop ~ community_type, data = wetlands@header,
        horizontal = TRUE)