Hy dar vision

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

src/000_set_environment.R

@@ -0,0 +1,80 @@
+# Set path ---------------------------------------------------------------------
+if(Sys.info()["sysname"] == "Windows"){
+ filepath_base = "D:/plygrnd/hySpecVisKili/"
+} else {
+ filepath_base = "/mnt/sd19006/data/users/tnauss/KI-Hyperspec/"
+}
+
+filepath_source = paste0(filepath_base, "hySpecVisKili/src/functions/001_functions.R")
+path_data = paste0(filepath_base, "/data/")
+path_biodiv = paste0(path_data, "/biodiv/")
+
+path_hyp_org = paste0(path_data, "/020_hypspec_org/")
+path_hyp_aio = paste0(path_data, "/025_hypspec_aio/")
+path_hyp_nrm = paste0(path_data, "/030_hypspec_nrm/")
+path_hyp_vegidcs = paste0(path_data, "/040_hypspec_vegidcs/")
+path_hyp_kmdc = paste0(path_data, "/050_hypspec_kmdc/")
+path_hyp_raoq = paste0(path_data, "/060_hypspec_raoq/")
+path_hyp_glcm = paste0(path_data, "/070_hypspec_glcm/")
+path_hyp_pred = paste0(path_data, "/090_hypspec_pred/")
+path_ldr_pred = paste0(path_data, "/095_lidar_pred/")
+path_comb_gpm_sr = paste0(path_data, "/100_comb_gpm_sr/")
+path_model_gpm_sr = paste0(path_data, "/110_model_gpm_sr/")
+path_compile_analysis_sr = paste0(path_data, "/120_compile_analysis_sr/")
+path_comb_gpm_sr_res = paste0(path_data, "/200_comb_gpm_sr_res/")
+path_model_gpm_sr_res = paste0(path_data, "/210_model_gpm_sr_res/")
+path_compile_analysis_sr_elev_res = paste0(path_data, "/220_compile_analysis_sr_elev_res/")
+path_model_gpm_sr_indp = paste0(path_data, "/300_model_gpm_sr_indp/")
+path_comb_gpm_sr_elev_res_indp = paste0(path_data, "/310_comb_gpm_sr_elev_res_indp/")
+path_model_gpm_sr_elev_res_indp = paste0(path_data, "/320_model_gpm_sr_elev_res_indp/")
+path_analysis_sr = paste0(path_data, "/500_analysis_sr/")
+path_analysis_sr_elev_res = paste0(path_data, "/510_analysis_sr_elev_res/")
+
+path_plots = paste0(path_data, "/plots/")
+path_rdata = paste0(path_data, "/rdata/")
+path_meta = paste0(path_data, "/meta/")
+path_temp = paste0(path_data, "/temp/")
+path_output = paste0(path_data, "/output/")
+path_vis = paste0(path_data, "/vis/")
+
+
+# Set libraries ----------------------------------------------------------------
+library(biodivTools) # devtools::install_github("environmentalinformatics-marburg/biodivTools")
+library(CAST)
+library(corrplot)
+library(doParallel)
+library(grid)
+library(gridExtra)
+library(gpm) # devtools::install_github("environmentalinformatics-marburg/gpm")
+library(ggplot2)
+# library(ggbiplot)
+library(hsdar)
+# library(lavaan)
+# library(rPointDB)
+library(rgeos)
+library(ggplot2)
+library(mapview)
+# library(metTools) # devtools::install_github("environmentalinformatics-marburg/metTools")
+library(raster)
+library(RStoolbox)
+library(reshape2)
+library(rgdal)
+# library(satellite)
+library(satelliteTools) # devtools::install_github("environmentalinformatics-marburg/satelliteTools")
+# library(semPlot)
+library(sp)
+library(spacetime)
+library(vegan)
+# library(yaml)
+
+
+# Other settings ---------------------------------------------------------------
+source(filepath_source)
+
+rasterOptions(tmpdir = path_temp)
+
+saga_cmd = "C:/OSGeo4W64/apps/saga-ltr/saga_cmd.exe"
+# initOTB("C:/OSGeo4W64/bin/")
+# initOTB("C:/OSGeo4W64/OTB-6.2.0-Win64/bin/")
+
+

src/000_set_environment_linux.R

