diff --git a/code/Trait-Environment-Relationships2.R b/code/Trait-Environment-Relationships2.R
deleted file mode 100644
index 6ad4f3259eda6931784cb2d8665e33a0586580c8..0000000000000000000000000000000000000000
--- a/code/Trait-Environment-Relationships2.R
+++ /dev/null
@@ -1,2008 +0,0 @@
-getwd()
-###### Traits #####
-# TRY 2.0
-load("data/TRY.all.mean.sd.2.Rdata")
-str(TRY.all.mean.sd.2)
-dim(TRY.all.mean.sd.2)
-# 40791 obs. of 38 variables
-# tests
-any(is.na(TRY.all.mean.sd.2$SLA.sd)) #T
-
-# TRY 3.0
-TRY3.0 <- read.csv("data/Export_sPlot_2016_09_12.csv")
-str(TRY3.0)
-levels(TRY3.0$Species)[100] # "Abronia villosa", just an example
-which(TRY3.0$Species==levels(TRY3.0$Species)[100]) # an example of one species with several (7) observations
-TRY3.0[which(TRY3.0$Species==levels(TRY3.0$Species)[100]),c(1:19)] # trait values differ between observations
-any(is.na(TRY3.0[,c(2:19)])) #F
-any(TRY3.0[,c(2:19)]==0) #F
-
-for (i in 2:19){
- TRY3.0[,i] <- log(TRY3.0[,i])
-}
-any(is.na(TRY3.0[,c(2:19)])) #F
-
-TRY.all.n.3 <- aggregate(TRY3.0[,1], by=list(TRY3.0$Species), FUN=length)
-head(TRY.all.n.3)
-str(TRY.all.n.3) #59319 obs. of 2 variables:
-names(TRY.all.n.3)
-names(TRY.all.n.3) <- c("StandSpeciesName","n")
-
-TRY.all.mean.3 <- aggregate(TRY3.0[,c(2:19)], by=list(TRY3.0$Species), FUN=mean)
-str(TRY.all.mean.3) #59319 obs. of 19 variables:
-names(TRY.all.mean.3)
-' [1] "Group.1" "X1" "X4" "X11" "X13" "X14" "X15" "X18" "X26" "X27" "X47"
-[12] "X50" "X56" "X78" "X138" "X163" "X169" "X237" "X282" '
-names(TRY.all.mean.3) <- c("StandSpeciesName","LeafArea.mean", "StemDens.mean", "SLA.mean", "LeafC.perdrymass.mean",
- "LeafN.mean","LeafP.mean", "PlantHeight.mean", "SeedMass.mean", "Seed.length.mean",
- "LDMC.mean","LeafNperArea.mean", "LeafNPratio.mean", "Leaf.delta.15N.mean",
- "Seed.num.rep.unit.mean", "Leaffreshmass.mean", "Stem.cond.dens.mean",
- "Disp.unit.leng.mean", "Wood.vessel.length.mean")
-any(is.na(TRY.all.mean.3$SLA.mean)) #F
-
-TRY.all.sd.3 <- aggregate(TRY3.0[,c(2:19)], by=list(TRY3.0$Species), FUN=sd)
-str(TRY.all.sd.3) #59319 obs. of 19 variables:
-names(TRY.all.sd.3)
-names(TRY.all.sd.3) <- c("StandSpeciesName","LeafArea.sd", "StemDens.sd", "SLA.sd", "LeafC.perdrymass.sd",
- "LeafN.sd","LeafP.sd", "PlantHeight.sd", "SeedMass.sd", "Seed.length.sd",
- "LDMC.sd","LeafNperArea.sd", "LeafNPratio.sd", "Leaf.delta.15N.sd",
- "Seed.num.rep.unit.sd", "Leaffreshmass.sd", "Stem.cond.dens.sd",
- "Disp.unit.leng.sd", "Wood.vessel.length.sd")
-any(is.na(TRY.all.sd.3$SLA.sd)) #T
-TRY.all.mean.sd.3 <- data.frame(TRY.all.n.3,TRY.all.mean.3[,c(2:19)],TRY.all.sd.3[,c(2:19)])
-str(TRY.all.mean.sd.3) #59319 obs. of 38 variables
-'data.frame: 59319 obs. of 38 variables:
-$ StandSpeciesName : Factor w/ 59319 levels "Aa sp","Abarema adenophora",..: 1 2 3 4 5 6 7 8 9 10 ...
-$ n : int 1 4 1 1 38 4 86 7 5 95 ...
-$ LeafArea.mean : num 6.61 6.96 7.04 7.11 6.75 ...
-$ StemDens.mean : num -0.822 -0.427 -0.57 -0.619 -0.543 ...
-$ SLA.mean : num 2.24 2.38 2.41 2.56 2.65 ...
-$ LeafC.perdrymass.mean : num 6.17 6.26 6.19 6.17 6.18 ...
-$ LeafN.mean : num 2.92 3.19 3.25 3.41 3.23 ...
-$ LeafP.mean : num -0.0143 0.3403 -0.1557 -0.0999 -0.217 ...
-$ PlantHeight.mean : num -0.498 2.874 2.544 2.293 2.394 ...
-$ SeedMass.mean : num -4.08 4.49 4.27 4.23 4.13 ...
-$ Seed.length.mean : num -0.317 2.021 1.88 1.92 1.912 ...
-$ LDMC.mean : num -1.471 -0.798 -0.938 -1.087 -1.111 ...
-$ LeafNperArea.mean : num 0.896 0.834 0.838 0.804 0.617 ...
-$ LeafNPratio.mean : num 2.54 3.26 3.49 3.55 3.57 ...
-$ Leaf.delta.15N.mean : num 0.38 1.258 0.925 -0.169 0.719 ...
-$ Seed.num.rep.unit.mean : num 9.56 1.64 1.39 1.18 1.32 ...
-$ Leaffreshmass.mean : num -1.3303 -0.2516 -0.0799 -0.0889 -0.4985 ...
-$ Stem.cond.dens.mean : num 3.67 2.14 1.94 2.08 2.06 ...
-$ Disp.unit.leng.mean : num -0.555 2.916 2.676 2.689 2.842 ...
-$ Wood.vessel.length.mean: num 6.13 5.73 6 6.06 6.21 ...
-$ LeafArea.sd : num NA 0.0303 NA NA 0.0556 ...
-$ StemDens.sd : num NA 0.1101 NA NA 0.0123 ...
-$ SLA.sd : num NA 0.1063 NA NA 0.0695 ...
-$ LeafC.perdrymass.sd : num NA 0.0086 NA NA 0.00511 ...
-$ LeafN.sd : num NA 0.0638 NA NA 0.0113 ...
-$ LeafP.sd : num NA 0.3936 NA NA 0.0232 ...
-$ PlantHeight.sd : num NA 0.0626 NA NA 0.1806 ...
-$ SeedMass.sd : num NA 0.0315 NA NA 0.0787 ...
-$ Seed.length.sd : num NA 0.0457 NA NA 0.0193 ...
-$ LDMC.sd : num NA 0.0553 NA NA 0.0573 ...
-$ LeafNperArea.sd : num NA 0.098 NA NA 0.068 ...
-$ LeafNPratio.sd : num NA 0.291 NA NA 0.0304 ...
-$ Leaf.delta.15N.sd : num NA 0.1651 NA NA 0.0398 ...
-$ Seed.num.rep.unit.sd : num NA 0.211 NA NA 0.259 ...
-$ Leaffreshmass.sd : num NA 0.0748 NA NA 0.0571 ...
-$ Stem.cond.dens.sd : num NA 0.1389 NA NA 0.0349 ...
-$ Disp.unit.leng.sd : num NA 0.0927 NA NA 0.0293 ...
-$ Wood.vessel.length.sd : num NA 0.1466 NA NA 0.0271 ...'
-
-plot(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.3)
-#summary(lm(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.2))
-#abline(lm(SLA.mean~LDMC.mean, data=TRY.all.mean.sd.2),col="blue", lwd=2)
-#sqrt(0.3415)
-
-#TRY.all.mean.sd.3 <- TRY.all.mean.sd.3[,-39]
-save(TRY.all.mean.sd.3, file="TRY.all.mean.sd.3.Rdata")
-
-###### Backbone #####
-load("data/backbone.v.2.splot.try3.Rdata")
-# alternatively:
-#sPlot3 <- read.csv(paste("data/", "backbone.v.2.splot.try3.csv",
-# sep = ""), sep=",", quote = "\"'")
-str(backbone.splot.try3)
-# data.frame': 122901 obs. of 33 variables:
-
-###### Harmonized plot data ######
-DT <- read.csv("data/plots_without_unified_cover_scales_20160120a.csv")
-DT <- read.csv("data/plots_without_unified_cover_scales_20160120a.csv")
-library(data.table)
-DT <- data.table(DT)
-colnames(DT)
-colnames(DT)[9] <- "Cover"
-setkey(DT,PlotObservationID)
-
-###### Header data #######
-splot.header <- read.csv("data/sPlot_14_04_2015_coord.csv")
-str(splot.header)
-'data.frame: 1117940 obs. of 3 variables:
-$ PlotObservationID: int 1 2 3 4 5 6 7 8 9 10 ...
-$ Longitude : num -154 -154 -154 -154 -154 ...
-$ Latitude : num 62.4 62.4 62.4 62.4 62.4 ...'
-
-##### Bioclim data ####
-splot.bioclim <- read.csv(paste("data/","sPlot.csv", sep = ""), sep = ";", na.strings=c("","NA"))
-str(splot.bioclim)
-#data.frame': 1117382 obs. of 52 variables:
-
-####################
-##### Matching #####
-#index6 <- match(TRY.all.mean.sd.2$StandSpeciesName,backbone.splot.try3$names.sPlot.TRY)
-index6 <- match(TRY.all.mean.sd.3$StandSpeciesName,backbone.splot.try3$names.sPlot.TRY)
-any(is.na(index6)) #T
-# not all species in the trait list are also in the backbone
-TRY.all.mean.sd.3$StandSpeciesName[is.na(index6)]
-# 16 species in TRY are not in the backbone
-'[1] Albizia greveanaÊ AlnusÊgymnothyrsus Carapanaœba sp
-[4] ChionanthusÊsp Clerodendron nudiflorumÊ Copa’ba sp
-[7] Copa’fera sp Corynocarpus laevigatusÊ Enterospermum madagascarienseÊ
-[10] Hierochlo‘ odorata Ing‡ sp Norrisia maiorÊ
-[13] NothofagusÊbrassospora QuercusÊsp ShoreaÊanthoshorea
-[16] ShoreaÊsp '
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Albizia greveanaÊ",] # n=10
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Albizia greveana",] # 0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Albizia greveanaÊ"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Albizia greveana"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="AlnusÊgymnothyrsus"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Alnus gymnothyrsus"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Carapanaœba sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Carapanauba sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Carapanaúba sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="ChionanthusÊsp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Chionanthus sp"] #1
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="ChionanthusÊsp",] # n=10
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Chionanthus sp",] # n=15
-DT$species[DT$species=="Chionanthus sp" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Clerodendron nudiflorumÊ"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Clerodendron nudiflorum"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Copa’ba sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Copaiba sp"] #1
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Copa’ba sp",] # n=0, strange, there should be at least one
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Copaiba sp",] # n=1
-DT$species[DT$species=="Copaiba sp" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Copa’fera sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Copaifera sp"] #1
-DT$species[DT$species=="Copaifera sp" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Corynocarpus laevigatusÊ"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Corynocarpus laevigatus"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Enterospermum madagascarienseÊ"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Enterospermum madagascariense"] #1
-DT$species[DT$species=="Enterospermum madagascariense" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Hierochlo‘ odorata"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Hierochloe odorata"] #1
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Hierochlo‘ odorata",] # n=0, strange, there should be at least one
-TRY.all.mean.sd.3[TRY.all.mean.sd.3$StandSpeciesName=="Hierochloe odorata",] # n=63
-DT$species[DT$species=="Hierochloe odorata" & !is.na(DT$species)] #545
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Ing‡ sp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Inga sp"] #1
-DT$species[DT$species=="Inga sp" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Norrisia maiorÊ"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Norrisia maior"] #1
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="NothofagusÊbrassospora"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Nothofagus brassospora"] #1
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="QuercusÊsp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Quercus sp"] #1
-DT$species[DT$species=="Quercus sp" & !is.na(DT$species)] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="ShoreaÊanthoshorea"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Shorea anthoshorea"] #1
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="ShoreaÊsp"] #0
-backbone.splot.try3$names.sPlot.TRY[backbone.splot.try3$names.sPlot.TRY=="Shorea sp"] #1
-DT$species[DT$species=="Shorea sp" & !is.na(DT$species)] #0
-# none of the species is in sPlot, thus we can take them out
-
-### data cleaning ###
-length(TRY.all.mean.sd.3$StandSpeciesName[is.na(TRY.all.mean.sd.3$StandSpeciesName)])
-# No record is NA!
-#TRY.all.mean.sd.2[which(is.na(TRY.all.mean.sd.2$StandSpeciesName)),]
-# record no. 7947 has StandSpeciesName = NA
-#TRY.all.mean.sd.2 <- TRY.all.mean.sd.2[-which(is.na(TRY.all.mean.sd.2$StandSpeciesName)),]
-
-TRY.all.mean.sd.3$species <- backbone.splot.try3$name.short.correct[index6]
-dim(TRY.all.mean.sd.3) # 59319 39
-
-
-index1 <- match(DT$PlotObservationID, splot.header$PlotObservationID)
-any(is.na(index1)) #F
-length(index1) # 24241941
-length(unique(index1[!is.na(index1)])) #1117940
-
-index2 <- match(splot.header$PlotObservationID,DT$PlotObservationID)
-any(is.na(index2)) #F
-length(index2) # 1117940
-length(unique(index2[!is.na(index2)])) # 1117940
-dim(splot.header)
-# 1117940 3
-identical(as.character(splot.header$PlotObservationID),as.character(DT$PlotObservationID[index2]))
-# T
-
-any(is.na(DT$species)) #T
-# that should not be the case!
-length(DT$species[!is.na(DT$species)]) #24221565
-length(DT$species) #24241941
-24241941-24221565 # 20376 NA names!!!
-# it gives nonsense to match them with traits
-index7 <- match(DT$species,TRY.all.mean.sd.3$species)
-length(index7) #24241941
-length(index7[!is.na(index7)])
-# 21839463, with TRY3.0
-# 21040927, with TRY2.0
-# 19841429, with first TRY version
-24241941 - 21839463 # 2402478 entries have no traits
-(24241941 - 2402478)/24241941*100
-# which are 90.08958% of all entries
-# previously: 86.79555% with TRY2.0
-# previously: 18.15247 % and with first TRY version
-
-# reduce DT, and splot.header
-#splot.header <- splot.header[!is.na(DT$species[index2]),]
-#dim(splot.header) #1117476 3
-DT <- DT[!is.na(DT$species),]
-index7 <- match(DT$species,TRY.all.mean.sd.3$species)
-length(index7) #24221565
-
-
-### CWM ###
-mean(TRY.all.mean.sd.3$SLA.mean) # 2.634988
-min(TRY.all.mean.sd.3$SLA.mean) # 0.03072317
-max(TRY.all.mean.sd.3$SLA.mean) # 5.182365
-# example
-DT$trait <- NA
-DT$trait <- TRY.all.mean.sd.3$SLA.mean[index7]
-length(DT$trait[!is.na(DT$trait)]) # 21819087
-length(DT$trait[is.na(DT$trait)]) # 2402478
-str(DT)
-
-colnames(TRY.all.mean.sd.3)
-CWM2 <- DT[,list(CWM.SLA = weighted.mean(trait,Cover,na.rm = T)),by=PlotObservationID]
-# I have checked that we do not need to prefilter only those entries that have a trait value
-# this works fine
-str(CWM2)
-dim(CWM2) #1117898
-which(colnames(TRY.all.mean.sd.3)=="LeafArea.mean") # 3
-which(colnames(TRY.all.mean.sd.3)=="Wood.vessel.length.mean") # 20
-CWM <- array(NA,c(dim(CWM2)[1],18),dimnames=list(CWM2$PlotObservationID,
- colnames(TRY.all.mean.sd.3)[3:20]))
-rm(CWM2)
-for (i in 1:18){
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.3[index7,i+2]
- CWM[,i] <- DT[,list(CWM.trait= weighted.mean(trait,Relative.cover,na.rm = T)),by=PlotObservationID]$CWM.trait
-
-}
-str(CWM) #num [1:1117898, 1:18]
-CWM[1:20,]
-tail(CWM)
-write.csv(CWM,file = "CWM_TRY3.csv", row.names = T)
-
-CWM <- read.csv("C://Daten//iDiv2//splot2//CWM_TRY3.csv")
-str(CWM) #1117898 obs. of 19 variables
-dimnames(CWM)[[1]] <- CWM$X
-CWM <- CWM[,-1]
-tail(CWM)
-
-##### MAP ####
-index3 <- match(dimnames(CWM)[[1]],splot.header$PlotObservationID)
-length(index3) #1117898
-any(is.na(index3)) # F
-#length(index3[is.na(index3)]) #422
-splot.header2 <- splot.header[index3,]
-'CWM1$Longitude <- splot.header$Longitude[index10]
-CWM1$Latitude <- splot.header$Latitude[index10]
-'
-str(CWM) # num [1:1117898, 1:18]
-CWM <- CWM[!(is.na(splot.header2$Longitude)|is.na(splot.header2$Latitude)|is.na(CWM[,"SLA.mean"])),]
-str(CWM) # num [1:1112287, 1:18]
-1112287/1117898*100 # 99.49808% of plots have traits, Long and Lat
-any(is.na(CWM)) #F
-# match header again
-index4 <- match(dimnames(CWM)[[1]],splot.header$PlotObservationID)
-length(index4) #1112287
-any(is.na(index4)) # F
-splot.header2 <- splot.header[index4,]
-# now in the same sequence as CWM
-
-library(raster)
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-coords <- cbind(splot.header2$Longitude, splot.header2$Latitude)
-coords <- SpatialPoints(coords, proj4string=CRSlonlat)
-map.CWM <- SpatialPointsDataFrame(coords, as.data.frame(CWM), proj4string=CRSlonlat)
-str(map.CWM)
-
-library(maptools)
-data(wrld_simpl)
-plot(wrld_simpl)
-points(coordinates(map.CWM), cex=0.1, col="red", pch=16)
-# sPlot2withSLA.png
-max(coordinates(map.CWM)[,1]) # 179.5901
-max(coordinates(map.CWM)[,2]) # 80.14912
-min(coordinates(map.CWM)[,1]) # -162.7414
-min(coordinates(map.CWM)[,2]) # -64.78
-
-
-r <- raster(nrows=180, ncols=360, xmn=-180, xmx=180, ymn=-90, ymx=90)
-str(r)
-raster.map.SLA <- rasterize(map.CWM,r,map.CWM$SLA.mean, run="mean")
-str(raster.map.SLA)
-plot(wrld_simpl)
-plot(raster.map.SLA, add=T)
-# SLAgridded1.png
-
-library(rworldmap)
-raster.map.SLA.grid <- as(raster.map.SLA, 'SpatialGridDataFrame')
-library(RColorBrewer)
-# Set a color ramp
-rf <- brewer.pal(8, 'Spectral')
-par(mar=c(0, 0, 0, 8) + 0.1)
-mapParams <- mapGriddedData(raster.map.SLA.grid,addLegend=F,
- numCats = 8, colourPalette=rf)
-str(mapParams)
-mapParams$cutVector
-# addMapLegend(plottedData=raster.map.SLA.grid,legendShrink=0.5)
-#do.call(addMapLegend, c( mapParams, legendLabels="all", legendWidth=0.5,
-# legendIntervals="data", legendMar = 2 ) )
-#do.call(addMapLegend, c( mapParams, legendLabels="all", legendWidth=0.3,
-# legendIntervals="page", legendMar = 4, legendShrink=1.7,horizontal=F))
-do.call(addMapLegend, c( mapParams, legendLabels="all", legendWidth=0.3,
- legendIntervals="page", legendMar =2, horizontal=F))
-
-#round(as.numeric(exp(mapParams$cutVector)),1)
-mtext(round(as.numeric(exp(mapParams$cutVector)),1),side=4,line=-1.5,
- at=seq(-200,230, by=51.3), las=2)
-#mtext("ln(SLA)",side=1,line=0, at=230, las=1)
-#mtext("ln(SLA)",side=4,line=-2,at=-170, las=2)
-# SLAgridded2.png
-par(mar=c(5, 4, 4, 2) + 0.1)
-
-############################################################
-# Regression with Bioclim
-############################################################
-'BIO1 = Annual Mean Temperature
-BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
-BIO3 = Isothermality (BIO2/BIO7) (* 100)
-BIO4 = Temperature Seasonality (standard deviation *100)
-BIO5 = Max Temperature of Warmest Month
-BIO6 = Min Temperature of Coldest Month
-BIO7 = Temperature Annual Range (BIO5-BIO6)
-BIO8 = Mean Temperature of Wettest Quarter
-BIO9 = Mean Temperature of Driest Quarter
-BIO10 = Mean Temperature of Warmest Quarter
-BIO11 = Mean Temperature of Coldest Quarter
-BIO12 = Annual Precipitation
-BIO13 = Precipitation of Wettest Month
-BIO14 = Precipitation of Driest Month
-BIO15 = Precipitation Seasonality (Coefficient of Variation)
-BIO16 = Precipitation of Wettest Quarter
-BIO17 = Precipitation of Driest Quarter
-BIO18 = Precipitation of Warmest Quarter
-BIO19 = Precipitation of Coldest Quarter'
-
-index11 <- match(splot.header2$PlotObservationID,splot.bioclim$PLOT_ID)
-length(index11) # 1112287
-length(index11[!is.na(index11)]) # 1112287
-
-#index12 <- match(CWM2$PlotObservationID, splot.header$PlotObservationID)
-#any(is.na(index12)) #F
-str(splot.bioclim$BIO_01)
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11])
-model2 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11]+I(splot.bioclim$BIO_01[index11]^2))
-summary(model1)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.929e+00 8.700e-04 3367 <2e-16 ***
-splot.bioclim$BIO_01[index11] -9.720e-04 8.603e-06 -113 <2e-16 ***
-Multiple R-squared: 0.01135, Adjusted R-squared: 0.01135
-F-statistic: 1.277e+04 on 1 and 1112285 DF, p-value: < 2.2e-16
-'
-summary(model2)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.791e+00 1.059e-03 2635.2 <2e-16 ***
-splot.bioclim$BIO_01[index11] 2.446e-03 1.765e-05 138.6 <2e-16 ***
-I(splot.bioclim$BIO_01[index11]^2) -1.673e-05 7.593e-08 -220.3 <2e-16 ***
-Multiple R-squared: 0.05269
-'
-abline(model1, lwd=3, col="blue")
-model2$coefficients
-x.new <- seq(-100,300,1)
-y.new <- model2$coefficients[1]+model2$coefficients[2]*x.new+model2$coefficients[3]*x.new^2
-lines(x.new,y.new, lwd=3, col="blue")
-# SLA_BIO1.png
-
-# devide in quantiles
-quantile(splot.bioclim$BIO_01[index11],probs = seq(0, 1, 1/3))
-' 0% 33.33333% 66.66667% 100%
- -227 79 99 303 '
-# -> most points (i.e. one third) are in the temperate zone between 7.9 and 9.9°C
-x1 <- splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=79]
-model1a <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]<=79]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=79])
-summary(model1a)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.688e+00 1.178e-03 2282 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] <= 79] 2.826e-03 2.004e-05 141 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.382 on 372215 degrees of freedom
-Multiple R-squared: 0.05069, Adjusted R-squared: 0.05068
-F-statistic: 1.987e+04 on 1 and 372215 DF, p-value: < 2.2e-16'
-model1b <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]<=99]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]<=99])
-summary(model1b)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.695e+00 1.112e-03 2422.8 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] <= 99] 2.709e-03 1.463e-05 185.1 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.3796 on 750485 degrees of freedom
-Multiple R-squared: 0.04367, Adjusted R-squared: 0.04367
-F-statistic: 3.427e+04 on 1 and 750485 DF, p-value: < 2.2e-16'
-model1c <- lm(CWM[,"SLA.mean"][splot.bioclim$BIO_01[index11]>99]~splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11]>99])
-summary(model1c)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 3.210e+00 2.425e-03 1323.3 <2e-16 ***
-splot.bioclim$BIO_01[index11][splot.bioclim$BIO_01[index11] > 99] -3.417e-03 1.739e-05 -196.5 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.4157 on 361798 degrees of freedom
-Multiple R-squared: 0.09646, Adjusted R-squared: 0.09646
-F-statistic: 3.862e+04 on 1 and 361798 DF, p-value: < 2.2e-16'
-
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-x.new <- seq(-100,79,1)
-y.new <- model1a$coefficients[1]+model1a$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="blue")
-x.new <- seq(79,99,1)
-y.new <- model1b$coefficients[1]+model1b$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="green")
-x.new <- seq(99,300,1)
-y.new <- model1c$coefficients[1]+model1c$coefficients[2]*x.new
-lines(x.new,y.new, lwd=3, col="red")
-
-### calculate a GAM
-library(gam)
-x <- splot.bioclim$BIO_01[index11]
-model3 <- gam(CWM[,"SLA.mean"]~lo(x, span=0.5))
-summary(model3)
-' Null Deviance: 207269.5 on 1111306 degrees of freedom
-Residual Deviance: 187452.7 on 1111302 degrees of freedom
-AIC: 1175898
-Anova for Parametric Effects
- Df Sum Sq Mean Sq F value Pr(>F)
-lo(x, span = 0.5) 1 3238 3238.0 19197 < 2.2e-16 ***
-Residuals 1111302 187453 0.2
-Anova for Nonparametric Effects
- Npar Df Npar F Pr(F)
-(Intercept)
-lo(x, span = 0.5) 2.7 37028 < 2.2e-16 ***
-'
-
-x.new <- seq(-100,300,1)
-y.new <- predict(model3, data.frame(x=x.new), se = F)
-# se=T does not work
-plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_1", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-lines(x.new,y.new, lwd=3, col="blue")
-
-### N to P ratio
-model4 <- lm(CWM[,"LeafNPratio.mean"]~splot.bioclim$BIO_01[index11])
-summary(model4)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.334e+00 4.243e-04 5501 <2e-16 ***
-splot.bioclim$BIO_01[index11] 1.313e-03 4.196e-06 313 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.1986 on 1112285 degrees of freedom
-Multiple R-squared: 0.08096, Adjusted R-squared: 0.08096
-F-statistic: 9.799e+04 on 1 and 1112285 DF, p-value: < 2.2e-16
-'
-plot(CWM[,"LeafNPratio.mean"]~splot.bioclim$BIO_01[index11],
- cex=0.05, xlab="BIO_01", ylab="CWM ln(Leaf N to P ratio)",cex.lab=1.5, cex.axis=1.5)
-abline(model4, lwd=3, col="blue")
-
-### turn the regression around
-# does not make sense
-model5 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"SLA.mean"])
-summary(model5)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 127.41676 0.28026 454.6 <2e-16 ***
-CWM[, "SLA.mean"] -12.97407 0.09769 -132.8 <2e-16 ***
-Multiple R-squared: 0.01562, Adjusted R-squared: 0.01562 '
-
-model6 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"LeafNPratio.mean"])
-summary(model6)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) -72.7479 0.4591 -158.5 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 67.2535 0.1883 357.2 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 42.46 on 1111305 degrees of freedom
-Multiple R-squared: 0.103, Adjusted R-squared: 0.103 '
-
-model7 <- lm(splot.bioclim$BIO_01[index11]~CWM[,"SLA.mean"]+CWM[,"LeafNPratio.mean"])
-summary(model7)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) -41.59378 0.55839 -74.49 <2e-16 ***
-CWM[, "SLA.mean"] -9.08978 0.09355 -97.17 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 65.04103 0.18885 344.40 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 42.28 on 1111304 degrees of freedom
-Multiple R-squared: 0.1106, Adjusted R-squared: 0.1106 '
-x <- splot.bioclim$BIO_01[index11]
-model8 <- lm(x~CWM[,"LeafNPratio.mean"]*CWM[,"SLA.mean"])
-summary(model8)
-'Coefficients:
- Estimate Std. Error t value Pr(>|t|)
-(Intercept) -432.2500 2.2645 -190.9 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 227.1681 0.9304 244.2 <2e-16 ***
-CWM[, "SLA.mean"] 135.8572 0.8202 165.6 <2e-16 ***
-CWM[, "LeafNPratio.mean"]:CWM[, "SLA.mean"] -60.2280 0.3386 -177.9 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 41.69 on 1111303 degrees of freedom
-Multiple R-squared: 0.1352, Adjusted R-squared: 0.1352 '
-
-library(scatterplot3d) # für 3D-Plots
-s3d1 <- scatterplot3d(CWM[,"SLA.mean"],CWM[,"LeafNPratio.mean"], x,
- xlab="CWM ln(Leaf N to P ratio)", ylab="CWM ln(SLA)", zlab="BIO_1",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p")
-#fit <- lm(AnzF ~ Prop.n+Prop.ba,data=data6[data6$Sp=="Qupe",])
-s3d1$plane3d(model8)
-# does not work
-### producing the plane with points
-min(CWM[,"SLA.mean"])
-max(CWM[,"SLA.mean"])
-x2 <- sort(rep(seq(0.1, 5, by=0.1),35))
-min(CWM[,"LeafNPratio.mean"])
-max(CWM[,"LeafNPratio.mean"])
-y2 <- rep(seq(1, 4.4, by=0.1),50)
-func <- function(x,y) { -432.2500 + 227.1681 * x + 135.8572*y -60.2280*x*y}
-z2 <- func(x2,y2)
-s3d1 <- scatterplot3d(CWM[,"SLA.mean"],CWM[,"LeafNPratio.mean"], x,
- xlab="CWM ln(Leaf N to P ratio)", ylab="CWM ln(SLA)", zlab="BIO_1",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=70,scale.y=0.5)
-s3d1$points3d(x2,y2,z2,col="blue", type="p", pch=20, cex=0.5)
-#
-x <- splot.bioclim$BIO_01[index11]
-model7 <- lm(CWM[,"SLA.mean"]~x*CWM[,"LeafNPratio.mean"])
-summary(model7)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 1.948e+00 9.377e-03 207.8 <2e-16 ***
-x 1.278e-02 7.985e-05 160.0 <2e-16 ***
-CWM[, "LeafNPratio.mean"] 3.916e-01 3.856e-03 101.6 <2e-16 ***
-x:CWM[, "LeafNPratio.mean"] -5.472e-03 3.166e-05 -172.8 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.4213 on 1111303 degrees of freedom
-Multiple R-squared: 0.04842, Adjusted R-squared: 0.04842 '
-min(x)
-max(x)
-x2 <- sort(rep(seq(-100, 290, by=10),35))
-min(CWM[,"LeafNPratio.mean"])
-max(CWM[,"LeafNPratio.mean"])
-y2 <- rep(seq(1, 4.4, by=0.1),40)
-func <- function(x,y) { 1.948e+00 + 1.278e-02 * x + 3.916e-01*y -5.472e-03*x*y}
-z2 <- func(x2,y2)
-s3d1 <- scatterplot3d(x,CWM[,"LeafNPratio.mean"], CWM[,"SLA.mean"],
- xlab="BIO_1", ylab="CWM ln(Leaf N to P ratio)", zlab="CWM ln(SLA)",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=50,scale.y=1)
-s3d1$points3d(x2,y2,z2,col="blue", type="p", pch=20, cex=0.5)
-# final plot
-# some BIO_01 and LeafNPratio combinations do not exist
-?cut
-x.cut <- cut(x, breaks=seq(-110, 290, by=10))
-y.cut <- cut(CWM[,"LeafNPratio.mean"],breaks=seq(0.9, 4.4, by=0.1))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #40 35
-z2.cut <- table(x2,y2)
-head(z2.cut)
-dim(z2.cut) #40 35
-
-#library(reshape2)
-z.cut <- melt(z.cut)
-str(z.cut) # 1400 obs. of 3 variables:
-head(z.cut)
-z2.cut <- melt(z2.cut)
-str(z2.cut) #1400 obs. of 3 variables:
-head(z2.cut)
-z2.cut$value <- 0
-z2.cut$value <- z.cut$value
-z2.cut$xy <- paste(z2.cut$x2,z2.cut$y2,sep="_")
-#z3.cut <- acast(z2.cut, x2~y2)
-#head(z3.cut)
-xy <- paste(x2,y2,sep="_")
-index20 <- match(xy,z2.cut$xy)
-values.present <- vector(mode="logical", length=length(x2))
-values.present <- z2.cut$value[index20] > 0
-table(values.present)
-'z.cut$value[z.cut$x.cut==x.cut[1] & z.cut$y.cut==y.cut[1]]
-values.present <- function(z.cut,x,y){
- z.cut$value[z.cut$x.cut==x & z.cut$y.cut==y]
-}
-x1 <- data.frame(x,x.cut, y=CWM[,"LeafNPratio.mean"],y.cut)
-head(x1)
-z1 <- apply(x1,1,FUN=values.present(z.cut, x1$x.cut, x1$y.cut))
-values.present(z.cut, x1$x.cut[1], x1$y.cut[1])
-'
-
-s3d1 <- scatterplot3d(x,CWM[,"LeafNPratio.mean"], CWM[,"SLA.mean"],
- xlab="BIO_1", ylab="CWM ln(Leaf N to P ratio)", zlab="CWM ln(SLA)",
- pch=1,cex.symbols=0.05,highlight.3d=T,type="p",angle=-10,scale.y=0.8)
-s3d1$points3d(x2[values.present],y2[values.present],z2[values.present],col="blue", type="p", pch=20, cex=0.5)
-
-
-#### Collect all CWM regressions ####
-names(splot.bioclim) # from 4 (BIO_01 to 49 T_DEC
-CWM.bioclim <- array(NA,c(46,18,5),dimnames=list(names(splot.bioclim)[4:49],dimnames(CWM)[[2]],
- c("p_lin","r2_lin","p_qua1", "p_qua2","r2_qua")))
-for (i in 1:46){
- for (j in 1:18){
- model1 <- lm(CWM[,j]~splot.bioclim[index11,i+3])
- CWM.bioclim[i,j,1] <- summary(model1)$coefficients[2,4]
- CWM.bioclim[i,j,2] <- summary(model1)$r.squared
- model2 <- lm(CWM[,j]~splot.bioclim[index11,i+3]+I(splot.bioclim[index11,i+3]^2))
- CWM.bioclim[i,j,3] <- summary(model2)$coefficients[2,4]
- CWM.bioclim[i,j,4] <- summary(model2)$coefficients[3,4]
- CWM.bioclim[i,j,5] <- summary(model2)$r.squared
- }
-}
-str(CWM.bioclim) #num [1:46, 1:18, 1:5]
-write.csv(data.frame(names(splot.bioclim)[4:49],CWM.bioclim[,,2]),file = "CWM_bioclim_r2_lin2.csv", row.names = FALSE)
-write.csv(data.frame(names(splot.bioclim)[4:49],CWM.bioclim[,,5]),file = "CWM_bioclim_r2_qua2.csv", row.names = FALSE)
-
-# for SLA
-which(CWM.bioclim[,3,2]==max(CWM.bioclim[,3,2])) # BIO_02
-which(CWM.bioclim[,3,5]==max(CWM.bioclim[,3,5])) # BIO_02
-which(CWM.bioclim[,,2]==max(CWM.bioclim[,,2]),arr.ind = TRUE)
-# row col
-#BIO_02 2 3
-CWM.bioclim[2,3,2] #0.1113908
-which(CWM.bioclim[,,5]==max(CWM.bioclim[,,5]),arr.ind = TRUE)
-# row col
-#PET 22 2
-CWM.bioclim[22,2,5] # 0.1407318
-names(splot.bioclim)[2+3] #"BIO_02"
-dimnames(CWM)[[2]][9] #"LeafNperArea.mean"
-plot(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per Area)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,9]~splot.bioclim[index11,2+3])
-abline(model1, lwd=3, col="blue")
-# MeanNperArea_BIO_02.png
-png(filename="MeanNperArea_BIO_02.png", width=800, height=800, units = "px", bg="white")
- plot(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per Area)",cex.lab=1.5, cex.axis=1.5)
- #model1 <- lm(CWM[,9]~splot.bioclim[index11,2+3])
- abline(model1, lwd=3, col="blue")
-dev.off() # creates a file
-png(filename="SLA_BIO_02.png", width=800, height=800, units = "px", bg="white")
- plot(CWM[,"SLA.mean"]~splot.bioclim$BIO_02[index11],
- cex=0.05, xlab="BIO_02", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
- model1 <- lm(CWM[,1]~splot.bioclim[index11,2+3])
-abline(model1, lwd=3, col="blue")
-dev.off() # creates a file
-
-#### PCA of all CWM and all Bioclim variables ####
-library(vegan)
-pca1 <- rda(CWM,scale=T)
-plot(pca1)
-head(summary(pca1),100)
-'Partitioning of correlations:
- Inertia Proportion
-Total 18 1
-Unconstrained 18 1
-
-Eigenvalues, and their contribution to the correlations
-
-Importance of components:
- PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13
-Eigenvalue 5.3194 3.5487 2.0227 1.37785 1.27883 1.00296 0.71963 0.66827 0.56198 0.49244 0.34211 0.26158 0.19518
-Proportion Explained 0.2955 0.1971 0.1124 0.07655 0.07105 0.05572 0.03998 0.03713 0.03122 0.02736 0.01901 0.01453 0.01084
-Cumulative Proportion 0.2955 0.4927 0.6050 0.68159 0.75264 0.80836 0.84834 0.88546 0.91668 0.94404 0.96305 0.97758 0.98842
- PC14 PC15 PC16 PC17 PC18
-Eigenvalue 0.07591 0.05484 0.05067 0.01766 0.009314
-Proportion Explained 0.00422 0.00305 0.00281 0.00098 0.000520
-Cumulative Proportion 0.99264 0.99569 0.99850 0.99948 1.000000'
-str(summary(pca1))
-#text(summary(pca1)$species[,c(1,2)],dimnames(summary(pca1)$species)[[1]])
-traitnames <- dimnames(summary(pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(pca1)$species[,c(1,2)]
-plotcoord <- summary(pca1)$sites[,c(1,2)]
-which(traitnames=="Seed.length") #9
-traitcoord[9,2] <- traitcoord[9,2]+0.3
-which(traitnames=="LDMC") #10
-traitcoord[10,2] <- traitcoord[10,2]-0.3
-ordiplot(pca1, type="n", ylim=c(-10,10))
-text(traitcoord,traitnames, col="red", cex=1)
-points(plotcoord*30, cex=0.05)
-#env1 <- envfit(pca1,splot.bioclim[index11,c(4:49)],permutations = 99)
-#with permutations = 999:
-#Error: cannot allocate vector of size 4.1 Gb
-plot(env1)
-png(filename="PCA1.png", width=1000, height=1000, units = "px", bg="white")
- ordiplot(pca1, type="n")
- text(traitcoord,traitnames)
- plot(env1)
-dev.off() # creates a file
-'rm(env1)
-rm(backbone.splot.try3)
-rm(model1)
-rm(model2)
-rm(pca1)
-rm(splot.header)
-rm(splot.header2)
-rm(traitdist)
-rm(coords)
-rm(TRY.all.mean.sd.2)
-rm(index7)
-rm(raster.map.SLA)
-rm(trait.vec)
-rm(FD)
-rm(splot.bioclim)
-'
-gc()
-
-### produce a PCA with kernel density
-# From Díaz et al. 2016, Nature
-library(ks)
-str(plotcoord) #num [1:1111307, 1:2]
-H <- Hpi(x=plotcoord) # optimal bandwidth estimation
-# does not work
-est<- kde(x=plotcoord, H=H, compute.cont=TRUE) # kernel density estimation
-
-library(raster)
-library(fBasics)
-# Color palette
-#plotcoord <- plotcoord*30
-min(plotcoord[,1]) #-0.3202043
-max(plotcoord[,1]) # 0.3078603
-min(plotcoord[,2]) #-0.5899272
-max(plotcoord[,2]) # 0.2530614
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #62 84
-z.cut <- melt(z.cut)
-str(z.cut) # 5208 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1112287
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30)) # T
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.992156
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-#values.present <- z.cut$value[index30]
-ordiplot(pca1, type="n")
-#points(plotcoord*30, cex=0.05, col=rainbowPalette(n=10,
-# s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1)[values.cut])
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-text(traitcoord,traitnames, col="red", cex=1)
-
-
-ext.factor <- 30
-r.pca <- raster(nrows=300, ncols=300, xmn=min(plotcoord[,1])*ext.factor, xmx=max(plotcoord[,2])*ext.factor,
- ymn=min(plotcoord[,1])*ext.factor, ymx=max(plotcoord[,2]*ext.factor))
-# attention, xmn and ymn as well as xmx and ymx have to be the same, otherwise the raster is shifted
-r.pca.raster <- rasterize(plotcoord*ext.factor, r.pca, fun="count")
-str(r.pca.raster)
-plot(log10(r.pca.raster), asp=0, col=seqPalette(n=10, name=c("Blues")),
- main="", cex.main=0.7, xlab="", ylab="", add=T)
-plot(log10(r.pca.raster), asp=0, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL)),
- main="", cex.main=0.7, xlab="", ylab="", add=T)
-plot(log10(r.pca.raster), asp=0, col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1),
- main="", cex.main=0.7, xlab="", ylab="", add=T, legend=F)
-
-'## graph for the PCA ##
-ordiplot(pca1, type="n", xlim=c(-10,13), ylim=c(-10,10))
-rasterImage(as.raster(r.pca.raster), -10, -10, 20, 10, interpolate = F,
- col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1))
-plot(log10(r.pca.raster), asp=0, col=rainbowPalette(n=10, s = 1, v = 1, start = 3/6, end = 5/6, alpha = 1),
- main="", cex.main=0.7, xlab="", ylab="", add=T, legend=F)
-text(traitcoord,traitnames, col="red", cex=1)
-'
-
-
-rda1 <- rda(CWM~.,splot.bioclim[index11,c(4:49)],scale=T)
-plot(rda1)
-head(summary(rda1), 100)
-rda1
-'Partitioning of correlations:
- Inertia Proportion
-Total 18.000 1.0000
-Constrained 2.321 0.1289
-Unconstrained 15.679 0.8711
-
-Eigenvalues, and their contribution to the correlations
-
-Importance of components:
-RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
-Eigenvalue 1.20610 0.57369 0.23626 0.15021 0.06339 0.02077 0.02014 0.01570 0.01164 0.008565 0.007286 0.002948
-Proportion Explained 0.06701 0.03187 0.01313 0.00835 0.00352 0.00115 0.00112 0.00087 0.00065 0.000480 0.000400 0.000160
-Cumulative Proportion 0.06701 0.09888 0.11200 0.12035 0.12387 0.12502 0.12614 0.12701 0.12766 0.128140 0.128540 0.128710
-RDA13 RDA14 RDA15 RDA16 RDA17 RDA18 PC1 PC2 PC3 PC4 PC5 PC6
-Eigenvalue 0.002155 0.00111 0.0002949 0.0001558 0.0000665 3.596e-05 4.2977 3.0445 1.8700 1.1826 1.10560 0.94049
-Proportion Explained 0.000120 0.00006 0.0000200 0.0000100 0.0000000 0.000e+00 0.2388 0.1691 0.1039 0.0657 0.06142 0.05225
-Cumulative Proportion 0.128830 0.12889 0.1289000 0.1289100 0.1289200 1.289e-01 0.3677 0.5368 0.6407 0.7064 0.76783 0.82008
-PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18
-Eigenvalue 0.70216 0.61116 0.50930 0.46537 0.32289 0.24648 0.1782 0.07293 0.05366 0.04986 0.01738 0.009213
-Proportion Explained 0.03901 0.03395 0.02829 0.02585 0.01794 0.01369 0.0099 0.00405 0.00298 0.00277 0.00097 0.000510
-Cumulative Proportion 0.85909 0.89304 0.92134 0.94719 0.96513 0.97882 0.9887 0.99277 0.99575 0.99852 0.99949 1.000000'
-#anova(rda1, perm=99)
-#anova(rda1, by="terms", permutations = how(nperm=9))
-library(parallel)
-detectCores() #32
-# Set up the cluster
-#clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
-#clust <- try(makeCluster(getOption("cl.cores", 2), type = clusterType))
-cl <- makeCluster(10)
-#registerDoParallel(cl)
-options(na.action = "na.pass")
-#clusterExport(cl, "results")
-
-#print(system.time(anova(rda1, by="terms", perm=99, parallel = getOption("mc.cores"))))
-print(system.time(out1 <- anova(rda1, by="margin", permutations =99, parallel = cl)))
-'user system elapsed
-5304.802 1311.088 14522.410
-'
-out1
-
-traitnames <- dimnames(summary(rda1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda1)$species[,c(1,2)]
-bioclimcoord <- summary(rda1)$biplot[,c(1,2)]
-which(traitnames=="LeafNperArea") #9
-traitcoord[9,1] <- traitcoord[9,1]+0.5
-which(traitnames=="StemDens") #5
-traitcoord[5,2] <- traitcoord[5,2]+0.3
-which(traitnames=="Disp.unit.leng") #18
-traitcoord[18,2] <- traitcoord[18,2]-0.3
-ordiplot(rda1, type="n", xlim=c(-7,7))
-arrows(rep(0,45),rep(0,45),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="cadetblue4", length=0.1)
-text(bioclimcoord*8,dimnames(bioclimcoord)[[1]], col="cadetblue4", cex=0.8)
-text(traitcoord,traitnames, col="red")
-# copied as first RDA graph into ppt
-
-#rda1 <- rda(CWM ~ ., splot.bioclim[index11,c(4:49)],scale=T)
-
-### Stepwise RDA ###
-rda2 <- rda(CWM ~ 1,splot.bioclim[index11,c(4:49)],scale=T)
-#rda3 <- step(rda2, scope = formula(rda1), test = "perm")
-rda3 <- ordiR2step(rda2, scope=formula(rda1), direction = "both", Pin = 0.05, R2scope = F, permutations = how(nperm = 99), trace = TRUE)
-#TRY2.0: CWM ~ GDD5 + BIO_02 + BIO_12 + T_MAY + P_MAY + T_AUG ....
-#TRY3.0: CWM ~ GDD5 + BIO_02 + BIO_12 + T_MAY + P_MAY + T_SEP + P_DEC + P_FEB
-#rm(rda1)
-
-# now rebuild the RDA sequentially
-pca1 <- rda(CWM,scale=T)
-pca1
-5.527+ 3.906+ 1.953+ 1.507+ 1.375+ 0.886+ 0.726+ 0.623+ 0.499 +0.360+ 0.295+ 0.131+ 0.087+ 0.077+ 0.048
-(5.527+ 3.906)/18 # 52.40556%
-
-traitnames <- dimnames(summary(pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(pca1)$species[,c(1,2)]
-which(traitnames=="SeedMass") #3
-traitcoord[3,2] <- traitcoord[3,2]-0.5
-ordiplot(pca1, type="n")
-text(traitcoord,traitnames, col="red")
-# copied as second RDA graph into ppt
-
-
-rda4 <- rda(CWM ~ GDD5,splot.bioclim[index11,c(4:49)],scale=T)
-rda4
-traitnames <- dimnames(summary(rda4)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda4)$species[,c(1,2)]
-plotcoord <- summary(rda4)$sites[,c(1,2)]
-str(plotcoord) #num [1:1112287, 1:2]
-bioclimcoord <- summary(rda4)$biplot[,c(1,2)]
-
-min(plotcoord[,1]) #-0.683271
-max(plotcoord[,1]) # 0.8519885
-min(plotcoord[,2]) #-0.3398233
-max(plotcoord[,2]) # 0.3047364
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #153 64
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1112287
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30)) #T
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.947336
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-
-ordiplot(rda4, type="n")
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-
-arrows(rep(0,1),rep(0,1),bioclimcoord[1]*8,bioclimcoord[2]*8, col="blue", length=0.1)
-text(bioclimcoord[1]*8+0.8,bioclimcoord[2]*8-0.2,"GDD5", col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as third RDA graph into ppt
-
-rda5 <- rda(CWM ~ GDD5+BIO_02,splot.bioclim[index11,c(4:49)],scale=T)
-rda5
-traitnames <- dimnames(summary(rda5)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(rda5)$species[,c(1,2)]
-bioclimcoord <- summary(rda5)$biplot[,c(1,2)]
-plotcoord <- summary(rda5)$sites[,c(1,2)]
-str(plotcoord) #num [1:1112287, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #141 205
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.228913
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-
-
-ordiplot(rda5, type="n", xlim=c(-10,10), ylim=c(-10,10))
-points(plotcoord*15, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,2),rep(0,2),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="blue", length=0.1)
-text(bioclimcoord*8,dimnames(bioclimcoord)[[1]], col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as fourth RDA graph into ppt
-
-
-plot(CWM[,"LeafNPratio.mean"]~splot.bioclim$GDD5[index11],
- cex=0.05, xlab="GDD5", ylab="CWM ln(leaf N to P ratio)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"LeafNPratio.mean"]~splot.bioclim$GDD5[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.334e+00 3.795e-04 6149.5 <2e-16 ***
-splot.bioclim$GDD5[index11] 5.757e-06 1.595e-08 360.9 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.196 on 1112285 degrees of freedom
-Multiple R-squared: 0.1048, Adjusted R-squared: 0.1048
-F-statistic: 1.303e+05 on 1 and 1112285 DF, p-value: < 2.2e-16
-'
-abline(model1, lwd=3, col="blue")
-# LeafNPRatio_GDD5
-
-plot(CWM[,"SLA.mean"]~jitter(splot.bioclim$BIO_02[index11]),
- cex=0.05, xlab="BIO_02", ylab="CWM ln(SLA)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"SLA.mean"]~splot.bioclim$BIO_02[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 3.399e+00 1.448e-03 2348.2 <2e-16 ***
-splot.bioclim$BIO_02[index11] -6.253e-03 1.551e-05 -403.3 <2e-16 ***
-
-Residual standard error: 0.4034 on 1111305 degrees of freedom
-Multiple R-squared: 0.1277, Adjusted R-squared: 0.1277
-'
-abline(model1, lwd=3, col="blue")
-# SLA_BIO_02
-
-plot(CWM[,"LeafNperArea.mean"]~jitter(splot.bioclim$BIO_02[index11]),
- cex=0.05, xlab="BIO_02", ylab="CWM ln(Leaf N per area)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(CWM[,"LeafNperArea.mean"]~splot.bioclim$BIO_02[index11])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) -5.989e-02 8.475e-04 -70.67 <2e-16 ***
-splot.bioclim$BIO_02[index11] 3.696e-03 9.079e-06 407.09 <2e-16 ***
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Residual standard error: 0.2362 on 1111305 degrees of freedom
-Multiple R-squared: 0.1298, Adjusted R-squared: 0.1298
-'
-abline(model1, lwd=3, col="blue")
-# LeafNper_BIO_02
-
-##### FD calculations ####
-j <- 3
-nam <- DT$species[DT$PlotObservationID==j]
-abu <- DT$Relative.cover[DT$PlotObservationID==j]
-trait <- DT$trait[DT$PlotObservationID==j]
-abundance.weighted=T
-
-FD.fun <- function(trait, abu){
- res <- as.double(NA)
- if (length(trait[!is.na(trait)])>0){
- res <- 0
- #nam <- as.character(nam[!is.na(trait)])
- abu <- as.data.frame(abu[!is.na(trait)])
- #rownames(abu) <- nam
- trait <- as.data.frame(trait[!is.na(trait)])
- #rownames(trait) <- nam
- dis <- dist(trait,method="euclidean")
- if (sum(dis)>1){
- sqrt.dis <- sqrt(dis)
- # taking the square root is required to obtain the same result from the divc function
- # as divc takes the square of distances values, as suggested by
- # Champely & Chessel (2002, Env. Ecol. Stat. 9: 167-177)
- #abu <- as.data.frame(abu)
- res <- unlist(divc(abu, sqrt.dis))
- #Rao's Q: sumi sumj (dist ij prop i prop j)'
- }
- }
- res
-}
-
-mpd.fun <- function(nam, trait, abu,abundance.weighted){
- #res <- numeric(1) # carries result from the function
- res <- as.double(NA)
- if (length(trait[!is.na(trait)])>0){
- # dis <- as.matrix(dist(t(t1),method="euclidean"))
- abu <- t(abu)
- colnames(abu) <- nam
- trait <- t(trait)
- colnames(trait) <- nam
- #dis <- as.matrix(vegdist(t(t1),method="euclidean", na.rm=T))
- dis <- as.matrix(dist(t(trait),method="euclidean"))
- res <- mpd(abu, dis, abundance.weighted=abundance.weighted)
- }
- res
-}
-
-# from mpd source code
-'N <- dim(samp)[1]
-mpd <- numeric(N)
-for (i in 1:N) {
- sppInSample <- names(samp[i, samp[i, ] > 0])
- if (length(sppInSample) > 1) {
- sample.dis <- dis[sppInSample, sppInSample]
- if (abundance.weighted) {
- sample.weights <- t(as.matrix(samp[i, sppInSample,
- drop = FALSE])) %*% as.matrix(samp[i, sppInSample,
- drop = FALSE])
- mpd[i] <- weighted.mean(sample.dis, sample.weights)
- }
- else {
- mpd[i] <- mean(sample.dis[lower.tri(sample.dis)])
- }
- }
- else {
- mpd[i] <- NA
- }
-}
-mpd
-'
-# from Oliver, who used picante
-# does not accept missing values in the function
-library(picante)
-# define function for mean pairwise distance
-'mpd.fun.dt <- function(abu, dis, nam, abundance.weighted){
- abu <- t(abu)
- colnames(abu) <- nam
- res <- mpd(abu, dis, abundance.weighted=abundance.weighted)
-}
-# define funtion of mean nearest neighbor distance
-mntd.fun.dt <- function(abu, dis, nam, abundance.weighted){
- abu <- t(abu)
- colnames(abu) <- nam
- res <- mntd(abu, dis, abundance.weighted=abundance.weighted)
- res
-}
-species.list <- unique(DT$species)
-any(is.na(species.list)) #T, there is one NA
-length(species.list) #60909
-species.list <- species.list[!is.na(species.list)]
-length(species.list) #60908
-index8 <- match(species.list,TRY.all.mean.sd.2$species)
-str(index8) # 60908
-length(index8[!is.na(index8)]) #21001
-#species.list2 <- TRY.all.mean.sd.2$species[index8]
-#species.list2 <- species.list2[!is.na(species.list2)]
-#length(species.list2) # 21001
-TRY.all.mean.sd.2.na <- TRY.all.mean.sd.2[index8,]
-str(TRY.all.mean.sd.2.na)
-TRY.all.mean.sd.2.na$species2 <- species.list
-rm(index9)
-'
-FD <- array(NA,c(dim(CWM)[1],18),dimnames=list(dimnames(CWM)[[1]],
- colnames(TRY.all.mean.sd.2)[3:20]))
-str(FD) #logi [1:1117898, 1:18]
-'for (i in 1:18){
- # create distance matrix
- #trait.vec <- TRY.all.mean.sd.2.na[,i+2]
- #names(trait.vec) <- TRY.all.mean.sd.2.na$species2
- #trait.dist.mat <- as.matrix(dist(trait.vec))
- # too large for my notebook
- #mpd[,i] <- DT[,list(mpd.trait= mpd.fun.dt(Relative.cover, trait.dist.mat, species, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
- print(i)
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.2[index7,i+2]
- FD[,i] <- DT[,list(mpd.trait= mpd.fun(species, trait, Relative.cover, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
-}
-write.csv(FD,file = "FD2.csv", row.names = T)
-'
-for (i in 1:18){
- # create distance matrix
- #trait.vec <- TRY.all.mean.sd.2.na[,i+2]
- #names(trait.vec) <- TRY.all.mean.sd.2.na$species2
- #trait.dist.mat <- as.matrix(dist(trait.vec))
- # too large for my notebook
- #mpd[,i] <- DT[,list(mpd.trait= mpd.fun.dt(Relative.cover, trait.dist.mat, species, abundance.weighted=TRUE)),by=PlotObservationID]$mpd.trait
- print(i)
- DT$trait <- NA
- DT$trait <- TRY.all.mean.sd.2[index7,i+2]
- # recalculate relative cover for this trait
- FD[,i] <- DT[,list(FD.trait= FD.fun(trait, Relative.cover)),by=PlotObservationID]$FD.trait
-}
-write.csv(FD,file = "FD3.csv", row.names = T)
-
-#FD <- read.csv("C://Daten//iDiv2//splot2//FD3.csv")
-FD <- read.csv("C://Daten//iDiv2//splot2//splot2.0//FD_TRY3.csv")
-str(FD)
-dimnames(FD)[[1]] <- FD$X
-FD <- FD[,-1]
-tail(FD)
-
-#### match FD to splot.header ####
-index12 <- match(dimnames(FD)[[1]],splot.header$PlotObservationID)
-length(index12) #1117898
-any(is.na(index12)) # F
-splot.header3 <- splot.header[index12,]
-str(FD) # num [1:1117898, 1:18]
-FD <- FD[!(is.na(splot.header3$Longitude)|is.na(splot.header3$Latitude)|is.na(FD[,"SLA.mean"])),]
-str(FD) # num [1:1112287 , 1:18]
-1112287 /1117898*100 # 99.49808 % of plots have traits, Long and Lat
-any(is.na(FD)) #F
-# match header again
-index13 <- match(dimnames(FD)[[1]],splot.header$PlotObservationID)
-length(index13) #1112287
-any(is.na(index13)) # F
-splot.header3 <- splot.header[index13,]
-# now in the same sequence as FD
-
-
-#### Collect all FD regressions ####
-names(splot.bioclim) # from 4 BIO_01 to 49 T_MAR
-index14 <- match(splot.header3$PlotObservationID,splot.bioclim$PLOT_ID)
-length(index14) # 1111307
-length(index14[!is.na(index14)]) # 1111307
-
-FD.bioclim <- array(NA,c(46,18,5),dimnames=list(names(splot.bioclim)[4:49],dimnames(FD)[[2]],
- c("p_lin","r2_lin","p_qua1", "p_qua2","r2_qua")))
-for (i in 1:46){
- for (j in 1:18){
- model1 <- lm(FD[,j]~splot.bioclim[index14,i+3])
- FD.bioclim[i,j,1] <- summary(model1)$coefficients[2,4]
- FD.bioclim[i,j,2] <- summary(model1)$r.squared
- model2 <- lm(FD[,j]~splot.bioclim[index14,i+3]+I(splot.bioclim[index14,i+3]^2))
- FD.bioclim[i,j,3] <- summary(model2)$coefficients[2,4]
- FD.bioclim[i,j,4] <- summary(model2)$coefficients[3,4]
- FD.bioclim[i,j,5] <- summary(model2)$r.squared
- }
-}
-str(FD.bioclim) #num [1:46, 1:18, 1:5]
-write.csv(data.frame(names(splot.bioclim)[4:49],FD.bioclim[,,2]),file = "FD2_bioclim_r2_lin.csv", row.names = FALSE)
-write.csv(data.frame(names(splot.bioclim)[4:49],FD.bioclim[,,5]),file = "FD2_bioclim_r2_qua.csv", row.names = FALSE)
-
-library(vegan)
-FD.pca1 <- rda(FD,scale=T)
-FD.pca1
-'
-Eigenvalues for unconstrained axes:
- PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.703 1.758 1.122 0.974 0.732 0.643 0.590 0.540
-(Showed only 8 of all 18 unconstrained eigenvalues)
-'
-plot(FD.pca1)
-head(summary(FD.pca1))
-#text(summary(FD.pca1)$species[,c(1,2)],dimnames(summary(FD.pca1)$species)[[1]])
-traitnames <- dimnames(summary(FD.pca1)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.pca1)$species[,c(1,2)]
-which(traitnames=="Seed.length") #17
-traitcoord[17,2] <- traitcoord[17,2]+0.3
-which(traitnames=="LeafC.perdrymass") #12
-traitcoord[12,2] <- traitcoord[12,2]-0.4
-which(traitnames=="Leaf.delta.15N") #13
-traitcoord[13,2] <- traitcoord[13,2]+0.4
-which(traitnames=="LeafNPratio") #11
-traitcoord[11,2] <- traitcoord[11,2]-0.3
-
-plotcoord <- summary(FD.pca1)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.pca1, type="n")
-#text(traitcoord,traitnames)
-points(plotcoord*30, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-
-#arrows(rep(0,45),rep(0,45),bioclimcoord[,1]*8,bioclimcoord[,2]*8, col="cadetblue4", length=0.1)
-text(traitcoord,traitnames, col="red")
-# copied as first RDA graph into ppt
-
-env1 <- envfit(FD.pca1,splot.bioclim[index14,c(4:49)],permutations = 0)
-env1
-#with permutations = 999:
-#Error: cannot allocate vector of size 4.1 Gb
-plot(env1)
-png(filename="PCA1.png", width=1000, height=1000, units = "px", bg="white")
-ordiplot(pca1, type="n")
-text(traitcoord,traitnames)
-plot(env1)
-dev.off() # creates a file
-
-### Stepwise RDA ###
-### RDA ###
-FD.rda1 <- rda(FD~.,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda1
-FD.rda2 <- rda(FD ~ 1,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda3 <- ordistep(FD.rda2, scope=formula(FD.rda1), direction = "both", Pin = 0.1,
- permutations = how(nperm = 10), trace = TRUE)
-'Start: FD ~ 1
-
-Df AIC F Pr(>F)
-+ BIO_14 1 3191679 20602.8 0.09091 .
-+ BIO_17 1 3191898 20379.5 0.09091 .
-+ BIO_12 1 3197070 15125.2 0.09091 .
-+ P_MAY 1 3198071 14110.8 0.09091 .
-+ P_APR 1 3198194 13986.4 0.09091 .
-+ BIO_18 1 3198767 13406.3 0.09091 .
-+ P_MAR 1 3199448 12717.1 0.09091 .
-+ P_FEB 1 3199451 12714.5 0.09091 .
-+ P_JUN 1 3199643 12519.9 0.09091 .
-+ BIO_05 1 3199804 12357.7 0.09091 .
-+ BIO_02 1 3200062 12096.5 0.09091 .
-+ P_NOV 1 3200602 11551.1 0.09091 .
-+ P_DEC 1 3200823 11328.1 0.09091 .
-+ P_JAN 1 3201335 10810.2 0.09091 .
-+ BIO_10 1 3202194 9943.1 0.09091 .
-+ T_JUL 1 3202226 9911.4 0.09091 .
-+ T_AUG 1 3202685 9448.0 0.09091 .
-+ GDD5 1 3202914 9217.7 0.09091 .
-+ T_JUN 1 3202934 9196.7 0.09091 .
-+ BIO_16 1 3203459 8668.3 0.09091 .
-+ BIO_19 1 3203685 8440.1 0.09091 .
-+ P_AUG 1 3204079 8043.8 0.09091 .
-+ BIO_13 1 3204383 7737.1 0.09091 .
-+ T_MAY 1 3204450 7669.7 0.09091 .
-+ P_OCT 1 3204755 7362.4 0.09091 .
-+ T_SEP 1 3204986 7130.0 0.09091 .
-+ P_SEP 1 3205174 6941.4 0.09091 .
-+ BIO_15 1 3205378 6735.8 0.09091 .
-+ P_JUL 1 3205600 6512.4 0.09091 .
-+ T_APR 1 3205746 6365.1 0.09091 .
-+ BIO_01 1 3205841 6270.3 0.09091 .
-+ BIO_03 1 3205882 6228.8 0.09091 .
-+ T_OCT 1 3206048 6062.0 0.09091 .
-+ AR 1 3206562 5545.3 0.09091 .
-+ BIO_07 1 3207063 5042.2 0.09091 .
-+ T_NOV 1 3207547 4555.5 0.09091 .
-+ T_MAR 1 3207659 4443.2 0.09091 .
-+ T_FEB 1 3208177 3923.3 0.09091 .
-+ T_DEC 1 3208260 3840.3 0.09091 .
-+ BIO_11 1 3208334 3765.7 0.09091 .
-+ T_JAN 1 3208465 3634.3 0.09091 .
-+ BIO_08 1 3208844 3254.4 0.09091 .
-+ BIO_06 1 3209374 2723.0 0.09091 .
-+ BIO_09 1 3209543 2553.3 0.09091 .
-+ PET 1 3209642 2454.5 0.09091 .
-+ BIO_04 1 3210074 2021.3 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Step: FD ~ BIO_14
-
-Df AIC F Pr(>F)
-- BIO_14 1 3212091 20603 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Df AIC F Pr(>F)
-+ BIO_02 1 3184224 7481.80 0.09091 .
-+ BIO_03 1 3184435 7269.37 0.09091 .
-+ GDD5 1 3187255 4435.29 0.09091 .
-+ BIO_05 1 3187772 3915.70 0.09091 .
-+ T_JAN 1 3188026 3661.00 0.09091 .
-+ BIO_01 1 3188040 3647.18 0.09091 .
-+ T_FEB 1 3188048 3638.72 0.09091 .
-+ T_MAR 1 3188075 3612.29 0.09091 .
-+ T_APR 1 3188082 3605.14 0.09091 .
-+ BIO_11 1 3188088 3598.96 0.09091 .
-+ T_DEC 1 3188123 3563.70 0.09091 .
-+ T_NOV 1 3188208 3478.25 0.09091 .
-+ T_OCT 1 3188225 3461.84 0.09091 .
-+ T_MAY 1 3188357 3329.42 0.09091 .
-+ BIO_06 1 3188491 3194.91 0.09091 .
-+ BIO_10 1 3188609 3076.07 0.09091 .
-+ T_SEP 1 3188617 3068.49 0.09091 .
-+ BIO_13 1 3188664 3021.74 0.09091 .
-+ T_JUN 1 3188669 3016.26 0.09091 .
-+ BIO_16 1 3188709 2976.11 0.09091 .
-+ P_SEP 1 3188724 2961.09 0.09091 .
-+ BIO_07 1 3188728 2956.84 0.09091 .
-+ BIO_12 1 3188782 2902.48 0.09091 .
-+ T_AUG 1 3188788 2897.30 0.09091 .
-+ T_JUL 1 3188879 2805.50 0.09091 .
-+ P_AUG 1 3189181 2503.15 0.09091 .
-+ BIO_18 1 3189256 2427.84 0.09091 .
-+ BIO_04 1 3189277 2406.93 0.09091 .
-+ P_JUL 1 3189398 2285.76 0.09091 .
-+ BIO_15 1 3189605 2078.03 0.09091 .
-+ P_APR 1 3189650 2032.61 0.09091 .
-+ P_JUN 1 3189703 1979.86 0.09091 .
-+ BIO_09 1 3189912 1770.71 0.09091 .
-+ P_MAR 1 3189965 1717.19 0.09091 .
-+ PET 1 3189984 1698.34 0.09091 .
-+ AR 1 3190003 1678.99 0.09091 .
-+ P_MAY 1 3190005 1677.55 0.09091 .
-+ P_OCT 1 3190063 1619.65 0.09091 .
-+ P_FEB 1 3190490 1192.28 0.09091 .
-+ BIO_17 1 3190497 1184.39 0.09091 .
-+ BIO_19 1 3190555 1126.59 0.09091 .
-+ P_NOV 1 3190570 1111.68 0.09091 .
-+ BIO_08 1 3190582 1100.05 0.09091 .
-+ P_JAN 1 3190893 787.98 0.09091 .
-+ P_DEC 1 3190955 726.52 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-
-Step: FD ~ BIO_14 + BIO_02
-
-Df AIC F Pr(>F)
-- BIO_02 1 3191679 7481.8 0.09091 .
-- BIO_14 1 3200062 15953.2 0.09091 .
----
-Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1'
-
-FD.rda4 <- rda(FD ~ BIO_14,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda4
-' Inertia Proportion Rank
-Total 18.0000 1.0000
-Constrained 0.3276 0.0182 1
-Unconstrained 17.6724 0.9818 18
-Inertia is correlations
-
-Eigenvalues for constrained axes:
-RDA1
-0.3276
-
-Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.307 1.850 1.378 0.884 0.789 0.780 0.578 0.551
-'
-traitnames <- dimnames(summary(FD.rda4)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.rda4)$species[,c(1,2)]
-bioclimcoord <- summary(FD.rda4)$biplot[,c(1,2)]
-plotcoord <- summary(FD.rda4)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.rda4, type="n", xlim=c(-10,10),ylim=c(-10,12))
-points(plotcoord*15, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,1),rep(0,1),bioclimcoord[1]*8,bioclimcoord[2]*8, col="blue", length=0.1)
-text(bioclimcoord[1]*8,bioclimcoord[2]*8+0.5,"BIO_14", col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as third RDA graph into ppt
-
-FD.rda5 <- rda(FD ~ BIO_14+BIO_02,splot.bioclim[index14,c(4:49)],scale=T)
-FD.rda5
-'
-Inertia Proportion Rank
-Total 18.00000 1.00000
-Constrained 0.44582 0.02477 2
-Unconstrained 17.55418 0.97523 18
-Inertia is correlations
-
-Eigenvalues for constrained axes:
-RDA1 RDA2
-0.3641 0.0817
-
-Eigenvalues for unconstrained axes:
-PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
-8.267 1.816 1.377 0.880 0.788 0.780 0.568 0.530 '
-traitnames <- dimnames(summary(FD.rda5)$species)[[1]]
-traitnames <- substr(traitnames,1,nchar(traitnames)-5)
-traitcoord <- summary(FD.rda5)$species[,c(1,2)]
-bioclimcoord <- summary(FD.rda5)$biplot[,c(1,2)]
-
-plotcoord <- summary(FD.rda5)$sites[,c(1,2)]
-str(plotcoord) #num [1:1111307, 1:2]
-x.cut <- cut(plotcoord[,1], breaks=seq(min(plotcoord[,1]), max(plotcoord[,1]), by=0.01))
-y.cut <- cut(plotcoord[,2], breaks=seq(min(plotcoord[,2]), max(plotcoord[,2]), by=0.01))
-z.cut <- table(x.cut,y.cut)
-head(z.cut)
-dim(z.cut) #131 50
-z.cut <- melt(z.cut)
-str(z.cut) # 6550 obs. of 3 variables:
-head(z.cut)
-z.cut$xy <- paste(z.cut$x.cut,z.cut$y.cut,sep="_")
-str(x.cut)
-length(x.cut) # 1111307
-xy <- paste(x.cut,y.cut,sep="_")
-index30 <- match(xy,z.cut$xy)
-length(index30)
-any(is.na(index30))
-head(log10(z.cut$value[index30]))
-values.present <- vector(mode="numeric", length=length(xy))
-min(log10(z.cut$value[index30]), na.rm=T) #0
-max(log10(z.cut$value[index30]), na.rm=T) # 3.399154
-values.cut <- cut(log10(z.cut$value[index30]), breaks=seq(min(log10(z.cut$value[index30]), na.rm=T),
- 4, by=0.4))
-
-ordiplot(FD.rda5, type="n", xlim=c(-10,10),ylim=c(-10,12))
-points(plotcoord*15, cex=0.05, col=rev(greyPalette(n=10, start = 0, end = 0.8,
- gamma = 2.2, alpha = NULL))[values.cut])
-arrows(rep(0,2),rep(0,2),bioclimcoord[,1]*4,bioclimcoord[,2]*4, col="blue", length=0.1)
-text(bioclimcoord*4,dimnames(bioclimcoord)[[1]], col="blue", cex=1)
-text(traitcoord,traitnames, col="red")
-# copied as fourth RDA graph into ppt
-
-plot(FD[,"Seed.length.mean"]~splot.bioclim$BIO_14[index14],
- cex=0.05, xlab="BIO_14", ylab="FD ln(seed length)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(FD[,"Seed.length.mean"]~splot.bioclim$BIO_14[index14])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 2.584e-01 3.193e-04 809.1 <2e-16 ***
-splot.bioclim$BIO_14[index14] 1.407e-04 6.813e-07 206.5 <2e-16 ***
-Multiple R-squared: 0.03697, Adjusted R-squared: 0.03697
-'
-abline(model1, lwd=3, col="blue")
-#
-
-plot(FD[,"Seed.num.rep.unit.mean"]~splot.bioclim$BIO_02[index14],
- cex=0.05, xlab="BIO_02", ylab="FD ln(Seed number of reproductive unit)",cex.lab=1.5, cex.axis=1.5)
-model1 <- lm(FD[,"Seed.num.rep.unit.mean"]~splot.bioclim$BIO_02[index14])
-summary(model1)
-' Estimate Std. Error t value Pr(>|t|)
-(Intercept) 1.251e+00 1.799e-03 695.3 <2e-16 ***
-splot.bioclim$BIO_02[index14] -3.751e-03 1.927e-05 -194.6 <2e-16 ***
-Multiple R-squared: 0.03296, Adjusted R-squared: 0.03296 '
-abline(model1, lwd=3, col="blue")
-
-
-
-
-
-
-
-
-
-
-
-
-# mpd presence-absence
-system.time(
- mpd.sla.no.abu <- DT[,list(mpd.sla = mntd.fun.dt(Cover.., SLA.dist.mat, StandSpeciesName, abundance.weighted=FALSE)),by=PlotObservationID]
-)
-
-
-###########
-levels(splot.header2$COUNTRY_CO)
-levels(splot.header$Country)
-CWM2$Country <- splot.header$Country[index12]
-data(countryExData)
-str(countryExData)
-sPDF <- joinCountryData2Map(splot.header, joinCode = "NAME",nameJoinColumn = "Country")
-?mapCountryData
-#mapCountryData(sPDF, nameColumnToPlot="BIODIVERSITY" )
-
-#### Climatic space, from Jonathon Lenoir ####
-library(fBasics)
-# Color palette
-
-# Import bio1 (MAT) of the world
-bio <- read.table(paste("/home/helge/sPlot/Bioclim/", "BIO.txt", sep = ""),sep=";", header=T)
-str(bio) #9033453 obs. of 22 variables
-r <- raster(nrows=33, ncols=33, xmn=min(bio$BIO_01), xmx=max(bio$BIO_01), ymn=min(bio$BIO_12), ymx=max(bio$BIO_12))
-mattar <- data.frame(bio$BIO_01,bio$BIO_12)
-str(mattar)
-names(mattar) <- c("bio1","bio12")
-mattar_r <- rasterize(mattar, r, fun="count")
-str(mattar_r)
-plot(log10(mattar_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="N 2.5 arc-minutes cells per bin (log10 scale)", cex.main=0.7,
- xlab="MAT (°C*10)", ylab="TAR (mm)")
-# Bio1_Bio12_world.png
-
-#plot(mattar_r/max(getValues(mattar_r), na.rm=TRUE), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))), main="Relative proportion of 10 arc-minutes cells per bin", xlab="MAR (°C*10)", ylab="TAR (mm)")
-# Plot the sampling effort of sPlot across the MAT-TAR climatic space
-str(splot.bioclim)
-# Make it a spatial data frame
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-splot.bioclim.coords <- cbind(splot.bioclim$X, splot.bioclim$Y)
-splot.bioclim.coords <- SpatialPoints(splot.bioclim.coords, proj4string=CRSlonlat)
-splot.bioclim2 <- data.frame(splot.bioclim$BIO_01,splot.bioclim$BIO_12)
-names(splot.bioclim2) <- c("bio1","bio12")
-#splot.bioclim2 <- SpatialPointsDataFrame(splot.bioclim.coords, splot.bioclim2, proj4string=CRSlonlat)
-#str(splot.bioclim2)
-#sPlot_mattar <- extract(bio1, coordinates(sPlot))
-#sPlot_mattar <- as.data.frame(sPlot_mattar)
-#names(sPlot_mattar) <- c("bio1")
-#sPlot_mattar$bio12 <- extract(bio12, coordinates(sPlot))
-sPlot_mattar_r <- rasterize(splot.bioclim2, r, fun="count")
-str(sPlot_mattar_r)
-plot(mattar_r, asp=0, col="grey", legend=FALSE, main="N vegetation plots per bin (log10 scale)",
- xlab="MAT (°C)", ylab="TAR (mm)")
-# Bio1_Bio12_world_grey.png
-plot(log10(sPlot_mattar_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- add=TRUE)
-# Bio1_Bio12_sPlot_SLA.png
-
-# Build the PC1-PC2 climatic space
-load("/home/helge/sPlot/Bioclim/pca_data.RData")
-str(pca.data)
-r <- raster(nrows=33, ncols=33, xmn=min(pca.data$PC1), xmx=max(pca.data$PC1), ymn=min(pca.data$PC2), ymx=max(pca.data$PC2))
-pca_r <- rasterize(pca.data[, c("PC1", "PC2")], r, fun="count")
-plot(log10(pca_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="N 2.5 arc-minutes cells per bin (log10 scale)",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-plot(pca_r/max(getValues(pca_r), na.rm=TRUE), asp=0,
- col=rev(divPalette(n=100, name=c("RdYlBu"))),
- main="Relative proportion of 2.5 arc-minutes cells per bin",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-# Plot the sampling effort of sPlot across the PC1-PC2 climatic space
-index13 <- match(splot.bioclim$RAST_ID,pca.data$RAST_ID)
-sPlot_pca <- pca.data[index13, c("PC1", "PC2")]
-sPlot_pca_r <- rasterize(sPlot_pca, r, fun="count")
-plot(pca_r, asp=0, col="grey", legend=FALSE,
- main="N vegetation plots per bin (log10 scale)",
- xlab="PC1 (cold-to-warm)", ylab="PC2 (dry-to-wet)")
-plot(log10(sPlot_pca_r), asp=0, col=rev(divPalette(n=100, name=c("RdYlBu"))),
- add=TRUE)
-
-#### For Marten Winter ####
-# species x country Matrix
-species.country <- table(DT$species,splot.header$Country[index30])
-str(species.country) # int [1:60908, 1:130]
-names(dimnames(species.country))
-names(dimnames(species.country)) <- c("species","country")
-
-library(reshape2)
-?melt.table
-species.country2 <- melt(species.country)
-str(species.country2) # 188470 obs. of 3 variables:
-names(species.country2)
-names(species.country2)[3] <- "p_a"
-table(species.country2$p_a)
-
-species.country2$p_a[species.country2$p_a>0] <- 1
-species.country2 <- species.country2[species.country2$p_a==1,]
-str(species.country2) # 188470 obs. of 3 variables:
-species.country2 <- species.country2[,-3]
-str(species.country2) # 188470
-write.csv(species.country2,file =
- paste("/home/oliver/shared/", "species.country2.csv", sep = ""), row.names = FALSE)
-levels(species.country2$country) # 130
-length(species.country2$country[species.country2$country=="Algeria"])
-# 97
-sum(species.country[,2])
-# 799, this is the amount of records for DZA = Algeria
-sum(species.country[,126])
-# 1467200, sum of records for United States
-
-TaxaNat <- read.csv(paste("/home/oliver/shared/Glonaf/", "TaxaNat.csv", sep = ""))
-str(TaxaNat) #178779 obs. of 2
-
-GlonafRegion <- read.csv(paste("/home/oliver/shared/Glonaf/", "GlonafRegion.csv", sep = ""))
-str(GlonafRegion) #843 obs. of 10
-trim.trailing <- function (x) sub("\\s+$", "", x)
-# returns string w/o leading or trailing whitespace
-GlonafRegion$Country_sPlot <- as.factor(trim.trailing(as.character(GlonafRegion$Country_sPlot)))
-index14 <- match(levels(species.country2$country),levels(GlonafRegion$Country_sPlot))
-index14
-levels(GlonafRegion$Country)
-
-data.frame(levels(species.country2$country),levels(GlonafRegion$Country_sPlot)[index14])
-## all correct
-
-index15 <- match(levels(GlonafRegion$Country_sPlot),GlonafRegion$Country_sPlot)
-index15
-index16 <- match(species.country2$country,levels(species.country2$country))
-any(is.na(index16)) # F
-length(levels(species.country2$species)) #60908
-'species.country2$species.country <- paste(species.country2$species,
-levels(species.country2$country)[index16],sep="_")
-species.country2$species.country <- paste(species.country2$species,
-levels(GlonafRegion$Country_sPlot)[index14][index16],sep="_")
-species.country2$species.country <- paste(species.country2$species,
-GlonafRegion$Country_sPlot[index15][index14][index16],sep="_")
-'
-species.country2$species.country <- paste(species.country2$species,
- GlonafRegion$ISO.country[index15][index14][index16],sep="_")
-str(species.country2) # 188470 obs. of 3 variables:
-head(species.country2)
-species.country2[c(160000:160100),]
-## fits!
-
-index17 <- match(TaxaNat$Region,GlonafRegion$ID.region)
-any(is.na(index17)) # F
-
-TaxaNat$species.country <-paste(TaxaNat$Taxon2,GlonafRegion$ISO.country[index17], sep="_")
-
-index18 <- match(species.country2$species.country,TaxaNat$species.country)
-length(index18) #188470
-length(index18[!is.na(index18)]) # 9283 out of 188470 species x country occurrences
-# are alien occurrences
-9283/188470*100 # 4.925452 % species x country occurrences
-species.country2$alien <- F
-species.country2$alien[!is.na(index18)] <- T
-species.country2$ISO.country <- GlonafRegion$ISO.country[index15][index14][index16]
-str(species.country2)
-head(species.country2)
-length(species.country2$ISO.country[species.country2$ISO.country=="DZA"])
-# 12939, this cannot be
-#species.country2$ISO.country[species.country2$ISO.country=="DZA"]
-sum(species.country[,2])
-# 799, this is the amount of records for DZA = Algeria
-length(species.country2$ISO.country[species.country2$ISO.country=="DZA"
- & !is.na(species.country2$ISO.country)])
-# 97!
-
-#table(species.country2$ISO.country)
-#species.country3 <- acast(species.country2, species~ISO.country, length)
-length(species.country2$alien[species.country2$alien])
-# 9283 alien records in sPlot
-length(species.country2$alien[species.country2$alien & species.country2$ISO.country=="DZA"])
-# 3 alien records in DZA = Algeria
-length(species.country[species.country[,2]>0,2])
-# 97, however there are only 97 records in total
-length(species.country2$alien[species.country2$alien & species.country2$ISO.country=="USA"])
-# 1911 alien species in USA
-length(species.country[species.country[,126]>0,2])
-# 9933 species in total in USA
-alien.richness <- as.matrix(table(species.country2$ISO.country,species.country2$alien))[,2]
-
-species.country3 <- species.country
-species.country3[species.country3>0] <- 1
-country.richness <- colSums(species.country3)
-country.richness
-names(country.richness)
-data.frame(names(country.richness),GlonafRegion$ISO.country[index15][index14])
-country.richness2 <- data.frame(GlonafRegion$ISO.country[index15][index14],as.numeric(country.richness))
-names(country.richness2) <- c("ISO.country","richness")
-country.richness2$ISO.country
-country.richness2
-
-index19 <- match(country.richness2$ISO.country,names(alien.richness))
-index19
-alien.richness
-country.richness2$alien.richness <- alien.richness[index19]
-country.richness2
-country.richness2$alien.prop <- country.richness2$alien.richness/country.richness2$richness
-country.richness2$alien.prop
-
-sPDF <- joinCountryData2Map(country.richness2,joinCode = "ISO3",
- nameJoinColumn = "ISO.country")
-'112 codes from your data successfully matched countries in the map
-18 codes from your data failed to match with a country code in the map
-132 codes from the map werent represented in your data
-'
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="richness",addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="richness",
- catMethod=seq(0,10000,by=1000),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.richness",
- catMethod=seq(0,2000,by=200),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-max(country.richness2$alien.prop, na.rm=T) #0.2823779
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.prop",addLegend=FALSE)
-mapParams <- mapCountryData(sPDF, nameColumnToPlot="alien.prop",
- catMethod=seq(0,0.3,by=0.02),addLegend=FALSE)
-do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2,legendLabels="all"))
-
-
-
-
-##################### OLD ###################
-library(spdep)
-library(maptools)
-library(hexbin)
-library(RColorBrewer)
-# Use a lighter version of sPlot focusing on mean annual temperature MAT (BIO1)
-# and total annual rainfall TAR (BIO12)
-sPlot <- splot.header[, c("LONGITUDE", "LATITUDE", "BIO_1", "BIO_12")]
-sPlot$LONGITUDE <- as.numeric(as.character(sPlot$LONGITUDE))
-sPlot$LATITUDE <- as.numeric(as.character(sPlot$LATITUDE))
-sPlot$BIO_1 <- as.numeric(sPlot$BIO_1)/10
-sPlot$BIO_12 <- as.numeric(sPlot$BIO_12)
-str(sPlot)
-# Make it a spatial data frame
-CRSlonlat <- CRS("+proj=longlat +datum=WGS84")
-coords <- cbind(sPlot$LONGITUDE, sPlot$LATITUDE)
-coords <- SpatialPoints(coords, proj4string=CRSlonlat)
-sPlot <- SpatialPointsDataFrame(coords, sPlot, proj4string=CRSlonlat)
-str(sPlot)
-# Plot the data on the world map to check the coordinates and projection
-data(wrld_simpl)
-plot(wrld_simpl)
-points(coordinates(sPlot), cex=0.1, col="red", pch=16)
-# Set a color ramp
-rf <- colorRampPalette(rev(brewer.pal(11, 'Spectral')))
-# Plot the plot distribution within the MAT-TAR bidimentionnal climatic space
-hexbinplot(BIO_12~BIO_1, data=sPlot@data, colramp=rf, trans=log, inv=exp,
- border="white", xlab="MAT", ylab="TAR", xbins=30)
-# Display only bins with a minimum of 10 plots and a maximum of 1000 plots
-hexbinplot(BIO_12~BIO_1, data=sPlot@data, colramp=rf, trans=log, inv=exp,
- border="white", xlab="MAT", ylab="TAR", xbins=30, mincnt=10, maxcnt=1000)
-#
-
-
-
-
-
-
-
-
-################ OLD ###################
-
-###
-load("C:\\Daten\\iDiv\\sPlot\\pca_data.Rdata")
-str(pca.data)
-any(is.na(pca.data$PC1)) #F
-
-load("C:\\Daten\\iDiv\\sPlot\\sPlot_header_pca.RData")
-str(splot.header2)
-any(is.na(pca.data$PC1)) #F
-
-library(raster)
-r <-
- splot_pca_r <- rasterize(splot_pca, r, run="count")
-plot(splot_pca_r, asp=0, col=rev(divPalette(n=100, name="RdY1Bu"))))
-plot(splot_r, asp=0, col=rev(divPalette(n=100, name="RdY1Bu"))))
-
-index <- which(getValues(splot_pca_r)>0)
-which(getValues(pca_r)>0)
-poinra
-getValues(splot_pca_r)[index] # 136 values
-?sample
-splo_pca$binid <- cellFromXY(pca_r, splot_pca)
-
-
-###########################
-load("C:\\Daten\\iDiv\\sPlot\\splot_pca2.Rdata")
-str(splot_pca)
-splot_pca <- data.table(splot_pca)
-setkey(splot_pca,binid)
-inter <- splot_pca[,list(count=length(PCA1)), by=binid]
-splot_pca2 <- splot_pca[,list(PCA1, PCA2,count=length(PCA1)), by=binid]
-
-str(splot_pca2)
-hist(splot_pca2$binid)
-inter2 <- splot_pca2[count>100,list(count=sum(PCA1)), by=binid]
-str(inter2)
-[,list(total.cover=sum(Cover..)), by=PlotObservationID]
-<- DT[,list(total.cover=sum(Cover..)), by=PlotObservationID]
-
-sample(splot_pca$PCA1,100)
-splot_pca2 <- splot_pca[table(splot_pca$binid)>100,sample(PCA1,100),by=binid]
-str(splot_pca2)
-
-resample <- function(PCA1_,PCA2_,binid_){
- if (length(PCA1)>100) {
-
- }
- else {
- 1
- }
- ]index10 <- match(species,labels(d))
- d3 <- as.dist(dist[index10,index10])
- if (sum(d3)==0){
- NA
- }
- divc(as.data.frame(rel.cover),d3)$diversity
-}
-splot_pca2 <- splot_pca[length(PCA1)<=100,res_ID := 1]
-head(splot_pca2,100)
-table(splot_pca2$binid,splot_pca2$count, exclude=NULL)
-table(splot_pca$binid)
-splot_pca2 <- splot_pca[,res_ID := 0]
-splot_pca2 <- splot_pca[table(splot_pca2$binid)<=100,res_ID := 1]
-splot.selection <- splot_pca[,list(res = resample(PCA1,PCA2,binid)),by=binid]
-#length(splot_pca[binid>100])
-hist(splot_pca$binid)
-splot.selection <- splot_pca[,list(length(binid)),by=binid]
-splot.selection <- splot_pca[,list(count=length(PCA1)),by=binid]
-splot.selection <- splot_pca[,,by=binid]
-hist(splot.selection$count)
-str(splot.selection)
-str(splot_pca2)
-splot.selection <- splot_pca[,list(nrow()),by=binid]
-splot.selection <- splot_pca[,count:=length(binid),by=binid]
-length(splot_pca,by=binid)
-########################################################################################
-
-### CWM ###
-#CWM1 <- DT[!is.na(plot.species.trait.value),list(CWM = weighted.mean(plot.species.trait.value,Relative.cover)),by=PlotObservationID]
-CWM1 <- DT[,list(CWM = weighted.mean(plot.species.trait.value,Relative.cover)),by=PlotObservationID]
-# see http://stackoverflow.com/questions/3367190/aggregate-and-weighted-mean-in-r
-CWM1[1:100,]
-hist(CWM1$CWM)
-write.csv((x, file = paste(path.sPlot, "CWM_SLA_all.csv", sep = ""), row.names = FALSE))
-
-### FD ###
-library(ade4) # f?r divc (Rao's Q)
-# this does not work: ~50 GB RAM required for 40k x 40k species matrix:
-'ktab1 <- ktab.list.df(list(as.data.frame(species.trait.value[!is.na(species.trait.value)])))
-d <- dist.ktab(ktab1, type=c("Q"), option = "scaledBYrange")
-#labels(d)
-#str(d)
-#attr(d,"Labels")
-attr(d,"Labels") <- dimnames(species.trait.value)[[1]][!is.na(species.trait.value)]
-labels(d)
-any(is.na(d)) # F
-dist <- as.matrix(d)
-
-#tree.trait.weighted <- t(t(veg.KHG3.prop) * all.traits.KHG[,j])
-#dist.sqrt <- sqrt(dist)
-d.sqrt <- sqrt(d)
-labels(d.sqrt)
-dist <- as.matrix(d.sqrt)
-'
-'Rao <- function(rel.cover,species){
-index10 <- match(species,labels(d))
-d3 <- as.dist(dist[index10,index10])
-if (sum(d3)==0){
-NA
-}
-divc(as.data.frame(rel.cover),d3)$diversity
-}
-FD.all <- DT[!is.na(plot.species.trait.value),list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]
-'
-#species.trait.value.scaled <- species.trait.value/(max(species.trait.value)-min(species.trait.value))
-species.trait.value.scaled <- TRY.mean.sd.splot$SLA.mean/(max(TRY.mean.sd.splot$SLA.mean)-min(TRY.mean.sd.splot$SLA.mean))
-
-# scale by range, which saves time for distance calculations
-
-library(reshape2)
-
-package.size = 1000
-#number.plots.total <- length(levels(DT[!is.na(plot.species.trait.value)]$PlotObservationID)) # 705272
-number.plots.total <- length(levels(DT$PlotObservationID)) # 701611 vs. 705272
-FD.all <- array(0,number.plots.total, dimnames=list(levels(DT$PlotObservationID)))
-number.packages <-ceiling(number.plots.total/package.size) #702
-pb <- winProgressBar(title = "progress bar", min = 0, max =number.packages, width = 300)
-for (k in 1:number.packages) {
- setWinProgressBar(pb, k, title=paste("Row ",k,"of ", number.packages, " is processed"))
- ifelse (k==number.packages, package<-number.plots.total%%package.size, package <- package.size)
- #index3 <- match(DT[!is.na(plot.species.trait.value)]$PlotObservationID,levels(DT[!is.na(plot.species.trait.value)]$PlotObservationID)[seq((k-1)*package.size+1,k*package.size,1)])
- index3 <- match(DT$PlotObservationID,levels(DT$PlotObservationID)[seq((k-1)*package.size+1,k*package.size,1)])
- #DT2 <- DT[!is.na(plot.species.trait.value) & !is.na(index3)]
- #DT2 <- DT2[!is.na(plot.species.trait.value)]
- #length(index3[!is.na(index3)])
- DT2 <- DT[!is.na(index3)]
- #str(DT2)
- #veg.matrix1 <- acast(DT2, Matched.concept ~ PlotObservationID, value.var="Relative.cover", fill=0)
- veg.matrix1 <- acast(DT2, names.stand ~ PlotObservationID, value.var="Relative.cover", fun.aggregate=sum, fill=0)
- # fun.aggregate = sum -> sums up multiple entries per species!
- # this can happen with the unified species nomenclature!
- #length(unique(dimnames(veg.matrix1)[[2]]))# 1000
- #length(unique(DT2$PlotObservationID)) # 1000
- #veg.matrix1[veg.matrix1 > 0] <- 1
- #any(colSums(veg.matrix1)!=1) #F
- index4 <- match(dimnames(veg.matrix1)[[1]],species.list)
- species.trait.value.scaled2 <- species.trait.value.scaled[index4]
- species.list2 <- species.list[index4]
- ktab1 <- ktab.list.df(list(as.data.frame(species.trait.value.scaled2)))
- d <- dist.ktab(ktab1, type=c("Q"), option = "noscale")
- #labels(d)
- #str(d)
- #attr(d,"Labels")
- attr(d,"Labels") <- species.list2
- #labels(d)
- #any(is.na(d)) # F
- #dist <- as.matrix(d)
- #as.matrix(d)[1:8,1:8]
- d <- sqrt(d)
- #dist <- as.matrix(d)
- #rm(d)
- if (k<number.packages){
- FD.all[seq((k-1)*package.size+1,k*package.size,1)] <-
- divc(as.data.frame(veg.matrix1),d)$diversity
- #DT2[,list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]$Rao.Q
- any(colSums(veg.matrix1)==0) # F o.K.
- any(colSums(veg.matrix1)!=1) # T ups....
- hist(colSums(veg.matrix1))
- any(rowSums(veg.matrix1)==0) # F o.K.
- } else {
- FD.all[seq((k-1)*package.size+1,number.plots.total,1)] <-
- divc(as.data.frame(veg.matrix1),d)$diversity
- #DT2[,list(Rao.Q = Rao(Relative.cover,Matched.concept)),by=PlotObservationID]
- }
-}
-close(pb)
-#FD.all
-write.csv(FD.all, file = paste(path.sPlot, "FDall.csv", sep = ""), row.names = FALSE)
-