@@ -0,0 +1,77 @@
+# Set path ---------------------------------------------------------------------
+if(Sys.info()["sysname"] == "Windows"){
+ filepath_base = "D:/plygrnd/hySpecVisKili/"
+} else {
+ filepath_base = "/mnt/sd19006/data/users/tnauss/KI-Hyperspec/"
+}
+
+filepath_source = paste0(filepath_base, "hySpecVisKili/src/functions/001_functions.R")
+path_data = paste0(filepath_base, "/data/")
+path_biodiv = paste0(path_data, "/biodiv/")
+
+path_hyp_org = paste0(path_data, "/020_hypspec_org/")
+path_hyp_aio = paste0(path_data, "/025_hypspec_aio/")
+path_hyp_nrm = paste0(path_data, "/030_hypspec_nrm/")
+path_hyp_vegidcs = paste0(path_data, "/040_hypspec_vegidcs/")
+path_hyp_kmdc = paste0(path_data, "/050_hypspec_kmdc/")
+path_hyp_raoq = paste0(path_data, "/060_hypspec_raoq/")
+path_hyp_glcm = paste0(path_data, "/070_hypspec_glcm/")
+path_hyp_pred = paste0(path_data, "/090_hypspec_pred/")
+path_ldr_pred = paste0(path_data, "/095_lidar_pred/")
+path_comb_gpm_sr = paste0(path_data, "/100_comb_gpm_sr/")
+path_model_gpm_sr = paste0(path_data, "/110_model_gpm_sr/")
+path_compile_analysis_sr = paste0(path_data, "/120_compile_analysis_sr/")
+path_comb_gpm_sr_res = paste0(path_data, "/200_comb_gpm_sr_res/")
+path_model_gpm_sr_res = paste0(path_data, "/210_model_gpm_sr_res/")
+path_compile_analysis_sr_elev_res = paste0(path_data, "/220_compile_analysis_sr_elev_res/")
+path_model_gpm_sr_indp = paste0(path_data, "/300_model_gpm_sr_indp/")
+path_comb_gpm_sr_elev_res_indp = paste0(path_data, "/310_comb_gpm_sr_elev_res_indp/")
+path_model_gpm_sr_elev_res_indp = paste0(path_data, "/320_model_gpm_sr_elev_res_indp/")
+path_analysis_sr = paste0(path_data, "/500_analysis_sr/")
+path_analysis_sr_elev_res = paste0(path_data, "/510_analysis_sr_elev_res/")
+
+path_plots = paste0(path_data, "/plots/")
+path_rdata = paste0(path_data, "/rdata/")
+path_meta = paste0(path_data, "/meta/")
+path_temp = paste0(path_data, "/temp/")
+path_output = paste0(path_data, "/output/")
+path_vis = paste0(path_data, "/vis/")
+
+
+# Set libraries ----------------------------------------------------------------
+#library(biodivTools) # devtools::install_github("environmentalinformatics-marburg/biodivTools")
+library(CAST)
+# library(corrplot)
+library(doParallel)
+library(grid)
+library(gridExtra)
+library(gpm) # devtools::install_github("environmentalinformatics-marburg/gpm")
+# library(hsdar)
+# library(lavaan)
+# library(rPointDB)
+# library(rgeos)
+# library(ggplot2)
+# library(mapview)
+# library(metTools) # devtools::install_github("environmentalinformatics-marburg/metTools")
+library(raster)
+# library(RStoolbox)
+library(reshape2)
+# library(rgdal)
+# library(satellite)
+# library(satelliteTools) # devtools::install_github("environmentalinformatics-marburg/satelliteTools")
+# library(semPlot)
+# library(sp)
+# library(spacetime)
+# library(vegan)
+# library(yaml)
+
+# Other settings ---------------------------------------------------------------
+source(filepath_source)
+
+rasterOptions(tmpdir = path_temp)
+
+# saga_cmd = "C:/OSGeo4W64/apps/saga-ltr/saga_cmd.exe"
+# initOTB("C:/OSGeo4W64/bin/")
+# initOTB("C:/OSGeo4W64/OTB-6.2.0-Win64/bin/")
+
+

src/000_setup_windows.R

@@ -0,0 +1,54 @@
+# Set environment for environmental information systems analysis
+require(envimaR)
+
+root_folder = path.expand("~/plygrnd/hySpecVisKili/")
+fcts_folder = file.path(root_folder, "hySpecVisKili/src/functions/")
+
+project_folders = c("data/",
+ "data/biodiv",
+ "data/020_hypspec_org/",
+ "data/025_hypspec_aio/",
+ "data/030_hypspec_nrm/",
+ "data/040_hypspec_vegidcs/",
+ "data/050_hypspec_kmdc/",
+ "data/060_hypspec_raoq/",
+ "data/070_hypspec_glcm/",
+ "data/090_hypspec_pred/",
+ "data/100_comb_gpm_sr/",
+ "data/110_model_gpm_sr/",
+ "data/120_compile_analysis_sr/",
+ "data/200_comb_gpm_sr_res/",
+ "data/210_model_gpm_sr_res/",
+ "data/220_compile_analysis_sr_elev_res/",
+ "data/300_model_gpm_sr_indp/",
+ "data/310_comb_gpm_sr_elev_res_indp/",
+ "data/320_model_gpm_sr_elev_res_indp/",
+ "data/500_analysis_sr/",
+ "data/510_analysis_sr_elev_res/",
+ "data/plots/",
+ "data/rdata/",
+ "data/meta/",
+ "data/output/",
+ "data/vis/",
+ "data/temp/")
+
+libs = c("biodivTools", "CAST", "corrplot", "doParallel", "grid", "gridExtra", 
+ "gpm", "ggplot2", "hsdar", "rgeos", "ggplot2", "mapview", 
+ "raster", "RStoolbox", "reshape2", "rgdal", "satelliteTools", "sp", 
+ "spacetime", "vegan")
+
+envrmt = createEnvi(root_folder = root_folder,
+ fcts_folder = fcts_folder,
+ folders = project_folders, 
+ path_prefix = "path_", libs = libs,
+ alt_env_id = "COMPUTERNAME", alt_env_value = "PCRZP",
+ alt_env_root_folder = "F:\\BEN\\edu")
+
+# More settings
+rasterOptions(tmpdir = envrmt$path_temp)
+mapviewOptions(basemaps = mapviewGetOption("basemaps")[c(3, 1:2, 4:5)])
+saga_cmd = "C:/OSGeo4W64/apps/saga-ltr/saga_cmd.exe"
+# initOTB("C:/OSGeo4W64/bin/")
+initOTB("C:/OSGeo4W64/OTB-6.2.0-Win64/bin/")
+
+

src/010_biodiv_preprocessing.R

@@ -0,0 +1,179 @@
+# Preprocess biodiversity observations.
+
+source("D:/plygrnd/hySpecVisKili/hySpecVisKili/src/000_set_environment.R")
+
+
+# Read species richness dataset (Peters et al. 2016)
+bd = read.table(paste0(path_biodiv, "Biodiversity_Data_Marcel.csv"),
+ header = TRUE, sep = ";", dec = ",")
+
+saveRDS(as.character(bd$plotID), file = paste0(path_biodiv, "biodiv_plots.rds"))
+saveRDS(bd, file = paste0(path_biodiv, "biodiv.rds"))
+
+# Read animal species dataset
+ad = read.table(paste0(path_biodiv, "animals_plotIDcomplete_Syn1.csv"),
+ header = TRUE, sep = ";", dec = ",")
+
+ad_diet = lapply(unique(ad$taxon), function(t){
+ data.frame(taxon = t,
+ ad$diet) 
+})
+
+
+
+# Set species number to 0/1 and reduce to complete cases
+ad[, 14:73][!is.na(ad[, 14:73]) & ad[, 14:73]>1] = 1
+adc = ad[complete.cases(ad[, 14:73]),]
+
+# Split into taxon groups
+adc_taxl = lapply(unique(adc$taxon), function(t){
+ act_level = adc[adc$taxon == t, ]
+ act_grp = act_level[, 14:73]
+ rownames(act_grp) = act_level$species
+ return(t(act_grp))
+})
+
+adc_taxl_all = adc[, 14:73]
+rownames(adc_all) = adc$species
+
+adc_taxl = c(list(t(adc_taxl_all)), adc_taxl)
+names(adc_taxl) = c("Animals", as.character(unique(adc$taxon)))
+
+saveRDS(adc_taxl, file = paste0(path_biodiv, "adc_taxl.rds"))
+
+
+
+# Split into trophic levels
+adc_tlevels = lapply(unique(adc$diet), function(d){
+ act_level = adc[adc$diet == d, ]
+ act_grp = act_level[, 14:73]
+ rownames(act_grp) = act_level$species
+ return(t(act_grp))
+})
+
+adc_all =adc[, 14:73]
+rownames(adc_all) = adc$species
+
+adc_tlevels = c(list(t(adc_all)), adc_tlevels)
+names(adc_tlevels) = c("Animals", as.character(unique(adc$diet)))
+
+saveRDS(adc_tlevels, file = paste0(path_biodiv, "adc_tlevels.rds"))
+
+
+
+# Compute species richness from taxon groups and tropich levels
+adc_taxl_sr = data.frame(plotID = rownames(adc_taxl[[1]]))
+for(i in seq(length(adc_taxl))){
+ adc_taxl_sr[, i+1] = rowSums(adc_taxl[[i]])
+ names(adc_taxl_sr)[i+1] = paste0("SR", tolower(names(adc_taxl[i])))
+}
+
+adc_tlevels_sr = data.frame(plotID = rownames(adc_tlevels[[1]]))
+for(i in seq(length(adc_tlevels))){
+ adc_tlevels_sr[, i+1] = rowSums(adc_tlevels[[i]])
+ names(adc_tlevels_sr)[i+1] = paste0("SR", tolower(names(adc_tlevels[i])))
+}
+
+adc_sr = merge(adc_tlevels_sr, 
+ adc_taxl_sr[, -grep("SRanimals", colnames(adc_taxl_sr))], 
+ by = "plotID")
+
+saveRDS(adc_sr, file = paste0(path_biodiv, "adc_sr.rds"))
+
+
+
+# Merge species richness data
+s_adc_sr = colnames(adc_sr[-grep("plotID", colnames(adc_sr))])
+
+colnames(bd)[which("SRsyrphids" == colnames(bd))] = "SRsyrphid_flies"
+colnames(bd)[which("SRbats" == colnames(bd))] = "SRinsectivorous_bats"
+colnames(bd)[which("SRmammals" == colnames(bd))] = "SRlarge_mammals"
+colnames(bd)[which("SRparasitoids" == colnames(bd))] = "SRparasitoid_wasps"
+colnames(bd)[which("SRotheraculeata" == colnames(bd))] = "SRaculeate_wasps"
+colnames(bd)[which("SRmillipedes" == colnames(bd))] = "SRmilipeds"
+colnames(bd)[which("SRsnails" == colnames(bd))] = "SRgastropods"
+
+s_bd = colnames(bd[, grep("SR", colnames(bd))])
+# which(tolower(s_adc_taxl_sr) %in% tolower(s_bd))
+
+species_richness = merge(bd, adc_sr, by = "plotID")
+species_richness = species_richness[, -grep("\\.x", colnames(species_richness))]
+colnames(species_richness)[grep("\\.y", colnames(species_richness))] = 
+ substr(colnames(species_richness)[grep("\\.y", colnames(species_richness))], 1, 
+ (nchar(colnames(species_richness)[grep("\\.y", colnames(species_richness))])-2))
+
+saveRDS(species_richness, file = paste0(path_biodiv, "species_richness.rds"))
+
+
+
+# Compute community composition using detrended correspondence analysis
+species_composition_dcor = lapply(adc_tlevels, function(l){
+ l = l[rowSums(l) > 0, ]
+ decorana(l)
+})
+names(species_composition_dcor) = names(adc_tlevels)
+# for(i in seq(5)) plot(species_composition_dcor[[i]], display = "sites")
+
+saveRDS(species_composition_dcor, file = paste0(path_biodiv, "species_composition_dcor.rds"))
+
+
+
+# Compute relative network using PCA
+tlevels = colnames(species_richness)[
+ seq(grep("SRanimals", colnames(species_richness)), 
+ grep("SRanimals", colnames(species_richness))+4)]
+
+adn_matrix = matrix(ncol = 5, nrow = 60)
+for(i in seq(5)){
+ adn_matrix[, i] = species_richness[, tlevels[i]]
+}
+rownames(adn_matrix) <- species_richness$plotID
+colnames(adn_matrix ) <- tlevels
+
+species_network_pca <- princomp(adn_matrix[,-1], cor=T)
+
+# biplot(species_network_pca, choices = 2:3)
+# summary(species_network_pca)
+
+saveRDS(species_network_pca, file = paste0(path_biodiv, "species_network_pca.rds"))
+
+
+# comb_grigusovaine diversity for Grigusova et al. 2019
+comb_grigusova = data.frame(species_network_pca$scores)
+colnames(comb_grigusova) = paste0("sn_pca", seq(4))
+comb_grigusova$plotID = rownames(species_network_pca$scores)
+
+comb_grigusova = merge(adc_sr, comb_grigusova, by = c("plotID"), all.x = TRUE, all.y = TRUE)
+
+for(i in seq(length(species_composition_dcor))){
+ act = data.frame(species_composition_dcor[[i]]$rproj)
+ colnames(act) = paste0("sn_dca", seq(4), "_", names(species_composition_dcor[i]))
+ act$plotID = rownames(act)
+ comb_grigusova = merge(comb_grigusova, act, by = c("plotID"), all.x = TRUE, all.y = TRUE)
+}
+
+comb_grigusova = droplevels(comb_grigusova)
+
+saveRDS(comb_grigusova, file = paste0(path_biodiv, "comb_grigusova.rds"))
+
+# Cross check
+# sort(colnames(adc_taxl_sr))
+# sort(colnames(bd[, c(1, grep("SR", colnames(bd)))]))
+# 
+# test = merge(adc_taxl_sr, bd, by = "plotID", all = TRUE)
+# 
+# test_df = lapply(grep("\\.x", colnames(test)), function(t){
+# x = colnames(test)[t]
+# y = paste0(substr(x, 1, (nchar(x)-1)), "y")
+# act_df = data.frame(plotID = test$plotID, test[, x], test[, y])
+# colnames(act_df) = c("plotId", x, y)
+# return(act_df)
+# })
+# test_df
+
+
+
+
+
+
+