diff --git a/.gitignore b/.gitignore
index 2ea4c35..5b6a065 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,19 +2,3 @@
 .Rhistory
 .RData
 .Ruserdata
-.Rproj
-.txt
-.xls
-.gitignore
-envdianthus.Rproj
-.~lock.rotation.tsv#
-environment.RData
-env.RData
-PCAphylo.RDS
-.RDS
-ouchtable.RDS
-ouchtree.RDS
-resultadosoverlapingcormethod.txt
-resultadosoverlapinggbif.txt
-maxent.jar
-MoranIPCAenv.jpeg
diff --git a/Between-species_analyses.R b/Between-species_analyses.R
deleted file mode 100644
index 4b92f28..0000000
--- a/Between-species_analyses.R
+++ /dev/null
@@ -1,216 +0,0 @@
-### Testing for climatic niche conservatism ###
-
-library(ade4)
-
-# pca1 <-dudi.pca(X)
-# Between-Class Analysis# En nuestro caso corresponde al BETWEEN-occ de Broennimann 2012
-# Nos da el % de varianza (inertia) explicado por los grupos
-Bocc <- bca(pca1, ploidyocc, scannf = FALSE)
-
-#Monte-Carlo Test on the between-groups inertia percentage
-rand1 <- rtest(bca, 99)
-rand1
-plot(rand1, main = "Monte-Carlo test")
-
-###  Within-Class Analysis# En nuestro caso corresponde al WITHIN-occ de Broennimann 2012
-witocc <- wca(pca1, ploidyocc, scan = FALSE, nf = 2)
-
-
-### WITHIN-env. En este caso como grupo habría que poner los backgrounds de las distintas especies
-
-witenc <- wca(pca1, backgroundploidy, scan = FALSE, nf = 2)
-
-
-
-#### Phylogenetic PCA (e.g., Revell 2009; Evolution)
-
-library(phytools)
-phyl.pca(tree, Y, method="BM", mode="cov")
-
-### Phylogenetical signal Enviromental variables Abouheif’s C statistic tested the null hypothesis
-# that traits did not experience phylogenetic autocorrelation (based on the topology)
-
-library(phytools)
-tree<-read.nexus("tree.nex")
-plot(tree)
-orden<-tree$tip.label # este es el orden de las poblaciones para los valores del PCA
-#Extraer los valorehansen(data = datos["PC2"], tree = ouchtree, regimes = datos["OU.4"],     sqrt.alpha = 1, sigma = 1)s del PCA -> PCA1values, PCA2values, PCA3values?
-PCAphylo<-readRDS("PCAphylo.RDS")
-rownames(PCAphylo)<-orden
-
-
-phylosig(tree,PCAphylo[,1],method="lambda",test=TRUE, nsim=100000)
-phylosig(tree,PCAphylo[,1],method="K",test=TRUE, nsim=100000)
-phylosig(tree,PCAphylo[,2],method="lambda",test=TRUE, nsim=100000)
-phylosig(tree,PCAphylo[,2],method="K",test=TRUE, nsim=100000)
-
-library(phangorn)
-tree2<-chronopl(tree, lambda = 1)
-phenogram(tree2,PCAphylo$Axis1)
-xa<-c(x,fastAnc(tree,PCAphylo$Axis1))
-H<-nodeHeights(tree)
-Y<-matrix(xa[tree$edge],nrow(tree$edge),2)
-plot.new()
-plot.window(xlim=range(H),ylim=range(Y))
-axis(1)
-axis(2)
-for(i in 1:nrow(tree$edge)) lines(H[i,],Y[i,],lwd=2)
-title(xlab="time",ylab="phenotype")
-
-
-
-#
-library(adephylo)
-library(phylobase)
-library(phylosignal)
-PCA<-phylo4d(tree,PCAphylo)
-
-correlo1<-phyloCorrelogram(PCA, trait = "Axis1")
-
-correlo2<-phyloCorrelogram(PCA, trait = "Axis2")
-par(mfrow=c(1,2))
-plot(correlo1, main= "PCA-env Axis 1", ylab="Moran's I")
-plot(correlo2, main= "PCA-env Axis 2", ylab="Moran's I")
-opographic wetness index
-table.phylo4d(PCA, ratio.tree = 1/3, cex.label=0.6, legend = F, cex.symbol=0.7,
-              treetype="cladogram", box =F, scale=T, show.var.label=T, symbol="circles", grid=F)
-x <- tdata(PCA, type="tip")
-
-W <- proxTips(tree, met="Abouheif")
-abouTests <- abouheif.moran(PCA)
-plot(abouTests)
-myProx <- vcv.phylo(tree)
-abouTestsBrown <- abouheif.moran(PCA, W=myProx)
-
-a1.ortgTest <- orthogram(x$Axis1, tree) # No señal en eje 1
-a2.ortgTest <- orthogram(x$Axis2, tree) # Sí en eje 2
-
-#Correlogram dist env dist phylo
-
-
-# OUCH
-ouchtree<-readRDS("ouchtree.RDS")
-datosouch<-readRDS("ouchtable.RDS")
-PCAphylo2 <- PCAphylo[as.character(datosouch$labels[-(1:16)]) , ] # Reorder to fix ouchtree
-rownames(PCAphylo2)==datosouch$labels[-(1:16)] # Check git reset --mergeorder
-
-datosouch$PC1[-(1:16)]<-PCAphylo2$Axis1
-datosouch$PC2[-(1:16)]<-PCAphylo2$Axis2
-
-library(ouch)
-#Brownian model
-brownPC1 <- brown(datosouch['PC1'],ouchtree)
-bootbrwonpc1<-bootstrap(brownPC1)
-
-brownPC2 <- brown(datosouch['PC2'],ouchtree)
-bootbrwonpc2<-bootstrap(brownPC2)
-
-#OU model 1 optima
-OU1PC1<-hansen(data = datosouch["PC1"], tree = ouchtree, regimes = datosouch["OU.1"],     sqrt.alpha = 1, sigma = 1)
-OU1PC1boot<-bootstrap(OU1PC1)
-
-OU1PC2<-hansen(data = datosouch["PC2"], tree = ouchtree, regimes = datosouch["OU.1"],     sqrt.alpha = 1, sigma = 1)
-OU1PC2boot<-bootstrap(OU1PC2)
-
-#OU model 4 optima == ploidy levels
-OU4PC1<-hansen(data = datosouch["PC1"], tree = ouchtree, regimes = datosouch["OU.4"],     sqrt.alpha = 1, sigma = 1)
-OU4PC1boot<-bootstrap(OU4PC1)
-
-OU4PC2<-hansen(data = datosouch["PC2"], tree = ouchtree, regimes = datosouch["OU.4"],     sqrt.alpha = 1, sigma = 1)
-OU4PC2boot<-bootstrap(OU4PC2)
-
-#OU model 5 optima == ploidy levels + northern 4x
-OU5PC1<-hansen(data = datosouch["PC1"], tree = ouchtree, regimes = datosouch["OU.5"],     sqrt.alpha = 1, sigma = 1)
-OU5PC1boot<-bootstrap(OU4PC1)
-
-OU5PC2<-hansen(data = datosouch["PC2"], tree = ouchtree, regimes = datosouch["OU.5"],     sqrt.alpha = 1, sigma = 1)
-OU5PC2boot<-bootstrap(OU4PC2)
-
-t.test(OU1PC1boot$aic.c,bootbrwonpc1$aic.c)
-t.test(OU4PC1boot$aic.c,bootbrwonpc1$aic.c)
-t.test(OU1PC1boot$aic.c,OU4PC1boot$aic.c)
-t.test(OU1PC1boot$aic.c,OU5PC1boot$aic.c)
-t.test(OU4PC1boot$aic.c,OU5PC1boot$aic.c)
-t.test(bootbrwonpc1$aic.c,OU5PC1boot$aic.c)
-
-# los test muestrasn que PC1 sigue model OU con 5 optimos
-
-t.test(OU1PC2boot$aic.c,bootbrwonpc2$aic.c)
-t.test(OU4PC2boot$aic.c,bootbrwonpc2$aic.c)
-t.test(OU1PC2boot$aic.c,OU4PC2boot$aic.c)
-t.test(OU1PC2boot$aic.c,OU5PC2boot$aic.c)
-t.test(OU4PC2boot$aic.c,OU5PC2boot$aic.c)
-t.test(bootbrwonpc2$aic.c,OU5PC2boot$aic.c)
-# los test muestran que PC2 sigue OU con 4 optima
-
-# PC1
-quantile(bootbrwonpc1$aic.c,c(0.05,0.5,0.95))
-quantile(OU1PC1boot$aic.c,c(0.05,0.5,0.95))
-quantile(OU4PC1boot$aic.c,c(0.05,0.5,0.95))
-quantile(OU5PC1boot$aic.c,c(0.05,0.5,0.95))
-
-#PC2
-quantile(bootbrwonpc2$aic.c,c(0.05,0.5,0.95))
-quantile(OU1PC2boot$aic.c,c(0.05,0.5,0.95))
-quantile(OU4PC2boot$aic.c,c(0.05,0.5,0.95))
-quantile(OU5PC2boot$aic.c,c(0.05,0.5,0.95))
-
-boxplot(datosouch$PC1~factor(datosouch$OU.4,levels=c("2x","4x","6x","12x")), col=rainbow(4))
-boxplot(datosouch$PC2~factor(datosouch$OU.4,levels=c("2x","4x","6x","12x")),col=rainbow(4))
-
-
-# El mismo analysis BM, OU1, OUM
-library(mvMORPH)
-state<-datosouch$OU.4
-names(state)<-datosouch$labels
-# Make the tree with mapped states using SIMMAP
-tree<-  make.simmap(tree, state, model="ER", nsim=1)
-plot(tree)
-
-trait1_BM<-  mvBM(tree, PCAphylo[,1], model="BMM")
-trait2_BM<-  mvBM(tree, PCAphylo[,2], model="BMM")
-
-trait1_OU1<-  mvOU(tree, PCAphylo[,1], model="OU1")
-trait2_OU1<-  mvOU(tree, PCAphylo[,2], model="OU1")
-
-trait1_OUM<-  mvOU(tree, PCAphylo[,1], model="OUM")
-trait2_OUM<-  mvOU(tree, PCAphylo[,2], model="OUM")
-
-AIC(trait1_BM);AIC(trait1_OUM);AIC(trait1_OU1)
-AIC(trait1_BM);AIC(trait2_OUM);AIC(trait2_OU1)
-
-OUM<-  mvOU(tree, PCAphylo, model="OUM")
-
-# Ajusta modelo BM con un único régimen o varias tasas de evolución
-library(phytools)
-fitBM_PCA1<-brownie.lite(tree,PCAphylo[,1], test="simulation")
-fitBM_PCA2<-brownie.lite(tree,PCAphylo[,2], test="simulation")
-
-# Ajusta OU
-library(OUwie)
-plotSimmap(tree,type="fan",fsize=0.8,ftype="i")
-tree$node.label<-getStates(tree,"nodes")
-dataPC1<-data.frame(Genus_species=PCA@label,Reg=factor(c(2,2,2,2,6,6,6,4,4,4,4,4,4,12,12,12,12,4,4,4,2,2)),X=PCA@data$Axis1)
-
-fitBMPC1<-OUwie(tree,dataPC1,model="BM1",simmap.tree=TRUE)
-fitOUMPC1<-OUwie(tree,dataPC1,model="OUM",simmap.tree=TRUE)
-fitOUMAPC1<-OUwie(tree,dataPC1,model="OUMA",simmap.tree=TRUE)
-fitOUMVAPC1<-OUwie(tree,dataPC1,model="OUMVA",simmap.tree=TRUE)
-
-fitBMPC1$AICc
-fitOUMPC1$AICc
-fitOUMAPC1$AICc
-fitOUMVAPC1$AICc
-
-
-dataPC2<-data.frame(Genus_species=PCA@label,Reg=factor(c(2,2,2,2,6,6,6,4,4,4,4,4,4,12,12,12,12,4,4,4,2,2)),X=PCA@data$Axis2)
-
-fitBMPC2<-OUwie(tree,dataPC2,model="BM1",simmap.tree=TRUE)
-fitOUMPC2<-OUwie(tree,dataPC2,model="OUM",simmap.tree=TRUE)
-fitOUMAPC2<-OUwie(tree,dataPC2,model="OUMA",simmap.tree=TRUE)
-fitOUMVAPC2<-OUwie(tree,dataPC2,model="OUMVA",simmap.tree=TRUE)
-
-fitBMPC2$AICc
-fitOUMPC2$AICc
-fitOUMAPC2$AICc
-fitOUMVAPC2$AICc
diff --git a/CSD.R b/CSD.R
deleted file mode 100644
index 086d80d..0000000
--- a/CSD.R
+++ /dev/null
@@ -1,21 +0,0 @@
-#Continuous species distribution
-cooDbroterigbifdata<-readRDS("cooDbroterigbifdata.RDS")
-
-rast <- raster(ncol = 100, nrow = 100)
-cooDbroterigbifdata$ploidy
-extent(rast) <- extent(cooDbroterigbifdata)
-x<-rasterize(cooDbroterigbifdata, rast, as.numeric(cooDbroterigbifdata$ploidy), fun = mean)
-plot(x)
-
-
-library(spatstat)
-crn.pp = as(cooDbroterigbifdata, 'ppp')
-
-q = quadratcount(crn.pp, 10, 10)
-q2 = quadratcount(crn.pp[crn.pp$marks=="2x"], 10, 10)
-q4 = quadratcount(crn.pp[crn.pp$marks=="4x"], 10, 10)
-q6 = quadratcount(crn.pp[crn.pp$marks=="6x"], 10, 10)
-q12 = quadratcount(crn.pp[crn.pp$marks=="12x"], 10, 10)
-plot(q2)
-plot(crn.pp, add=T)
-coa<-dudi.coa(as.data.frame(matrix(q, ncol=10)))
diff --git a/NicheComparisonTests.R b/NicheComparisonTests.R
new file mode 100644
index 0000000..cc9f7e6
--- /dev/null
+++ b/NicheComparisonTests.R
@@ -0,0 +1,425 @@
+library (dismo)
+library ('rgbif')
+library (rgdal)
+library (raster)
+library (gtools)
+library (rgeos)
+library (rjson)
+library (sp)
+library (GSIF)
+library (plotrix)
+library (ggmap)
+library (spatstat)
+library (ecospat)
+library (ggbiplot)
+library (spatialEco)
+library (fmsb)
+library (Hmisc)
+library (devtools)
+library (ENMTools)
+library (vioplot)
+library (adegraphics)
+library (lattice)
+
+
+###Downloading GBIF data
+Dbroterigbif <- occ_search (scientificName = "Dianthus broteri",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
+Dinoxgbif <- occ_search (scientificName = "Dianthus inoxianus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
+Dhinoxgbif <- occ_search (scientificName = "Dianthus hinoxianus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
+Dvalengbif <- occ_search (scientificName = "Dianthus valentinus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
+
+
+###Cleaning of location data
+Dbroterigbif <- Dbroterigbif [Dbroterigbif$year>2004,]
+Dbroterigbif <- Dbroterigbif [-c(129:189),]
+Dhinoxgbif <- Dhinoxgbif [Dhinoxgbif$year>2004,]
+Dhinoxgbif <- Dhinoxgbif [1,]
+Dvalengbif <- Dvalengbif [1,]
+
+Dbroterigbif <- Dbroterigbif [, c(1:3)]
+Dinoxgbif <- Dinoxgbif [, c(1:3)]
+Dhinoxgbif <- Dhinoxgbif [, c(1:3)]
+Dvalengbif <- Dvalengbif [, c(1:3)]
+
+todo <- rbind (Dbroterigbif,Dinoxgbif,Dhinoxgbif,Dvalengbif)
+todo$name <- "Dianthus broteri"
+gbifmap (todo, mapdatabase = "world", region = "Spain")
+
+
+###Cleaning of duplicate records
+todo_cleaned <- unique(todo)
+todo_cleaned <- todo_cleaned [,-1]
+
+
+###Map
+e <- extent (-10,3.5,35.5,44)
+coordinates(todo_cleaned) <- ~decimalLongitude+ decimalLatitude
+crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
+proj4string (todo_cleaned) <- crs.geo
+map <- plot (gmap (e, type = "satellite"))
+points <- points (Mercator(todo_cleaned), col = "red", pch=20, cex = 1.5)
+
+
+###Thinning the records to one per 1 km2 
+r <- raster(todo_cleaned)
+res(r) <- 0.008333333
+r <- extend(r, extent(r)+1)
+cooDbroterigbif <- as.data.frame(gridSample(todo_cleaned, r, n=1))
+
+
+###Additional populations from our records
+e1 <- extent (-7,-6.4,36.8,37.5)
+
+dianthuspops <- read.table ("FILE.txt", header = T, sep = ",")
+dianthuspops <- dianthuspops [,-c(1,4)]
+dianthuspops <- dianthuspops [,c(2,1)]
+colnames (dianthuspops) <- c("long", "lat")
+
+coordinates(dianthuspops) <- ~long+ lat
+proj4string(dianthuspops) <- crs.geo
+r1 <- raster(dianthuspops)
+res(r1) <- 0.008333333
+r1 <- extend(r1, extent(r1)+1)
+dianthuspops_fran <- as.data.frame(gridSample(dianthuspops, r1, n=1))
+
+coordinates(dianthuspops_fran) <- ~long+ lat
+proj4string(dianthuspops_fran) <- crs.geo
+plot(gmap(e1, type = "satellite"))
+points(Mercator(dianthuspops_fran), col="red" , pch=20, cex=1)
+
+dianthuspops_fran <- as.data.frame (dianthuspops_fran)
+colnames(cooDbroterigbif) <- colnames(dianthuspops_fran)
+cooDbroterigbif.def <- rbind (cooDbroterigbif, dianthuspops_fran)
+cooDbroterigbif.def2 <- cooDbroterigbif.def
+
+coordinates(cooDbroterigbif.def2) <- ~long+ lat
+proj4string (cooDbroterigbif.def2) <- crs.geo
+
+r <- raster(cooDbroterigbif.def2)
+res(r) <- 0.008333333
+r <- extend(r, extent(r)+1)
+cooDbroterigbif.def.data <- as.data.frame(gridSample(cooDbroterigbif.def2, r, n=1))
+
+cooDbroterigbif.def.data <- cbind (cooDbroterigbif.def.data, 1)
+colnames(cooDbroterigbif.def.data)[3]<-"ploidy"
+
+###Ploidy estimation based on geographic coordinates
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38 & cooDbroterigbif.def.data$long > -9.5 & cooDbroterigbif.def.data$lat < 39 & cooDbroterigbif.def.data$long < -8] <- "4x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 37 & cooDbroterigbif.def.data$long > -9 & cooDbroterigbif.def.data$lat < 37.5 & cooDbroterigbif.def.data$long < -7.5] <- "2x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36.9 & cooDbroterigbif.def.data$long > -7.4 & cooDbroterigbif.def.data$lat < 37.6 & cooDbroterigbif.def.data$long < -6] <- "12x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36 & cooDbroterigbif.def.data$long > -6.3 & cooDbroterigbif.def.data$lat < 37.14 & cooDbroterigbif.def.data$long < -4.72] <- "4x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36.5 & cooDbroterigbif.def.data$long > -4.7 & cooDbroterigbif.def.data$lat < 38.1 & cooDbroterigbif.def.data$long < -1.67] <- "2x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 37.3 & cooDbroterigbif.def.data$long > -1.7 & cooDbroterigbif.def.data$lat < 38.3 & cooDbroterigbif.def.data$long < -0.5] <- "6x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38.2 & cooDbroterigbif.def.data$long > -2.33 & cooDbroterigbif.def.data$lat < 38.83 & cooDbroterigbif.def.data$long < 0.33] <- "6x"
+cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38.83 & cooDbroterigbif.def.data$long > -1.7 & cooDbroterigbif.def.data$lat < 40.9 & cooDbroterigbif.def.data$long < 0.38] <- "4x"
+cooDbroterigbif.def.data$ploidy[c(58,62:63,120:123,125:130)]<- "4x"
+cooDbroterigbif.def.data <- cooDbroterigbif.def.data[cooDbroterigbif.def.data$ploidy!="1",]
+
+
+ploidy<-factor(cooDbroterigbif.def.data$ploidy,levels = c("2x","4x","6x","12x"),ordered = TRUE)
+
+coordinates(cooDbroterigbif.def.data) <- ~long+ lat
+proj4string(cooDbroterigbif.def.data) <- crs.geo
+plot(gmap(e, type = "satellite"))
+points(Mercator(cooDbroterigbif.def.data), col=ploidy, pch=20, cex=1)
+
+plot(gmap(e1, type = "satellite"))
+points(Mercator(cooDbroterigbif.def.data), col=ploidy, pch=20, cex=1)
+
+
+cooDbroterigbif.def.vars <- as.data.frame (cooDbroterigbif.def.data)
+cooDbroterigbif.def.vars <- cooDbroterigbif.def.vars [,c(3,1,2)]
+coordinates(cooDbroterigbif.def.vars) <- ~long+ lat
+proj4string(cooDbroterigbif.def.vars) <- crs.geo
+
+
+###Downloading of predictor variables (chelsa, envirem, elevation, SoilGrids)
+###Data extraction and merging into the same extent
+
+#cooDbroterigbif.def.vars <- read.csv2 ("FILE.csv")
+ploidy<-factor(cooDbroterigbif.def.vars$ploidy,levels = c("2x","4x","6x","12x"),ordered = TRUE)
+
+crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
+coordinates(cooDbroterigbif.def.vars) <- ~long+ lat
+proj4string(cooDbroterigbif.def.vars) <- crs.geo
+
+e <- extent (-10,4.5,35.5,44)
+
+chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/chelsa", pattern = ".tif", full.names = TRUE))
+chelsa <- stack (chelsafiles)
+che.c <- crop (chelsa,e)
+
+enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/envirem", pattern = ".bil", full.names = TRUE))
+envirem <- stack (enviremfiles)
+env.c <- crop (envirem, e)
+
+alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
+alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
+alt15 <- stack (alt15files)
+alt16 <- stack (alt16files)
+alt.m <- merge (alt15, alt16, ext=e)
+
+variables <- stack (che.c, env.c, alt.m)
+names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","climaticMoistureIndex","continentality","growingDegDays5","maxTempColdest","minTempWarmest","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","topoWet", "elevation")
+
+soilgrids<-extract.list(cooDbroterigbif.def.vars, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas", ID = "ploidy")
+colnames (soilgrids) <- c("ploidy","AWCh2","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT")
+soilgrids$ploidy<-factor(soilgrids$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="1","clay")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="2","silty clay")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="3","sandy clay")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="4","clay loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="5","silty clay loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="6","sandy clay loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="7","loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="8","silty loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="9","sandy loam")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="10","silt")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="11","loamy sand")
+# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="12","sand")
+# soilgrids$TEXMHT<-factor(soilgrids$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
+# soilgrids<-cbind(soilgrids,apply(soilgrids[,c(2:4)], 1, mean))
+# soilgrids<-soilgrids[,-c(1:4)]
+colnames(soilgrids)[2]<-"AWC"
+
+
+presvals <- extract (variables, cooDbroterigbif.def.vars)
+presvals <- cbind (ploidy, presvals, soilgrids) 
+presvals$PHIHOX <- presvals$PHIHOX/10
+presvals <- presvals [,-34]
+
+###tri calculation from SpatilEco package
+tri.ext <- tri(alt.m)
+projection(tri.ext) <- crs.geo 
+trivalues<-extract(tri.ext,cooDbroterigbif.def.vars)
+presvals <- cbind (presvals, trivalues)
+presvals <- presvals[,c(1:31,40,32:39)]
+colnames(presvals)[32] <- "tri"
+presvals2 <- na.omit (presvals)
+ploidy <- as.data.frame(presvals2[,1])
+colnames (ploidy) <- "ploidy"
+
+
+###Removing correlated variables
+###PCA only with presence points using non-correlated variables
+presvals.pca <- presvals2[,-1]
+selected <- vif_func(presvals.pca) 
+presvals.pca.2 <- presvals.pca[,c(selected)]
+
+ploidy <- ploidy$ploidy
+ploidy <- factor (ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
+
+pca <- prcomp(presvals.pca.2, scale. = TRUE, retx = T)
+ggbiplot(pca, obs.scale = 1,var.scale = 1,
+         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
+  scale_color_discrete(name = '') +
+  geom_point(aes(colour=ploidy), size = 3) +
+  theme(legend.direction = 'vertical', legend.position = 'right')
+
+
+###Creating the general background
+dbroteridata <- as.data.frame(cooDbroterigbif.def.vars)
+dbroteridata <- dbroteridata [,-1]
+presvalsdata <- cbind (dbroteridata, presvals)
+presvalsdata <- na.omit (presvalsdata)
+presvals2 <- presvalsdata
+coordinates(presvals2) <- ~long+ lat
+proj4string(presvals2) <- crs.geo
+
+backgroundcoord <- presvalsdata [,c(1,2)]
+coordinates(backgroundcoord) <- ~long+ lat
+proj4string(backgroundcoord) <- crs.geo
+
+mask <- raster(backgroundcoord)
+res(mask) <- 0.008333333
+x <- circles(backgroundcoord, d=150000, lonlat=TRUE)
+pol <- gUnaryUnion(x@polygons)
+samp <- spsample(pol, 10000, type='random', iter=25)
+extent(mask) <- extent(pol)
+cells <- cellFromXY(mask, samp)
+length(cells)
+cells <- unique(cells)
+length(cells)
+xy <- xyFromCell(mask, cells)
+xy <- as.data.frame(xy)
+colnames(xy) <- c("long", "lat")
+
+###Thinning the records to one per 1 km2 
+coordinates(xy) <- ~long+ lat
+proj4string(xy) <- crs.geo
+r <- raster(xy)
+res(r) <- 0.008333333
+r <- extend(r, extent(r)+1)
+backcoord_sel <- as.data.frame(gridSample(xy, r, n=1))
+coordinates(backcoord_sel) <- ~long+ lat
+proj4string(backcoord_sel) <- crs.geo
+
+###Data extraction of background points and modification of the dataset
+backgroundclim <- extract(variables,backcoord_sel)
+backgroundsoil <- extract.list(backcoord_sel, list.files("FILE", ID = "ploidy"))
+backgrounddat <- cbind("background",as.data.frame(backcoord_sel),backgroundclim, backgroundsoil)
+backgrounddat <- backgrounddat[,-36]
+
+coordinates(backgrounddat) <- ~long+ lat
+proj4string(backgrounddat) <- crs.geo
+trivalues1 <- extract(tri.ext,backgrounddat)
+backgrounddat <- as.data.frame(backgrounddat)
+backgrounddat <- cbind (backgrounddat, trivalues1)
+backgrounddat <- backgrounddat[,c(1:33,42,34:41)]
+colnames(backgrounddat)[34] <- "tri"
+
+backgrounddat.c <- na.omit(backgrounddat)
+backgrounddat.c <- backgrounddat.c[,c(2,3,1,4:42)]
+colnames(backgrounddat.c) <- colnames(presvalsdata)
+
+#Plotting a map with presence and background points
+coordinates(backgrounddat.c) <- ~long+ lat
+proj4string(backgrounddat.c) <- crs.geo
+plot(gmap(e, type = "satellite"))
+points(Mercator(backgrounddat.c), col = 'orange', pch=20, cex=1)
+points(Mercator(presvals2), col=presvals2$ploidy, pch=20, cex=1)
+
+
+###PCA with presence and background points
+tododef <- rbind (presvalsdata, as.data.frame(backgrounddat.c))
+
+#tododef <- read.csv2 ("tododef.csv")
+
+todo.pca <- tododef[,-c(1:4)]
+selected2<-vif_func(todo.pca)
+todo.pca.2 <- todo.pca[,c(selected2)]
+
+todoploidy <- factor (tododef$ploidy, levels = c("2x", "4xs", "4xe", "6x", "12x", "background"), ordered = TRUE)
+w<-c(rep(0,nrow(presvalsdata)),rep(1,5155))
+todoploidy2<-todoploidy[todoploidy!="background"]
+todoploidy2 <- factor (todoploidy2, levels = c("2x", "4xs", "4xe", "6x", "12x"), ordered = T)
+
+library (adegraphics)
+frame.plot=FALSE
+
+pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
+gcol = c("green", "darkred", "red", "blue", "purple")
+adegpar(pellipses.alpha=0.3,ppoints.cex=1.5,plabels =list(cex=1.5, optim=T, alpha=0.8,boxes=list(draw=T, lwd=1.5, alpha=0.8)), plines.lwd = 2)
+g1<-s.corcircle(pcaback$co/80, label = row.names(pcaback$co), fullcircle = TRUE)
+g1_1<-update (g1, xlim = c(-2.4,2.2), ylim = c(-1,1.01))
+adegpar(pellipses.alpha=0.3,ppoints.cex=1.5,plabels =list(cex=1.4, optim=T, alpha=1,boxes=list(draw=T, lwd=2, alpha=1)), plines.lwd = 2)
+g2<-s.class(pcaback$li[1:143,], todoploidy2, col = gcol, porigin.origin=c(-0.5,-5), porigin.draw=T)
+g2_2<-update (g2, xlim = c(-4.2,4.4), ylim = c(-9.2,-0.8))
+
+
+###Niche comparison tests: cytotypes and 4x lineages (package ECOSPAT)
+row.di<-which(tododef[,3] == "2x")
+row.tes<-which(tododef[,3] == "4xs")
+row.tee<-which(tododef[,3] == "4xe")
+row.te<-which(tododef[,4] == "4x")
+row.he<-which(tododef[,3] == "6x")
+row.do<-which(tododef[,3] == "12x")
+
+
+scores.clim<- pcaback$li 
+scores.di<- pcaback$li[row.di,]
+scores.te<- pcaback$li[row.te,]
+scores.tes<- pcaback$li[row.tes,]	
+scores.tee<- pcaback$li[row.tee,]	
+scores.he<- pcaback$li[row.he,]	
+scores.do<- pcaback$li[row.do,]	
+
+zdi<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.di, R=500)
+zte<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.te, R=500)
+ztes<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.tes, R=500)
+ztee<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.tee, R=500)
+zhe<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.he, R=500)
+zdo<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.do, R=500)
+
+ecospat.plot.niche (zdi)
+ecospat.plot.niche (zte)
+ecospat.plot.niche (ztes)
+ecospat.plot.niche (ztee)
+ecospat.plot.niche (zhe)
+ecospat.plot.niche (zdo)
+
+###Equivalency tests
+equivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 100, alternative = "lower")
+equivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 100, alternative = "lower")
+equivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 100, alternative = "lower")
+equivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 100, alternative = "lower")
+equivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 100, alternative = "lower")
+equivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 100, alternative = "lower")
+
+equivalency.test.dites<-ecospat.niche.equivalency.test (zdi, ztes, 100, alternative = "lower")
+equivalency.test.ditee<-ecospat.niche.equivalency.test (zdi, ztee, 100, alternative = "lower")
+equivalency.test.hetes<-ecospat.niche.equivalency.test (zhe, ztes, 100, alternative = "lower")
+equivalency.test.hetee<-ecospat.niche.equivalency.test (zhe, ztee, 100, alternative = "lower")
+equivalency.test.dotes<-ecospat.niche.equivalency.test (zdo, ztes, 100, alternative = "lower")
+equivalency.test.dotee<-ecospat.niche.equivalency.test (zdo, ztee, 100, alternative = "lower")
+
+###Overlap tests
+overlap.test.dite<-ecospat.niche.overlap (zdi, zte, cor=FALSE)
+overlap.test.dihe<-ecospat.niche.overlap (zdi, zhe, cor=FALSE)
+overlap.test.dido<-ecospat.niche.overlap (zdi, zdo, cor=FALSE)
+overlap.test.tehe<-ecospat.niche.overlap (zte, zhe, cor=FALSE)
+overlap.test.tedo<-ecospat.niche.overlap (zte, zdo, cor=FALSE)
+overlap.test.hedo<-ecospat.niche.overlap (zhe, zdo, cor=FALSE)
+
+###Similarity tests
+similarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 100, alternative = "greater")
+similarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 100, alternative = "greater")
+similarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 100, alternative = "greater")
+similarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 100, alternative = "greater")
+similarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 100, alternative = "greater")
+similarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 100, alternative = "greater")
+similarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 100, alternative = "greater")
+similarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 100, alternative = "greater")
+similarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 100, alternative = "greater")
+similarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 100, alternative = "greater")
+similarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 100, alternative = "greater")
+similarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 100, alternative = "greater")
+
+
+###Vioplots
+par(cex.lab=1.5, cex.axis=1.2, cex.main = 2)
+vio1<-vioplot_fun (scores.di$Axis1, scores.te$Axis1, scores.tes$Axis1, scores.tee$Axis1, scores.he$Axis1, scores.do$Axis1, names = c("2x","4x","4xs","4xe","6x","12x"), col = c("green", "orange", "darkred", "red", "blue", "purple"))
+title(xlab="Ploidy level", main = "Axis 1")
+vioplot_fun (scores.di$Axis2, scores.te$Axis2, scores.tes$Axis2, scores.tee$Axis2, scores.he$Axis2, scores.do$Axis2, names = c("2x","4x","4xs","4xe","6x","12x"), col = c("green", "orange", "darkred", "red", "blue", "purple"))
+title(xlab="Ploidy level", main = "Axis 2")
+
+
+###Niche breadth
+bredi<-raster.breadth (zdi$w)
+brete<-raster.breadth (zte$w)
+brehe<-raster.breadth (zhe$w)
+bredo<-raster.breadth (zdo$w)
+bretes<-raster.breadth (ztes$w)
+bretee<-raster.breadth (ztee$w)
+
+
+###Niche divergence test
+#PC1
+PC1dn<-t.test(scores.di[,1], scores.tes[,1]) 
+PC1dn1<-t.test(scores.di[,1], scores.tee[,1]) 
+PC1dn2<-t.test(scores.tes[,1], scores.tee[,1]) 
+PC1dn3<-t.test(scores.he[,1], scores.tes[,1]) 
+PC1dn4<-t.test(scores.he[,1], scores.tee[,1]) 
+PC1dn5<-t.test(scores.do[,1], scores.tes[,1]) 
+PC1dn6<-t.test(scores.do[,1], scores.tee[,1])
+PC1dn7<-t.test(scores.di[,1], scores.te[,1]) 
+PC1dn8<-t.test(scores.di[,1], scores.he[,1]) 
+PC1dn9<-t.test(scores.di[,1], scores.do[,1]) 
+PC1dn10<-t.test(scores.te[,1], scores.he[,1]) 
+PC1dn11<-t.test(scores.te[,1], scores.do[,1]) 
+PC1dn12<-t.test(scores.he[,1], scores.do[,1]) 
+
+#PC2
+PC2dn<-t.test(scores.di[,2], scores.tes[,2]) 
+PC2dn1<-t.test(scores.di[,2], scores.tee[,2]) 
+PC2dn2<-t.test(scores.tes[,2], scores.tee[,2]) 
+PC2dn3<-t.test(scores.he[,2], scores.tes[,2]) 
+PC2dn4<-t.test(scores.he[,2], scores.tee[,2]) 
+PC2dn5<-t.test(scores.do[,2], scores.tes[,2]) 
+PC2dn6<-t.test(scores.do[,2], scores.tee[,2])
+PC2dn7<-t.test(scores.di[,2], scores.te[,2]) 
+PC2dn8<-t.test(scores.di[,2], scores.he[,2]) 
+PC2dn9<-t.test(scores.di[,2], scores.do[,2]) 
+PC2dn10<-t.test(scores.te[,2], scores.he[,2]) 
+PC2dn11<-t.test(scores.te[,2], scores.do[,2]) 
+PC2dn12<-t.test(scores.he[,2], scores.do[,2])
diff --git a/NicheModels.R b/NicheModels.R
new file mode 100644
index 0000000..0ef8be7
--- /dev/null
+++ b/NicheModels.R
@@ -0,0 +1,59 @@
+library (zoon)
+library (future)
+Sys.setenv(NOAWT=TRUE)
+options(java.parameters="-Xmx120g")
+plan(multicore)
+LoadModule('PrintMap1.R')
+vars.stack.gbif <- stack("stack_zoon_gbif.grd")
+workgbif1 <- function(){ workflow (occurrence = LocalOccurrenceData (filename = "OCURRENCEFILE.csv", occurrenceType = 'presence/absence'),
+                                   covariate  = LocalRaster(vars.stack.gbif),
+                                   process = Chain (StandardiseCov, Crossvalidate (k = 10)),
+                                   model = MaxEnt,
+                                   output = Chain (PrintMap1 (dir = "DIR"), PerformanceMeasures))}
+
+workresgbif1 <- workgbif1()
+save (workresgbif1, file = 'workflow_FILE.RData')
+
+
+###Testing significant deviation from random expectation
+load ("workflow_FILE.RData")
+library (dismo)
+
+back<-workresgbif1[["process.output"]][[1]][["df"]]
+back<-back[,-c(1:5)]
+back<-back[,c(1:12,19,13:18)]
+
+nullModel <- function (x, n = 10, rep = 99)
+{
+  e <- list()
+  #stopifnot(n < nrow(x))
+  for (r in 1:rep) {
+    z <- sample(nrow(x), n)
+    pres <- x[z, ] # randomly drawn presences
+    absc <- x[-z, ] # remaining points set as absences
+    rs10k <- sample(1:nrow(absc), 9500) # random select 10k rows
+    absc <- absc[rs10k,]
+    d <- rbind(pres, absc)
+    v <- c(rep(1, nrow(pres)), rep(0, nrow(absc)))
+    m.null <- maxent(d,v, args = c("noproduct", "nothreshold", "nohinge"))
+    
+    e[[r]] <- evaluate(pres, absc, m.null)
+    cat("-")
+    if (r%%50 == 0) 
+      cat(" ", r, "\n")
+    flush.console()
+  }
+  if (r%%50 != 0) {
+    cat(" ", r, "\n")
+  }
+  else {
+    cat("\n")
+  }
+  e
+}
+
+nulltest <- nullModel (back, n = X, rep = 99)
+nulltest
+
+mean(sapply(nulltest, function(x)x@auc))
+max(sapply(nulltest, function(x)x@auc))
\ No newline at end of file
diff --git a/Nichedianthus.R b/Nichedianthus.R
deleted file mode 100644
index f655509..0000000
--- a/Nichedianthus.R
+++ /dev/null
@@ -1,201 +0,0 @@
-#R script dianthus polyploids model divergence
-# This script downloads species coordinates from gbif. 
-# it clips the pH information from SoilGrids1km
-# it clip the wordclim data
-
-library(dismo)
-library('rgbif')
-library(rgdal)
-library(raster)
-library(gtools)
-library(rjson)
-library(sp)
-library(GSIF)
-library(maps)
-library(rgeos)
-library(ade4)
-library(caret)
-
-#Dbroteri data gbif
-Dbroteriloc<-occ_search(scientificName = "Dianthus broteri",hasGeospatialIssue=FALSE,fields=c('name','decimalLatitude','decimalLongitude', 'country'), limit=10000 ,hasCoordinate=TRUE, basisOfRecord="PRESERVED_SPECIMEN")
-dups <- duplicated(Dbroteriloc[[3]][,c(2,3)]) #Look for duplicates
-brot_cleaned<- Dbroteriloc[[3]][!dups, ] #Remove duplicates
-
-summary(as.factor(fuc_cleaned$country)) #Number ocurrences by country
-
-
-
-# #Dbroteri owndata
-# owndata<-read.table("datos.txt", header=F)
-# colnames(owndata)<- colnames(brot_cleaned)
-
-
-
-
-coordinates(brot_cleaned)<- ~decimalLongitude+ decimalLatitude
-crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string(brot_cleaned) <- crs.geo 
-
-# #Remove ocurrences within 10km2
-# r <- raster(coofuc_final)
-# res(r) <- 0.08333333
-# r <- extend(r, extent(r)+1)
-# coofuc_final_sel <- as.data.frame(gridSample(coofuc_final, r, n=1))
-# coordinates(coofuc_final_sel)<- ~decimalLongitude+ decimalLatitude
-# proj4string(coofuc_final_sel) <- crs.geo 
-# # p <- rasterToPolygons(r)
-# # plot(p, border='gray')
-# # points(coofuc_final)
-# # points(acsel, cex=1, col='red', pch='x')
-
-
-
-#MAP
-map('world', xlim=c(-11,6), ylim=c(35,45))
-box()
-title(main="D. broteri")
-points(brot_cleaned, pch=16, cex=0.4, col="orange")
-
-
-#BIOCLIM DATA
-#Read the raster, stack the raster, crop and merge
-list.ras <- mixedsort(list.files("/home/fbalao/Datos/Dactylorhiza/BDGeoDact/wc0.5/", full.names = T, pattern = ".bil"))
-bio <- stack(list.ras)
-projection(bio)<-crs.geo
-
-#Extract the bioclim variables
-envbrot<-extract(bio,brot_cleaned)
-
-  
-#Tree coverage
-vcf<-raster("/home/fbalao/MODIS/out.tif")
-vcfbrot<-extract(vcf,brot_cleaned)
-vcfbrot[vcfbrot==200]<-0
-
-envbrot<-cbind(envbrot,vcfbrot)
-envbrot<-na.omit(envbrot)
-
-#PCA
-library(ade4)
-pca1 <- dudi.pca(envbrot, scann = FALSE, nf = 3)
-scatter(pca1, clab.row = 0, posieig = "none", add.plot=T, clab.col = 0.8, col=2)
-gcol<-c("blue", "red", "green", "yellow")
-s.class(pca1$li, as.factor(data[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-7,-1,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1.8)
-
-# MANOVA on PC Axis y plidy levels
-
-summary(manova(as.matrix(pca1$li)~as.factor(ploidy)))
-
-boxplot(pca1$li$Axis1~as.factor(ploidy), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(ploidy), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(ploidy), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-#Figura completa
-
-biolegend<-read.table("/home/fbalao/Datos/Proyectos/ProyectoExcelencia2016/Figuras/PCA_NichoEcologico_Dianthusbroteri/biolab.txt", sep="@")
-layout(matrix(c(1,2,2,3,4,5),nrow=2, byrow =T ))
-plot.new()
-addtable2plot(-0.7,0,biolegend,display.rownames=F,hlines=F, cex=0.5, display.colnames = F, bty="n")
-scatter(pca1, clab.row = 0, posieig = "none",clab.col = 1, col=2, xlim=c(10,10))
-s.class(pca1$li, as.factor(data[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-6.8,-0.1,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1, cex=0.7)
-boxplot(pca1$li$Axis1~as.factor(ploidy), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(ploidy), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(ploidy), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-# Kernel plots
-# plot(density(data$bio5[data$ploidy==2]), col="blue", xlim=c(200,500), ylim=c(0,0.25))
-# lines(density(data$bio5[data$ploidy==4]), col="red")
-# lines(density(data$bio5[data$ploidy==6]), col="green")
-# lines(density(data$bio5[data$ploidy==12]), col="yellow")
-
-#Soil
-soilgrids.r <- REST.SoilGrids(c("ORCDRC","PHIHOX"))
-ovfuc <- over(soilgrids.r, coofuc)
-ovinc <- over(soilgrids.r, cooinc)
-#pH
-library(lattice)
-library(aqp)
-data(soil.legends)
-plot(density(na.omit(ovfuc$PHIHOX.M.sd1)/10), col=3)
-lines(density(na.omit(ovinc$PHIHOX.M.sd1)/10), col=2)
-
-PHIHOX.range = range(soil.legends[["PHIHOX"]]$MIN, soil.legends[["PHIHOX"]]$MAX)
-
-
-
-#BIOCLIM DATA
-list.ras <- mixedsort(list.files("/home/fbalao/Datos/Dactylorhiza/BDGeoDact/wc2-5/", full.names = T, pattern = ".bil"))
-bio <- stack(list.ras)
-bio.brick <- brick(bio)
-newext <- drawExtent()
-bio.c<-crop(bio, newext)
-projection(bio.c)<-crs.geo
-
-prevfuc<-extract(bio.brick,coofuc)
-previnc<-extract(bio.brick,cooinc)
-
-biofuc<-cbind(as.data.frame(coofuc),prevfuc)
-bioinc<-cbind(as.data.frame(cooinc),previnc)
-bioinc<-na.omit(bioinc)
-library('caret')
-rem<-findCorrelation(cor(biofuc[,-c(1,2)]), cutoff = .90, verbose = F, names=F)
-rem<-rem+2
-biofuc.c<-biofuc[, -c(1,2,rem)]
-
-pca <- prcomp(biofuc.c)
-plot(pca$x[,1:2], pch=16, main=biofuc.c[1,"Species"]) # 2 first principal components
-z
-
-merged<-rbind(biofuc,bioinc)
-sp<-c(rep('fuc',348), rep('inc',924))
-merged<-cbind(sp,merged)
-mergedtable<-merged[,-c(2,3)]
-mergedtable
-
-#PCA on the climatic variables
-library(ade4)
-pca1 <- dudi.pca(mergedtable[,-c(1,rem+1)], scann = FALSE, nf = 3)
-barplot(pca1$eig)
-gcol = c("blue", "red", "green", "red")
-
-# s.class(dfxy = pca1$li, fac = as.factor(biofuc[,2]), col = gcol, xax = 1, yax = 2)
-# s.corcircle(pca1$co, xax = 1, yax = 2)
-# scatter(pca1, clab.row = 0, posieig = "none")
-# s.class(pca1$li, mergedtable[,1], col = gcol, add.plot = TRUE, cstar = 1.5, clabel = 0, cellipse = 1.5)
-
-
-biolegend<-read.table("/home/fbalao/Datos/Proyectos/ProyectoExcelencia2016/Figuras/PCA_NichoEcologico_Dianthusbroteri/biolab.txt", sep="@")
-layout(matrix(c(1,2,2,3,4,5),nrow=2, byrow =T ))
-plot.new()
-addtable2plot(-0,0,biolegend,display.rownames=F,hlines=F, cex=0.5, display.colnames = F, bty="n")
-scatter(pca1, clab.row = 0, posieig = "none",clab.col = 1, col=2, xlim=c(10,10))
-s.class(pca1$li, as.factor(biofuc[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-8,2,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1, cex=1)
-boxplot(pca1$li$Axis1~as.factor(biofuc[,2]), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(biofuc[,2]), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(biofuc[,2]), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-
-#MANOVA
-rem<-findCorrelation(cor(mergedtable[,-1]), cutoff = .70, verbose = F, names=F)
-manovabio<-manova(as.matrix(mergedtable[,-c(1,rem+1)]) ~ as.factor(mergedtable[,1]) )
-summary(manovabio)
-
-
-
-
-#Soil
-soilgrids.r <- REST.SoilGrids(c("ORCDRC","PHIHOX"))
-pHbrot <- over(soilgrids.r, brot_cleaned)
-
-#pH
-library(lattice)
-library(aqp)
-data(soil.legends)
-plot(density(na.omit(ovfuc$PHIHOX.M.sd1)/10), col=3)
-lines(density(na.omit(ovinc$PHIHOX.M.sd1)/10), col=2)
-
-PHIHOX.range = range(soil.legends[["PHIHOX"]]$MIN, soil.legends[["PHIHOX"]]$MAX)
-
diff --git a/Nichoecologico_Dianthus.R b/Nichoecologico_Dianthus.R
deleted file mode 100644
index 47402ab..0000000
--- a/Nichoecologico_Dianthus.R
+++ /dev/null
@@ -1,145 +0,0 @@
-#R script Dactylorhiza Niche model divergence
-# This script downloads species coordinates from gbif. 
-# it clips the pH information from SoilGrids1km
-# it clip the wordclim data
-
-library(dismo)
-library('rgbif')
-library(rgdal)
-library(raster)
-library(gtools)
-library(rjson)
-library(sp)
-library(GSIF)
-library(plotrix)
-
-
-#GBIF data
-Dfucloc<-occ_search(scientificName = "Dianthus broteri",hasGeospatialIssue=FALSE,fields='minimal', limit=50000 ,hasCoordinate=TRUE)
-gbifmap(Dfucloc[[3]])
-dups <- duplicated(Dfucloc[[3]][,c(3,4)])
-fuc_cleaned<- Dfucloc[[3]][dups, ]
-fuc_cleaned
-coofuc<-fuc_cleaned[,c(4,3)]
-coordinates(coofuc)<- ~decimalLongitude+ decimalLatitude
-crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string(coofuc) <- crs.geo 
-
-data.shape<-readOGR(dsn="/home/fbalao/poligonodianthus.shp", layer="poligonodianthus")
-#Own data
-dianthus<-read.table("/home/fbalao/Datos/Proyectos/ProyectoExcelencia2016/Figuras/PCA_NichoEcologico_Dianthusbroteri/coordinatespopulationDianthus.csv",header=T)
-coordinates(dianthus)<- ~Longitude+ Latitude
-crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string(dianthus) <- crs.geo 
-
-
-#PCA
-data<-read.table("/home/fbalao/Datos/R/Analisis/Ecologicalniche/Datosbioclim.txt",header=T)
-attach(data)
-clima<-cbind(bio1,bio2,bio3,bio4,bio5,bio6,bio7,bio8,bio9,bio10,bio11,bio12,bio13,bio14,bio15,bio16,bio17,bio18,bio19)
-library(ade4)
-pca1 <- dudi.pca(clima, scann = FALSE, nf = 3)
-scatter(pca1, clab.row = 0, posieig = "none", add.plot=T, clab.col = 0.8, col=2)
-gcol<-c("blue", "red", "green", "yellow")
-s.class(pca1$li, as.factor(data[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-7,-1,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1.8)
-
-# MANOVA on PC Axis y plidy levels
-
-summary(manova(as.matrix(pca1$li)~as.factor(ploidy)))
-
-boxplot(pca1$li$Axis1~as.factor(ploidy), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(ploidy), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(ploidy), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-#Figura completa
-
-biolegend<-read.table("/home/fbalao/Datos/Proyectos/ProyectoExcelencia2016/Figuras/PCA_NichoEcologico_Dianthusbroteri/biolab.txt", sep="@")
-layout(matrix(c(1,2,2,3,4,5),nrow=2, byrow =T ))
-plot.new()
-addtable2plot(-0.7,0,biolegend,display.rownames=F,hlines=F, cex=0.5, display.colnames = F, bty="n")
-scatter(pca1, clab.row = 0, posieig = "none",clab.col = 1, col=2, xlim=c(10,10))
-s.class(pca1$li, as.factor(data[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-6.8,-0.1,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1, cex=0.7)
-boxplot(pca1$li$Axis1~as.factor(ploidy), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(ploidy), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(ploidy), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-# Kernel plots
-# plot(density(data$bio5[data$ploidy==2]), col="blue", xlim=c(200,500), ylim=c(0,0.25))
-# lines(density(data$bio5[data$ploidy==4]), col="red")
-# lines(density(data$bio5[data$ploidy==6]), col="green")
-# lines(density(data$bio5[data$ploidy==12]), col="yellow")
-
-#Soil
-soilgrids.r <- REST.SoilGrids(c("ORCDRC","PHIHOX"))
-ovfuc <- over(soilgrids.r, coofuc)
-ovinc <- over(soilgrids.r, cooinc)
-#pH
-library(lattice)
-library(aqp)
-data(soil.legends)
-plot(density(na.omit(ovfuc$PHIHOX.M.sd1)/10), col=3)
-lines(density(na.omit(ovinc$PHIHOX.M.sd1)/10), col=2)
-
-PHIHOX.range = range(soil.legends[["PHIHOX"]]$MIN, soil.legends[["PHIHOX"]]$MAX)
-
-
-
-#BIOCLIM DATA
-list.ras <- mixedsort(list.files("/home/fbalao/Datos/Dactylorhiza/BDGeoDact/wc2-5/", full.names = T, pattern = ".bil"))
-bio <- stack(list.ras)
-bio.brick <- brick(bio)
-newext <- drawExtent()
-bio.c<-crop(bio, newext)
-projection(bio.c)<-crs.geo
-
-prevfuc<-extract(bio.brick,coofuc)
-previnc<-extract(bio.brick,cooinc)
-
-biofuc<-cbind(as.data.frame(coofuc),prevfuc)
-bioinc<-cbind(as.data.frame(cooinc),previnc)
-bioinc<-na.omit(bioinc)
-library('caret')
-rem<-findCorrelation(cor(biofuc[,-c(1,2)]), cutoff = .90, verbose = F, names=F)
-rem<-rem+2
-biofuc.c<-biofuc[, -c(1,2,rem)]
-
-pca <- prcomp(biofuc.c)
-plot(pca$x[,1:2], pch=16, main=biofuc.c[1,"Species"]) # 2 first principal components
-z
-
-merged<-rbind(biofuc,bioinc)
-sp<-c(rep('fuc',348), rep('inc',924))
-merged<-cbind(sp,merged)
-mergedtable<-merged[,-c(2,3)]
-mergedtable
-
-#PCA on the climatic variables
-library(ade4)
-pca1 <- dudi.pca(mergedtable[,-c(1,rem+1)], scann = FALSE, nf = 3)
-barplot(pca1$eig)
-gcol = c("blue", "red", "green", "red")
-
-# s.class(dfxy = pca1$li, fac = as.factor(biofuc[,2]), col = gcol, xax = 1, yax = 2)
-# s.corcircle(pca1$co, xax = 1, yax = 2)
-# scatter(pca1, clab.row = 0, posieig = "none")
-# s.class(pca1$li, mergedtable[,1], col = gcol, add.plot = TRUE, cstar = 1.5, clabel = 0, cellipse = 1.5)
-
-
-biolegend<-read.table("/home/fbalao/Datos/Proyectos/ProyectoExcelencia2016/Figuras/PCA_NichoEcologico_Dianthusbroteri/biolab.txt", sep="@")
-layout(matrix(c(1,2,2,3,4,5),nrow=2, byrow =T ))
-plot.new()
-addtable2plot(-0,0,biolegend,display.rownames=F,hlines=F, cex=0.5, display.colnames = F, bty="n")
-scatter(pca1, clab.row = 0, posieig = "none",clab.col = 1, col=2, xlim=c(10,10))
-s.class(pca1$li, as.factor(biofuc[,2]), col = gcol, add.plot = TRUE, cstar = 0, clabel = 1, cellipse = 1.5, cpoint = 2)
-legend(-8,2,legend = c("2x","4x","6x", "12x"),  pch=21, pt.bg = gcol, pt.cex=1, cex=1)
-boxplot(pca1$li$Axis1~as.factor(biofuc[,2]), col=gcol, ylab="Componente 1",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis2~as.factor(biofuc[,2]), col=gcol, ylab="Componente 2",xlab="Nivel de ploidía")
-boxplot(pca1$li$Axis3~as.factor(biofuc[,2]), col=gcol,  ylab="Componente 3",xlab="Nivel de ploidía")
-
-
-#MANOVA
-rem<-findCorrelation(cor(mergedtable[,-1]), cutoff = .70, verbose = F, names=F)
-manovabio<-manova(as.matrix(mergedtable[,-c(1,rem+1)]) ~ as.factor(mergedtable[,1]) )
-summary(manovabio)
diff --git a/PhylogeneticAnalyses.R b/PhylogeneticAnalyses.R
new file mode 100644
index 0000000..21c599c
--- /dev/null
+++ b/PhylogeneticAnalyses.R
@@ -0,0 +1,116 @@
+###Phylogenetic signal
+library(phytools)
+library(phangorn)
+library(ape)
+
+###Loading tree and PCA values
+tree<-read.tree("dianthustree25.nw")
+tree2<-root(tree,outgroup = "BArportel")
+tree2<-chronos(tree2)
+tree2<-multi2di(tree2)
+plot(tree2)
+
+###PCA order
+orden<-tree2$tip.label
+PCAvalues<-read.csv2("scoresdata_gbif.csv", row.names = 4)
+PCAphylo<-PCAvalues[orden,]
+
+PCA1<-PCAphylo$Axis1
+names(PCA1)<-tree2$tip.label
+PCA2<-PCAphylo$Axis2
+names(PCA2)<-tree2$tip.label
+
+###Blomberg's K
+phylosig(tree2,PCAphylo[,1],method="K",test=TRUE, nsim=100000)
+phylosig(tree2,PCAphylo[,2],method="K",test=TRUE, nsim=100000)
+
+########################################################################
+
+###Phenograms
+PCA1<-PCAphylo$Axis1
+names(PCA1)<-tree2$tip.label
+PCA2<-PCAphylo$Axis2
+names(PCA2)<-tree2$tip.label
+
+###Phylogeny with ploidy levels mapped
+Ploidy<-PCAphylo$Ploidy #ploidy level assign
+names(Ploidy)<-tree2$tip.label
+trees <- make.simmap(tree2, Ploidy, model = "SYM", nsim=1)
+plot(trees)
+
+###Phylogeny with ploidy levels and lineages mapped
+Ploidy2<-factor(Ploidy, levels=c("2x","4x","4xn","6x","12x"))
+Ploidy2[20:25]<-"4xn"
+trees2 <- make.simmap(tree2, Ploidy2, model = "ER", nsim=1)
+plot(trees2)
+
+###Figure
+pdf(file="Phenograms.pdf", height = 10, width = 6, pagecentre=T)
+par(mfcol=c(2,1), mai = c(0.8, 0.8, 0.1, 0.1))
+phenogram(trees,PCA1, xlab="Relative time", ylab="Enviromental Axis 1", fsize=0.8,
+          lwd=3, ylim=c(-3.3,2.8),cex.lab=2)
+phenogram(trees,PCA2, xlab="Relative time", ylab="Enviromental Axis 1", fsize=0.8,
+          lwd=3, ylim=c(-7,-0.7))
+
+dev.off()
+
+
+#################################################################################
+
+###Niche evolution models
+
+library(mvMORPH)
+
+lambdaTree<-function(tree,lambda){
+  ii<-which(tree$edge[,2]>length(tree$tip.label))
+  H1<-nodeHeights(tree)
+  tree$edge.length[ii]<-lambda*tree$edge.length[ii]
+  H2<-nodeHeights(tree)
+  tree$edge.length[-ii]<-tree$edge.length[-ii]+     H1[-ii,2]-H2[-ii,2]
+  tree
+}
+
+treestar<-lambdaTree(tree, 0)
+
+
+###Null Model
+trait1_white<-mvBM(treestar, PCAphylo[,1], model="BM1")
+trait2_white<-mvBM(treestar, PCAphylo[,2], model="BM1")
+
+###Brownian model
+trait1_BM<-mvBM(trees, PCAphylo[,1], model="BM1")
+trait2_BM<-mvBM(trees, PCAphylo[,2], model="BM1")
+
+###OU model 1 optimum
+trait1_OU1<-mvOU(trees, PCAphylo[,1], model="OU1")
+trait2_OU1<-mvOU(trees, PCAphylo[,2], model="OU1")
+
+###OU model 4 optima
+trait1_OUM4<-mvOU(trees, PCAphylo[,1], model="OUM")
+trait2_OUM4<-mvOU(trees, PCAphylo[,2], model="OUM")
+
+###OU model 5 optima
+trait1_OUM5<-mvOU(trees2, PCAphylo[,1], model="OUM")
+trait2_OUM5<-mvOU(trees2, PCAphylo[,2], model="OUM")
+
+
+###Models comparison for PCA1
+resultsPC1<-list(trait1_white,trait1_BM,trait1_OU1,trait1_OUM4,trait1_OUM5)
+
+mvMORPH::aicw(resultsPC1, aicc=F)
+
+###Models comparison for PCA2
+resultsPC2<-list(trait2_white,trait2_BM,trait2_OU1,trait2_OUM4,trait2_OUM5)
+
+mvMORPH::aicw(resultsPC2, aicc=F)
+
+
+###Multivariate enviromental evolution models (both PCA axes)
+Mwhite<-mvBM(treestar, PCAphylo[,1:2], model="BM1")
+MBM1<-mvBM(trees, PCAphylo[,1:2], model="BM1")
+MOU1<-mvOU(trees, PCAphylo[,1:2], model="OU1")
+MOU4<-mvOU(trees, PCAphylo[,1:2], model="OUM")
+MOU5<-mvOU(trees2, PCAphylo[,1:2], model="OUM")
+
+results<-list(Mwhite,MBM1,MOU1,MOU4 ,MOU5)
+aicw(results)
diff --git a/SCRIPT_NICHE_GBIF.R b/SCRIPT_NICHE_GBIF.R
deleted file mode 100644
index 6ffec52..0000000
--- a/SCRIPT_NICHE_GBIF.R
+++ /dev/null
@@ -1,453 +0,0 @@
-library (dismo)
-library ('rgbif')
-library (rgdal)
-library (raster)
-library (gtools)
-library (rgeos)
-library (rjson)
-library (sp)
-library (GSIF)
-library (plotrix)
-library (ggmap)
-library (spatstat)
-library (ecospat)
-library (ggbiplot)
-library (spatialEco)
-library (fmsb)
-library (Hmisc)
-library (devtools)
-library (ENMTools)
-library (vioplot)
-
-
-#DATOS DE GBIF
-# Dbroterigbif <- occ_search (scientificName = "Dianthus broteri",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
-# Dinoxgbif <- occ_search (scientificName = "Dianthus inoxianus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
-# Dhinoxgbif <- occ_search (scientificName = "Dianthus hinoxianus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
-# Dvalengbif <- occ_search (scientificName = "Dianthus valentinus",hasGeospatialIssue =FALSE, limit=50000, fields = c("name", "decimalLatitude", "decimalLongitude", "coordinateUncertaintyInMeters", "year", "eventDate", "locality", "occurrenceRemarks", "recordedBy"), hasCoordinate=TRUE, basisOfRecord = "PRESERVED_SPECIMEN", return = "data")
-# 
-# 
-# #LIMPIEZA PARA QUEDARNOS CON UBICACIONES CON EXACTITUD (A PARTIR DE 2005)  
-# Dbroterigbif <- Dbroterigbif [Dbroterigbif$year>2004,]
-# Dbroterigbif <- Dbroterigbif [-c(129:189),]
-# Dhinoxgbif <- Dhinoxgbif [Dhinoxgbif$year>2004,]
-# Dhinoxgbif <- Dhinoxgbif [1,]
-# Dvalengbif <- Dvalengbif [1,]
-#   
-# Dbroterigbif <- Dbroterigbif [, c(1:3)]
-# Dinoxgbif <- Dinoxgbif [, c(1:3)]
-# Dhinoxgbif <- Dhinoxgbif [, c(1:3)]
-# Dvalengbif <- Dvalengbif [, c(1:3)]
-# 
-# todo <- rbind (Dbroterigbif,Dinoxgbif,Dhinoxgbif,Dvalengbif)
-# todo$name <- "Dianthus broteri"
-# gbifmap (todo, mapdatabase = "world", region = "Spain")
-# 
-# 
-# #LIMPIEZA DE PUNTOS DUPLICADOS
-# todo_cleaned <- unique(todo)
-# todo_cleaned <- todo_cleaned [,-1]
-# 
-# 
-# #MAPA
-# e <- extent (-10,3.5,35.5,44)
-# coordinates(todo_cleaned) <- ~decimalLongitude+ decimalLatitude
-# crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-# proj4string (todo_cleaned) <- crs.geo 
-# map <- plot (gmap (e, type = "satellite"))
-# points <- points (Mercator(todo_cleaned), col = "red", pch=20, cex = 1.5)
-# 
-# 
-# # Quitamos los puntos en el mismo km^2
-# r <- raster(todo_cleaned)
-# res(r) <- 0.008333333
-# r <- extend(r, extent(r)+1)
-# cooDbroterigbif <- as.data.frame(gridSample(todo_cleaned, r, n=1))
-# 
-# 
-# #===============PUNTOS EN DOÑANA (4x y 12x) DE FRAN===========#
-# 
-# e1 <- extent (-7,-6.4,36.8,37.5)
-# 
-# dianthuspops <- read.table ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/DianthusPoints.txt", header = T, sep = ",")
-# dianthuspops <- dianthuspops [,-c(1,4)]
-# dianthuspops <- dianthuspops [,c(2,1)]
-# colnames (dianthuspops) <- c("long", "lat")
-# 
-# coordinates(dianthuspops) <- ~long+ lat
-# proj4string(dianthuspops) <- crs.geo
-# r1 <- raster(dianthuspops)
-# res(r1) <- 0.008333333
-# r1 <- extend(r1, extent(r1)+1)
-# dianthuspops_fran <- as.data.frame(gridSample(dianthuspops, r1, n=1))
-# 
-# coordinates(dianthuspops_fran) <- ~long+ lat
-# proj4string(dianthuspops_fran) <- crs.geo
-# plot(gmap(e1, type = "satellite"))
-# points(Mercator(dianthuspops_fran), col="red" , pch=20, cex=1)
-# 
-# 
-# #===============PUNTOS EN DOÑANA (4x y 12x) DE FRAN===========#
-# 
-# dianthuspops_fran <- as.data.frame (dianthuspops_fran)
-# colnames(cooDbroterigbif) <- colnames(dianthuspops_fran)
-# cooDbroterigbif.def <- rbind (cooDbroterigbif, dianthuspops_fran)
-# cooDbroterigbif.def2 <- cooDbroterigbif.def
-# 
-# coordinates(cooDbroterigbif.def2) <- ~long+ lat
-# proj4string (cooDbroterigbif.def2) <- crs.geo 
-# 
-# r <- raster(cooDbroterigbif.def2)
-# res(r) <- 0.008333333
-# r <- extend(r, extent(r)+1)
-# cooDbroterigbif.def.data <- as.data.frame(gridSample(cooDbroterigbif.def2, r, n=1))
-# 
-# cooDbroterigbif.def.data <- cbind (cooDbroterigbif.def.data, 1)
-# colnames(cooDbroterigbif.def.data)[3]<-"ploidy"
-# 
-# #Estimacion del nivel de ploidia por la ubicacion geografica (coordenadas)
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38 & cooDbroterigbif.def.data$long > -9.5 & cooDbroterigbif.def.data$lat < 39 & cooDbroterigbif.def.data$long < -8] <- "4x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 37 & cooDbroterigbif.def.data$long > -9 & cooDbroterigbif.def.data$lat < 37.5 & cooDbroterigbif.def.data$long < -7.5] <- "2x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36.9 & cooDbroterigbif.def.data$long > -7.4 & cooDbroterigbif.def.data$lat < 37.6 & cooDbroterigbif.def.data$long < -6] <- "12x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36 & cooDbroterigbif.def.data$long > -6.3 & cooDbroterigbif.def.data$lat < 37.14 & cooDbroterigbif.def.data$long < -4.72] <- "4x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 36.5 & cooDbroterigbif.def.data$long > -4.7 & cooDbroterigbif.def.data$lat < 38.1 & cooDbroterigbif.def.data$long < -1.67] <- "2x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 37.3 & cooDbroterigbif.def.data$long > -1.7 & cooDbroterigbif.def.data$lat < 38.3 & cooDbroterigbif.def.data$long < -0.5] <- "6x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38.2 & cooDbroterigbif.def.data$long > -2.33 & cooDbroterigbif.def.data$lat < 38.83 & cooDbroterigbif.def.data$long < 0.33] <- "6x"
-# cooDbroterigbif.def.data$ploidy[cooDbroterigbif.def.data$lat > 38.83 & cooDbroterigbif.def.data$long > -1.7 & cooDbroterigbif.def.data$lat < 40.9 & cooDbroterigbif.def.data$long < 0.38] <- "4x"
-# cooDbroterigbif.def.data$ploidy[c(58,62:63,120:123,125:130)]<- "4x"
-# cooDbroterigbif.def.data <- cooDbroterigbif.def.data[cooDbroterigbif.def.data$ploidy!="1",]
-# 
-# 
-# ploidy<-factor(cooDbroterigbif.def.data$ploidy,levels = c("2x","4x","6x","12x"),ordered = TRUE)
-# 
-# coordinates(cooDbroterigbif.def.data) <- ~long+ lat
-# proj4string(cooDbroterigbif.def.data) <- crs.geo
-# plot(gmap(e, type = "satellite"))
-# points(Mercator(cooDbroterigbif.def.data), col=ploidy, pch=20, cex=1)
-# 
-# plot(gmap(e1, type = "satellite"))
-# points(Mercator(cooDbroterigbif.def.data), col=ploidy, pch=20, cex=1)
-# 
-# 
-# cooDbroterigbif.def.vars <- as.data.frame (cooDbroterigbif.def.data)
-# cooDbroterigbif.def.vars <- cooDbroterigbif.def.vars [,c(3,1,2)]
-# coordinates(cooDbroterigbif.def.vars) <- ~long+ lat
-# proj4string(cooDbroterigbif.def.vars) <- crs.geo
-
-
-#carga de variables predictoras y union con mismos limites (chelsa, envirem, altitud, SoilGrids)
-#extraccion de datos de las variables predictoras en las poblaciones
-
-cooDbroterigbif.def.vars <- read.csv2 ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/resultados_analisis/occurrences_gbif.csv")
-ploidy<-factor(cooDbroterigbif.def.vars$ploidy,levels = c("2x","4x","6x","12x"),ordered = TRUE)
-
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-coordinates(cooDbroterigbif.def.vars) <- ~long+ lat
-proj4string(cooDbroterigbif.def.vars) <- crs.geo
-
-e <- extent (-10,4.5,35.5,44)
-
-chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chelsafiles)
-che.c <- crop (chelsa,e)
-
-enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (enviremfiles)
-env.c <- crop (envirem, e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-variables <- stack (che.c, env.c, alt.m)
-names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","climaticMoistureIndex","continentality","growingDegDays5","maxTempColdest","minTempWarmest","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","topoWet", "elevation")
-
-soilgrids<-extract.list(cooDbroterigbif.def.vars, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas", ID = "ploidy")
-colnames (soilgrids) <- c("ploidy","AWCh2","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT")
-soilgrids$ploidy<-factor(soilgrids$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="1","clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="2","silty clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="3","sandy clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="4","clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="5","silty clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="6","sandy clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="7","loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="8","silty loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="9","sandy loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="10","silt")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="11","loamy sand")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="12","sand")
-# soilgrids$TEXMHT<-factor(soilgrids$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
-# soilgrids<-cbind(soilgrids,apply(soilgrids[,c(2:4)], 1, mean))
-# soilgrids<-soilgrids[,-c(1:4)]
-colnames(soilgrids)[2]<-"AWC"
-
-
-presvals <- extract (variables, cooDbroterigbif.def.vars)
-presvals <- cbind (ploidy, presvals, soilgrids) 
-presvals$PHIHOX <- presvals$PHIHOX/10
-presvals <- presvals [,-34]
-
-#calculo del tri a partir de altitud con libreria spatialEco
-
-tri.ext <- tri(alt.m)
-projection(tri.ext) <- crs.geo 
-trivalues<-extract(tri.ext,cooDbroterigbif.def.vars)
-
-
-#sustitucion de los valores tri por los del dataframe con NAs
-presvals <- cbind (presvals, trivalues)
-presvals <- presvals[,c(1:31,40,32:39)]
-colnames(presvals)[32] <- "tri"
-presvals2 <- na.omit (presvals)
-ploidy <- as.data.frame(presvals2[,1])
-colnames (ploidy) <- "ploidy"
-
-
-#analisis para descartar variables muy correlacionadas
-#PCA de puntos de presencia con variables seleccionadas
-
-presvals.pca <- presvals2[,-1]
-selected <- vif_func(presvals.pca) 
-presvals.pca.2 <- presvals.pca[,c(selected)]
-
-ploidy <- ploidy$ploidy
-ploidy <- factor (ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-
-pca <- prcomp(presvals.pca.2, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#PCA con todas las variables
-pca <- prcomp(presvals.pca, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#===============GENERAL BACKGROUND=============#
-
-dbroteridata <- as.data.frame(cooDbroterigbif.def.vars)
-dbroteridata <- dbroteridata [,-1]
-presvalsdata <- cbind (dbroteridata, presvals)
-presvalsdata <- na.omit (presvalsdata)
-presvals2 <- presvalsdata
-coordinates(presvals2) <- ~long+ lat
-proj4string(presvals2) <- crs.geo
-
-backgroundcoord <- presvalsdata [,c(1,2)]
-coordinates(backgroundcoord) <- ~long+ lat
-proj4string(backgroundcoord) <- crs.geo
-
-
-mask <- raster(backgroundcoord)
-res(mask) <- 0.008333333
-x <- circles(backgroundcoord, d=150000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-pol <- gUnaryUnion(x@polygons)
-samp <- spsample(pol, 10000, type='random', iter=25)
-extent(mask) <- extent(pol) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cells <- cellFromXY(mask, samp)
-length(cells)
-cells <- unique(cells)
-length(cells)
-xy <- xyFromCell(mask, cells)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xy <- as.data.frame(xy)
-colnames(xy) <- c("long", "lat")
-
-# Quitamos los puntos en el mismo km^2
-coordinates(xy) <- ~long+ lat
-proj4string(xy) <- crs.geo
-r <- raster(xy)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-backcoord_sel <- as.data.frame(gridSample(xy, r, n=1))
-coordinates(backcoord_sel) <- ~long+ lat
-proj4string(backcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-backgroundclim <- extract(variables,backcoord_sel)
-backgroundsoil <- extract.list(backcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas", ID = "ploidy")
-backgrounddat <- cbind("background",as.data.frame(backcoord_sel),backgroundclim, backgroundsoil)
-backgrounddat <- backgrounddat[,-36]
-
-coordinates(backgrounddat) <- ~long+ lat
-proj4string(backgrounddat) <- crs.geo
-trivalues1 <- extract(tri.ext,backgrounddat)
-backgrounddat <- as.data.frame(backgrounddat)
-backgrounddat <- cbind (backgrounddat, trivalues1)
-backgrounddat <- backgrounddat[,c(1:33,42,34:41)]
-colnames(backgrounddat)[34] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-backgrounddat.c <- na.omit(backgrounddat)
-backgrounddat.c <- backgrounddat.c[,c(2,3,1,4:42)]
-colnames(backgrounddat.c) <- colnames(presvalsdata)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(backgrounddat.c) <- ~long+ lat
-proj4string(backgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(backgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvals2), col=presvals2$ploidy, pch=20, cex=1)
-
-#===============GENERAL BACKGROUND=============#
-
-#===============PCA BACKGROUND=============#
-
-todo <- rbind (presvalsdata, as.data.frame(backgrounddat.c))
-
-todo.pca <- todo[,-c(1:3)]
-selected2<-vif_func(todo.pca)
-todo.pca.2 <- todo.pca[,c(selected2)]
-
-todoploidy <- factor (todo$ploidy, levels = c("2x", "4x", "6x", "12x", "background"), ordered = TRUE)
-w<-c(rep(0,nrow(presvalsdata)),rep(1,nrow(as.data.frame(backgrounddat.c))))
-
-pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
-gcol = c("blue", "red", "green", "purple", "black")
-s.label(pcaback$li, clabel = 0.1)
-scatter(pcaback, clab.row = 0, posieig = "none", cex=0.1, clab.col = 0.5)
-s.class(pcaback$li, todoploidy, col = gcol, add.plot = TRUE, cstar = 0, clabel = 0, cellipse = 1.5, pch = 16)
-legend (8.5,-3,c("2x", "4x", "6x", "12x","Background"), col = gcol, pch =19, text.width = 1.8, y.intersp = 0.5, cex = 0.8)
-
-
-#===============ECOSPAT=============#
-
-row.di<-which(todo[,3] == "2x")
-row.te<-which(todo[,3] == "4x")
-row.he<-which(todo[,3] == "6x")
-row.do<-which(todo[,3] == "12x")
-
-
-scores.clim<- pcaback$li 
-scores.di<- pcaback$li[row.di,]		
-scores.te<- pcaback$li[row.te,]	
-scores.he<- pcaback$li[row.he,]	
-scores.do<- pcaback$li[row.do,]	
-
-
-zdi<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.di, R=500)
-zte<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.te, R=500)
-zhe<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.he, R=500)
-zdo<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.do, R=500)
-
-ecospat.plot.niche (zdi)
-ecospat.plot.niche (zte)
-ecospat.plot.niche (zhe)
-ecospat.plot.niche (zdo)
-ecospat.plot.niche.dyn (zdi, zte, quant = 0.75)
-
-
-#EQUIVALENCY TEST
-
-equivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 100, alternative = "lower")
-equivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 100, alternative = "lower")
-equivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 100, alternative = "lower")
-equivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 100, alternative = "lower")
-equivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 100, alternative = "lower")
-equivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 100, alternative = "lower")
-
-
-#OVERLAP TEST
-
-overlap.test.dite<-ecospat.niche.overlap (zdi, zte, cor=FALSE)
-overlap.test.dihe<-ecospat.niche.overlap (zdi, zhe, cor=FALSE)
-overlap.test.dido<-ecospat.niche.overlap (zdi, zdo, cor=FALSE)
-overlap.test.tehe<-ecospat.niche.overlap (zte, zhe, cor=FALSE)
-overlap.test.tedo<-ecospat.niche.overlap (zte, zdo, cor=FALSE)
-overlap.test.hedo<-ecospat.niche.overlap (zhe, zdo, cor=FALSE)
-
-
-#SIMILARITY TEST
-
-similarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 100, alternative = "greater")
-similarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 100, alternative = "greater")
-similarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 100, alternative = "greater")
-similarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 100, alternative = "greater")
-similarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 100, alternative = "greater")
-similarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 100, alternative = "greater")
-similarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 100, alternative = "greater")
-similarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 100, alternative = "greater")
-similarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 100, alternative = "greater")
-similarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 100, alternative = "greater")
-similarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 100, alternative = "greater")
-similarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 100, alternative = "greater")
-
-
-# TABLA
-
-overlap<-rbind(overlap.test.dite$D,overlap.test.dihe$D, overlap.test.dido$D,overlap.test.tehe$D, overlap.test.tedo$D, overlap.test.hedo$D)
-
-similarityab<-rbind(similarity.testdite$p.D, similarity.testdihe$p.D, similarity.testdido$p.D, similarity.testtehe$p.D, similarity.testtedo$p.D, similarity.testhedo$p.D)
-similarityba<-rbind(similarity.testtedi$p.D, similarity.testhedi$p.D, similarity.testdodi$p.D, similarity.testhete$p.D, similarity.testdote$p.D, similarity.testdohe$p.D)
-
-equivalency<-rbind(equivalency.test.dite$p.D, equivalency.test.dihe$p.D, equivalency.test.dido$p.D, equivalency.test.tehe$p.D, equivalency.test.tedo$p.D, equivalency.test.hedo$p.D)
-
-tablaresul<-data.frame(overlap,similarityab,similarityba,equivalency)
-
-write.table (tablaresul, "results_gbif_vif.txt", sep = "\t")
-
-
-#VIOPLOTS
-
-vioplot_fun (scores.di$Axis1, scores.te$Axis1, scores.he$Axis1, scores.do$Axis1, names = c("2x","4x","6x","12x"), col = c("green", "red", "blue", "purple"))
-vioplot_fun (scores.di$Axis2, scores.te$Axis2, scores.he$Axis2, scores.do$Axis2, names = c("2x","4x","6x","12x"), col = c("green", "red", "blue", "purple"))
-
-
-#DENSITY PLOTS
-
-ggplot(presvalsdata, aes(x = PETWettestQuarter, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = AWC, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = PHIHOX, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = aridityIndexThornthwaite, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = BLDFIE, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = SNDPPT, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = ORCDRC, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio3, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio5, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio8, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio9, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio12, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio18, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-
-#NICHE BREADTH
-
-raster.breadth (zdi$w)
-raster.breadth (zte$w)
-raster.breadth (zhe$w)
-raster.breadth (zdo$w)
diff --git a/SCRIPT_NICHE_POPS.R b/SCRIPT_NICHE_POPS.R
deleted file mode 100644
index c4125a0..0000000
--- a/SCRIPT_NICHE_POPS.R
+++ /dev/null
@@ -1,358 +0,0 @@
-library (dismo)
-library (raster)
-# library (rJava)
-# library (rgdal)
-library (rgeos)
-# library (rgbif)
-# library (rjson)
-library (gtools)
-# library (maps)
-# library (ggmap)
-# library (fuzzySim)
-# library (ade4)
-# library (pcaMethods)
-library (ecospat)
-library (sp)
-library (GSIF)
-# library (caret)
-# library (RCurl)
-# library (gdalUtils)
-# library (plotKML)
-# library (XML)
-# library (lattice)
-# library (aqp)
-# library (soiltexture)
-library (ggplot2)
-library (ggbiplot)
-# library (psych)
-library (spatialEco)
-library (fmsb)
-library (Hmisc)
-library (devtools)
-library (ENMTools)
-library (animation)
-library (vioplot)
-
-
-#dataset de poblaciones con coordenadas
-
-dbroteri <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/Poblaciones_Nicho.csv", sep = ";", fileEncoding = "latin1", colClasses = c("factor", "factor", "numeric", "numeric"))
-ploidy <- dbroteri[,2]
-ploidy <- as.data.frame (ploidy)
-rownames (ploidy) <- c("Monte Clerigo","Sao Bras de Alportel","Zafarraya 1","Zafarraya 2","Orgiva","Laroles","Lliria","Troia","Comporta","Donana (Peladillo)","Albufera de Valencia","Albufera de Valencia 2","Alcublas","Azuebar","Sierra de Espadan","Chiclana","Ronda","Calblanque (Cabezo de la Fuente)","Socovos","Cartagena","San Miguel de Salinas","Penon de Ifach","Valverde del Camino","Moguer","Hinojos","Donana (Acebron)","Donana (Puntal)","Huertos del Batan","Isla Cristina (Las Palmeritas)")
-coordinates (dbroteri) <- ~long + lat
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string (dbroteri) <- crs.geo
-
-
-#carga de variables predictoras y union con mismos limites (chelsa, envirem, altitud, SoilGrids)
-#extraccion de datos de las variables predictoras en las poblaciones
-
-e <- extent (-10,4.5,35.5,44)
-
-chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chelsafiles)
-che.c <- crop (chelsa,e)
-
-enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (enviremfiles)
-env.c <- crop (envirem, e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-variables <- stack (che.c, env.c)
-names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","climaticMoistureIndex","continentality","growingDegDays5","maxTempColdest","minTempWarmest","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","topoWet")
-
-soilgrids<-extract.list(dbroteri, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas", ID = "ploidy")
-colnames (soilgrids) <- c("ploidy","AWCh2","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT")
-soilgrids$ploidy<-factor(soilgrids$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="1","clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="2","silty clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="3","sandy clay")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="4","clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="5","silty clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="6","sandy clay loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="7","loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="8","silty loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="9","sandy loam")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="10","silt")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="11","loamy sand")
-# soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="12","sand")
-# soilgrids$TEXMHT<-factor(soilgrids$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
-# soilgrids<-cbind(soilgrids,apply(soilgrids[,c(2:4)], 1, mean))
-# soilgrids<-soilgrids[,-c(1:4)]
-colnames(soilgrids)[2]<-"AWC"
-
-
-presvals <- extract (variables, dbroteri)
-presvals <- cbind (ploidy, presvals, soilgrids) 
-presvals$PHIHOX <- presvals$PHIHOX/10
-presvals <- presvals [,-33]
-
-#calculo del tri a partir de altitud con libreria spatialEco
-
-tri.ext <- tri(alt.m)
-projection(tri.ext) <- crs.geo 
-trivalues<-extract(tri.ext,dbroteri)
-
-
-#sustitucion de los valores tri por los del dataframe con NAs
-presvals <- cbind (presvals, trivalues)
-presvals <- presvals[,c(1:31,39,32:38)]
-colnames(presvals)[32] <- "tri"
-
-
-#analisis para descartar variables muy correlacionadas
-#PCA de puntos de presencia con variables seleccionadas
-
-presvals.pca <- presvals[,-1]
-selected <- vif_func(presvals.pca) 
-presvals.pca.2 <- presvals.pca[,c(selected)]
-
-# presvals.pca.corselect <- cbind (presvals.pca,1)
-# correlations <- corSelect (presvals.pca.corselect, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.75)
-# presvals.pca.2.corselect <- presvals.pca.corselect[,correlations$selected.var.cols]
-
-
-ploidy <- ploidy$ploidy
-ploidy <- factor (ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-
-pca <- prcomp(presvals.pca.2, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#PCA con todas las variables
-pca_all <- prcomp(presvals.pca, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#===============GENERAL BACKGROUND=============#
-
-dbroteridata <- as.data.frame(dbroteri)
-presvalsdata <- cbind (dbroteridata[,c(3,4)], presvals)
-presvalsdata <- presvalsdata [,c(2,1,3:41)]
-presvals2 <- presvalsdata
-coordinates(presvals2) <- ~long+ lat
-proj4string(presvals2) <- crs.geo
-
-backgroundcoord <- presvalsdata [,c(1,2)]
-coordinates(backgroundcoord) <- ~long+ lat
-proj4string(backgroundcoord) <- crs.geo
-
-
-mask <- raster(backgroundcoord)
-res(mask) <- 0.008333333
-x <- circles(backgroundcoord, d=150000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-pol <- gUnaryUnion(x@polygons)
-samp <- spsample(pol, 10000, type='random', iter=25)
-extent(mask) <- extent(pol) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cells <- cellFromXY(mask, samp)
-length(cells)
-cells <- unique(cells)
-length(cells)
-xy <- xyFromCell(mask, cells)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xy <- as.data.frame(xy)
-colnames(xy) <- c("long", "lat")
-
-# Quitamos los puntos en el mismo km^2
-coordinates(xy) <- ~long+ lat
-proj4string(xy) <- crs.geo
-r <- raster(xy)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-backcoord_sel <- as.data.frame(gridSample(xy, r, n=1))
-coordinates(backcoord_sel) <- ~long+ lat
-proj4string(backcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-backgroundclim <- extract(variables,backcoord_sel)
-backgroundsoil <- extract.list(backcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/soilgrids/capas", ID = "ploidy")
-backgrounddat <- cbind("background",as.data.frame(backcoord_sel),backgroundclim, backgroundsoil)
-backgrounddat <- backgrounddat[,-35]
-
-coordinates(backgrounddat) <- ~long+ lat
-proj4string(backgrounddat) <- crs.geo
-trivalues1 <- extract(tri.ext,backgrounddat)
-backgrounddat <- as.data.frame(backgrounddat)
-backgrounddat <- cbind (backgrounddat, trivalues1)
-backgrounddat <- backgrounddat[,c(1:33,41,34:40)]
-colnames(backgrounddat)[34] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-backgrounddat.c <- na.omit(backgrounddat)
-backgrounddat.c <- backgrounddat.c[,c(2,3,1,4:41)]
-colnames(backgrounddat.c) <- colnames(presvalsdata)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(backgrounddat.c) <- ~long+ lat
-proj4string(backgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(backgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvals2), col=presvals2$ploidy, pch=20, cex=1)
-
-#===============GENERAL BACKGROUND=============#
-
-#===============PCA BACKGROUND=============#
-
-todo <- rbind (presvalsdata, as.data.frame(backgrounddat.c))
-
-todo.pca <- todo[,-c(1:3)]
-selected2 <- vif_func(todo.pca)
-todo.pca.2 <- todo.pca[,c(selected2)]
-
-# presvals.pca.corselect2 <- cbind (todo.pca,1)
-# correlations2 <- corSelect (presvals.pca.corselect2, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.75)
-# presvals.pca.2.corselect2 <- presvals.pca.corselect2[,correlations2$selected.var.cols]
-
-todoploidy <- factor (todo$ploidy, levels = c("2x", "4x", "6x", "12x", "background"), ordered = TRUE)
-w<-c(rep(0,nrow(presvalsdata)),rep(1,nrow(as.data.frame(backgrounddat.c))))
-
-pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
-gcol = c("blue", "red", "green", "purple", "black")
-s.label(pcaback$li, clabel = 0.1)
-scatter(pcaback, clab.row = 0, posieig = "none", cex=0.1, clab.col = 0.5)
-s.class(pcaback$li, todoploidy, col = gcol, add.plot = TRUE, cstar = 0, clabel = 0, cellipse = 1.5, pch = 16)
-legend (-12,2,c("2x", "4x", "6x", "12x","Background"), col = gcol, pch =19, text.width = 1.8, y.intersp = 0.5, cex = 0.8)
-
-
-#===============ECOSPAT=============#
-
-row.di<-which(todo[,3] == "2x")
-row.te<-which(todo[,3] == "4x")
-row.he<-which(todo[,3] == "6x")
-row.do<-which(todo[,3] == "12x")
-
-
-scores.clim<- pcaback$li 
-scores.di<- pcaback$li[row.di,]		
-scores.te<- pcaback$li[row.te,]	
-scores.he<- pcaback$li[row.he,]	
-scores.do<- pcaback$li[row.do,]	
-
-
-zdi<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.di, R=500)
-zte<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.te, R=500)
-zhe<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.he, R=500)
-zdo<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.do, R=500)
-
-ecospat.plot.niche (zdi)
-ecospat.plot.niche (zte)
-ecospat.plot.niche (zhe)
-ecospat.plot.niche (zdo)
-ecospat.plot.niche.dyn (zdi, zte, quant = 0.75)
-
-
-#EQUIVALENCY TEST
-equivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 100, alternative = "lower")
-equivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 100, alternative = "lower")
-equivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 100, alternative = "lower")
-equivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 100, alternative = "lower")
-equivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 100, alternative = "lower")
-equivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 100, alternative = "lower")
-
-
-#OVERLAP TEST
-overlap.test.dite<-ecospat.niche.overlap (zdi, zte, cor=FALSE)
-overlap.test.dihe<-ecospat.niche.overlap (zdi, zhe, cor=FALSE)
-overlap.test.dido<-ecospat.niche.overlap (zdi, zdo, cor=FALSE)
-overlap.test.tehe<-ecospat.niche.overlap (zte, zhe, cor=FALSE)
-overlap.test.tedo<-ecospat.niche.overlap (zte, zdo, cor=FALSE)
-overlap.test.hedo<-ecospat.niche.overlap (zhe, zdo, cor=FALSE)
-
-
-#SIMILARITY TEST greater alternative hypothesis more similar than random
-similarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 100, alternative = "greater")
-similarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 100, alternative = "greater")
-similarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 100, alternative = "greater")
-similarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 100, alternative = "greater")
-similarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 100, alternative = "greater")
-similarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 100, alternative = "greater")
-similarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 100, alternative = "greater")
-similarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 100, alternative = "greater")
-similarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 100, alternative = "greater")
-similarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 100, alternative = "greater")
-similarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 100, alternative = "greater")
-similarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 100, alternative = "greater")
-
-
-# TABLA
-
-overlap<-rbind(overlap.test.dite$D,overlap.test.dihe$D, overlap.test.dido$D,overlap.test.tehe$D, overlap.test.tedo$D, overlap.test.hedo$D)
-
-similarityab<-rbind(similarity.testdite$p.D, similarity.testdihe$p.D, similarity.testdido$p.D, similarity.testtehe$p.D, similarity.testtedo$p.D, similarity.testhedo$p.D)
-similarityba<-rbind(similarity.testtedi$p.D, similarity.testhedi$p.D, similarity.testdodi$p.D, similarity.testhete$p.D, similarity.testdote$p.D, similarity.testdohe$p.D)
-
-equivalency<-rbind(equivalency.test.dite$p.D, equivalency.test.dihe$p.D, equivalency.test.dido$p.D, equivalency.test.tehe$p.D, equivalency.test.tedo$p.D, equivalency.test.hedo$p.D)
-
-tablaresul<-data.frame(overlap,similarityab,similarityba,equivalency)
-
-write.table (tablaresul, "results_pops_vif.txt", sep = "\t")
-
-
-#VIOPLOTS
-
-vioplot_fun (-scores.di$Axis1, -scores.te$Axis1, -scores.he$Axis1, -scores.do$Axis1, names = c("2x","4x","6x","12x"), col = c("green", "red", "blue", "purple"))
-vioplot_fun (-scores.di$Axis2, -scores.te$Axis2, -scores.he$Axis2, -scores.do$Axis2, names = c("2x","4x","6x","12x"), col = c("green", "red", "blue", "purple"))
-
-
-#DENSITY PLOTS
-
-ggplot(presvalsdata, aes(x = PETWettestQuarter, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = AWC, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = PHIHOX, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = aridityIndexThornthwaite, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = BLDFIE, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = SNDPPT, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = ORCDRC, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio3, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio5, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio8, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio9, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio12, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-ggplot(presvalsdata, aes(x = bio18, fill = ploidy)) + 
-  geom_density(alpha=0.3)
-
-
-#NICHE BREADTH
-raster.breadth (zdi$w)
-raster.breadth (zte$w)
-raster.breadth (zhe$w)
-raster.breadth (zdo$w)
-
diff --git a/divergence_test.R b/divergence_test.R
deleted file mode 100644
index 9c01c29..0000000
--- a/divergence_test.R
+++ /dev/null
@@ -1,47 +0,0 @@
-
-PC1dn<-t.test(scores.di[,1], scores.te[,1]) 
-PC1dn1<-t.test(scores.di[,1], scores.he[,1]) 
-PC1dn2<-t.test(scores.di[,1], scores.do[,1]) 
-PC1dn3<-t.test(scores.te[,1], scores.he[,1]) 
-PC1dn4<-t.test(scores.te[,1], scores.do[,1]) 
-PC1dn5<-t.test(scores.he[,1], scores.do[,1]) 
-
-PC2dn<-t.test(scores.di[,2], scores.te[,2]) 
-PC2dn1<-t.test(scores.di[,2], scores.he[,2]) 
-PC2dn2<-t.test(scores.di[,2], scores.do[,2]) 
-PC2dn3<-t.test(scores.te[,2], scores.he[,2]) 
-PC2dn4<-t.test(scores.te[,2], scores.do[,2]) 
-PC2dn5<-t.test(scores.he[,2], scores.do[,2]) 
-
-dbPC1<-vector()
-for (i in 1:9999){
-  dbiPC1<-mean(sample(scores.clim[,1],500, replace=T))
-  dbfPC1<-mean(sample(scores.clim[,1],500, replace=T))
-  dbPC1[i]<-dbfPC1-dbiPC1
-}
-
-dbPC2<-vector()
-for (i in 1:9999){
-  dbiPC2<-mean(sample(scores.clim[,2],500, replace=T))
-  dbfPC2<-mean(sample(scores.clim[,2],500, replace=T))
-  dbPC2[i]<-dbfPC2-dbiPC2
-}
-
-dn1<-mean(scores.di[,1])- mean(scores.te[,1])
-dn1_1<-mean(scores.di[,1])- mean(scores.he[,1])
-dn1_2<-mean(scores.di[,1])- mean(scores.do[,1])
-dn1_3<-mean(scores.te[,1])- mean(scores.he[,1])
-dn1_4<-mean(scores.te[,1])- mean(scores.do[,1])
-dn1_5<-mean(scores.he[,1])- mean(scores.do[,1])
-db1<- quantile(dbPC1,c(0.025, .975))
-
-dn2<-mean(scores.di[,2])- mean(scores.te[,2])
-dn2_1<-mean(scores.di[,2])- mean(scores.he[,2])
-dn2_2<-mean(scores.di[,2])- mean(scores.do[,2])
-dn2_3<-mean(scores.te[,2])- mean(scores.he[,2])
-dn2_4<-mean(scores.te[,2])- mean(scores.do[,2])
-dn2_5<-mean(scores.he[,2])- mean(scores.do[,2])
-db2<-quantile(dbPC2,c(0.025, .975))
-
-#dbIC= (-0.011362342  0.002173957)
-#daPC1=0.256439
diff --git a/envdianthus.Rproj b/envdianthus.Rproj
new file mode 100644
index 0000000..8e3c2eb
--- /dev/null
+++ b/envdianthus.Rproj
@@ -0,0 +1,13 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
diff --git a/maxent.R b/maxent.R
deleted file mode 100644
index ba2be56..0000000
--- a/maxent.R
+++ /dev/null
@@ -1,62 +0,0 @@
-library (dismo)
-library (ENMeval)
-
-dbroteri <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/Poblaciones_Nicho.csv", sep = ";", fileEncoding = "latin1", colClasses = c("factor", "factor", "numeric", "numeric"))
-coordinates (dbroteri) <- ~long + lat
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string (dbroteri) <- crs.geo
-
-dbroteri_back <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/background_pops.csv", sep = ";")
-coordinates (dbroteri_back) <- ~long + lat
-proj4string (dbroteri_back) <- crs.geo
-
-vars.stack <- stack ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/stack_zoon.grd")
-
-dbroteridata <- as.data.frame(dbroteri)
-dbroteridata <- dbroteridata [,c(4,3)]
-dbroteri_backdata <- as.data.frame (dbroteri_back)
-
-dbroteri_di <- dbroteridata [c(1:6),]
-dbroteri_te <- dbroteridata [c(7:17),]
-dbroteri_he <- dbroteridata [c(18:22),]
-dbroteri_do <- dbroteridata [c(23:29),]
-
-coordinates (dbroteri_di) <- ~long + lat
-proj4string (dbroteri_di) <- crs.geo
-coordinates (dbroteri_te) <- ~long + lat
-proj4string (dbroteri_te) <- crs.geo
-coordinates (dbroteri_he) <- ~long + lat
-proj4string (dbroteri_he) <- crs.geo
-coordinates (dbroteri_do) <- ~long + lat
-proj4string (dbroteri_do) <- crs.geo
-
-mod_maxent <- maxent (x=vars.stack,p=dbroteri,a=dbroteri_back,path="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/maxent", responsecurves = T)
-mod_predict <- predict (mod_maxent, vars.stack, c(-10,1,36,40.5))
-plot (mod_predict)
-
-mod_maxent_di <- maxent (x=vars.stack,p=dbroteri_di,a=dbroteri_back,path="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/maxent", responsecurves = T)
-mod_predict_di <- predict (mod_maxent_di, vars.stack, c(-10,1,36,40.5))
-plot (mod_predict_di)
-
-mod_maxent_te <- maxent (x=vars.stack,p=dbroteri_te,a=dbroteri_back,path="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/maxent", responsecurves = T)
-mod_predict_te <- predict (mod_maxent_te, vars.stack, c(-10,1,36,40.5))
-plot (mod_predict_te)
-
-mod_maxent_he <- maxent (x=vars.stack,p=dbroteri_he,a=dbroteri_back,path="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/maxent", responsecurves = T)
-mod_predict_he <- predict (mod_maxent_he, vars.stack, c(-10,1,36,40.5))
-plot (mod_predict_he)
-
-mod_maxent_do <- maxent (x=vars.stack,p=dbroteri_do,a=dbroteri_back,path="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/maxent", responsecurves = T)
-mod_predict_do <- predict (mod_maxent_do, vars.stack, c(-10,1,36,40.5))
-plot (mod_predict_do)
-
-dbroteri_di <- as.data.frame (dbroteri_di)
-dbroteri_te <- as.data.frame (dbroteri_te)
-dbroteri_he <- as.data.frame (dbroteri_he)
-dbroteri_do <- as.data.frame (dbroteri_do)
-
-evaluation <- ENMevaluate (dbroteridata, vars.stack, dbroteri_backdata, method = "jackknife", rasterPreds = T, bin.output=TRUE)
-evaluation_di <- ENMevaluate (dbroteri_di, vars.stack, dbroteri_backdata, method = "jackknife", rasterPreds = T, bin.output=TRUE)
-evaluation_te <- ENMevaluate (dbroteri_te, vars.stack, dbroteri_backdata, method = "jackknife", rasterPreds = T, bin.output=TRUE)
-evaluation_he <- ENMevaluate (dbroteri_he, vars.stack, dbroteri_backdata, method = "jackknife", rasterPreds = T, bin.output=TRUE)
-evaluation_do <- ENMevaluate (dbroteri_do, vars.stack, dbroteri_backdata, method = "jackknife", rasterPreds = T, bin.output=TRUE)
diff --git a/nichobroteri_new.R b/nichobroteri_new.R
deleted file mode 100644
index fbb82c7..0000000
--- a/nichobroteri_new.R
+++ /dev/null
@@ -1,512 +0,0 @@
-library (dismo)
-library (raster)
-# library (rJava)
-# library (rgdal)
-library (rgeos)
-# library (rgbif)
-# library (rjson)
-library (gtools)
-# library (maps)
-# library (ggmap)
-library (fuzzySim)
-# library (ade4)
-# library (pcaMethods)
-library (ecospat)
-library (sp)
-library (GSIF)
-# library (caret)
-# library (RCurl)
-# library (gdalUtils)
-# library (plotKML)
-# library (XML)
-# library (lattice)
-# library (aqp)
-# library (soiltexture)
-library (ggbiplot)
-# library (psych)
-library (spatialEco)
-library (spThin)
-
-
-#dataset de poblaciones con coordenadas
-
-dbroteri <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/Poblaciones_Nicho.csv", sep = ";", fileEncoding = "latin1", colClasses = c("factor", "factor", "numeric", "numeric"))
-ploidy <- dbroteri[,2]
-ploidy <- as.data.frame (ploidy)
-rownames (ploidy) <- c("Monte Clerigo","Sao Bras de Alportel","Zafarraya 1","Zafarraya 2","Orgiva","Laroles","Lliria","Troia","Comporta","Donana (Peladillo)","Albufera de Valencia","Alcublas","Azuebar","Sierra de Espadan","Chiclana","Ronda","Calblanque (Cabezo de la Fuente)","Socovos","Cartagena","San Miguel de Salinas","Penon de Ifach","Valverde del Camino","Moguer","Hinojos","Donana (Acebron)","Donana (Puntal)","Huertos del Batan","Isla Cristina (Las Palmeritas)")
-coordinates (dbroteri) <- ~long + lat
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string (dbroteri) <- crs.geo
-
-
-#carga de variables predictoras y union con mismos limites (chelsa, envirem, altitud, SoilGrids)
-#extraccion de datos de las variables predictoras en las poblaciones
-
-e <- extent (-10,3,35,42)
-
-chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chelsafiles)
-che.c <- crop (chelsa,e)
-
-enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (enviremfiles)
-env.c <- crop (envirem, e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-variables <- stack (che.c, env.c, alt.m)
-names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","aridityIndexThornthwaite","climaticMoistureIndex","continentality","embergerQ","growingDegDays0","growingDegDays5","maxTempColdest","minTempWarmest","monthCountByTemp10","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","thermicityIndex","topoWet","tri","elevation")
-
-soilgrids<-extract.list(dbroteri, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-colnames (soilgrids) <- c("ploidy","AWCh1","AWCh2","AWCh3","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT","TEXMHT")
-soilgrids$ploidy<-factor(soilgrids$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="1","clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="2","silty clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="3","sandy clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="4","clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="5","silty clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="6","sandy clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="7","loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="8","silty loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="9","sandy loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="10","silt")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="11","loamy sand")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="12","sand")
-soilgrids$TEXMHT<-factor(soilgrids$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
-soilgrids<-cbind(soilgrids,apply(soilgrids[,c(2:4)], 1, mean))
-soilgrids<-soilgrids[,-c(1:4)]
-colnames(soilgrids)[7]<-"AWC"
-soilgrids <- soilgrids[,c(7,1,2,3,4,5,6)]
-
-
-presvals <- extract (variables, dbroteri)
-presvals <- cbind (ploidy, presvals, soilgrids) 
-presvals$PHIHOX <- presvals$PHIHOX/10
-
-#calculo del tri a partir de altitud con libreria spatialEco
-
-tri.ext <- tri(alt.m)
-projection(tri.ext) <- crs.geo 
-trivalues<-extract(tri.ext,dbroteri)
-
-
-#sustitucion de los valores tri por los del dataframe con NAs
-presvals <- presvals[,-38]
-presvals <- cbind (presvals, trivalues)
-presvals <- presvals[,c(1:37,46,38:45)]
-colnames(presvals)[38] <- "tri"
-
-
-#analisis para descartar variables muy correlacionadas
-#PCA de puntos de presencia con variables seleccionadas
-
-colvar <- presvals[c(2:45)]
-x <- cor(colvar, method="pearson")
-ecospat.npred (x, th=0.70)
-
-presvals.pca <- presvals[,-c(1,46)]
-presvals.pca <- cbind (presvals.pca,1)
-correlations <- corSelect (presvals.pca, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.70)
-selected <- correlations$selected.var.cols 
-presvals.pca.2 <- presvals.pca[,c(selected)]
-
-ploidy <- ploidy$ploidy
-ploidy <- factor (ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-
-pca <- prcomp(presvals.pca.2, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#PCA con todas las variables
-pca <- prcomp(presvals[,-c(1,46)], scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#===============PREPARATION OF DATAFRAMES=============#
-
-dbroteridata <- as.data.frame(dbroteri)
-presvalsdata <- cbind (dbroteridata[,c(3,4)], presvals)
-presvalsdata <- presvalsdata [,-48]
-presvals2 <- presvalsdata
-presvals2 <- presvals2 [,c(2,1,3:47)]
-coordinates(presvalsdata)<- ~long+ lat
-proj4string(presvalsdata) <- crs.geo 
-
-diploid <- as.data.frame(dbroteri[-c(7:28),])
-diploid <- diploid [,-c(1,2)]
-tetraploid <- as.data.frame(dbroteri[c(7:16),])
-tetraploid <- tetraploid [,-c(1,2)]
-hexaploid <- as.data.frame(dbroteri[c(17:21),])
-hexaploid <- hexaploid [,-c(1,2)]
-dodecaploid <- as.data.frame(dbroteri[c(22:28),])
-dodecaploid <- dodecaploid [,-c(1,2)]  
-
-coordinates(diploid)<- ~long+ lat
-proj4string(diploid) <- crs.geo 
-coordinates(tetraploid)<- ~long+ lat
-proj4string(tetraploid) <- crs.geo 
-coordinates(hexaploid)<- ~long+ lat
-proj4string(hexaploid) <- crs.geo 
-coordinates(dodecaploid)<- ~long+ lat
-proj4string(dodecaploid) <- crs.geo 
-
-#===============PREPARATION OF DATAFRAMES=============#
-
-
-#===============DIPLOID=============#
-
-maskdi <- raster(diploid)
-res(maskdi) <- 0.008333333
-xdi <- circles(diploid, d=50000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-poldi <- gUnaryUnion(xdi@polygons)
-sampdi <- spsample(poldi, 1000, type='random', iter=25)
-extent(maskdi) <- extent(poldi) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellsdi <- cellFromXY(maskdi, sampdi)
-length(cellsdi)
-cellsdi <- unique(cellsdi)
-length(cellsdi)
-xydi <- xyFromCell(maskdi, cellsdi)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xydi <- as.data.frame(xydi)
-diploid <- as.data.frame(diploid)
-colnames(xydi) <- colnames(diploid)
-dibackcoord <- rbind(xydi, diploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(dibackcoord) <- ~long+ lat
-proj4string(dibackcoord) <- crs.geo
-r <- raster(dibackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-dibackcoord_sel <- as.data.frame(gridSample(dibackcoord, r, n=1))
-coordinates(dibackcoord_sel) <- ~long+ lat
-proj4string(dibackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-dibackgroundclim <- extract(variables,dibackcoord_sel)
-dibackgroundsoil <- extract.list(dibackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-dibackgrounddat <- cbind("dibackground",as.data.frame(dibackcoord_sel),dibackgroundclim, dibackgroundsoil)
-dibackgrounddat <- dibackgrounddat[,-42]
-
-coordinates(dibackgrounddat) <- ~long+ lat
-proj4string(dibackgrounddat) <- crs.geo
-trivalues1 <- extract(tri.ext,dibackgrounddat)
-dibackgrounddat <- as.data.frame(dibackgrounddat)
-dibackgrounddat <- dibackgrounddat[,-40]
-dibackgrounddat <- cbind (dibackgrounddat, trivalues1)
-dibackgrounddat <- dibackgrounddat[,c(1:39,50,40:49)]
-colnames(dibackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-dibackgrounddat.c <- na.omit(dibackgrounddat)
-dibackgrounddat.c <- cbind(dibackgrounddat.c,apply(dibackgrounddat.c[,c(42:44)], 1, mean))
-dibackgrounddat.c <- dibackgrounddat.c[,-c(42:44)]
-colnames(dibackgrounddat.c)[48] <- "AWC"
-dibackgrounddat.c <- dibackgrounddat.c[,c(1:41,48,42:47)]
-dibackgrounddat.c <- dibackgrounddat.c[,-48]
-dibackgrounddat.c <- dibackgrounddat.c[,c(2,3,1,4:47)]
-colnames(dibackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(dibackgrounddat.c) <- ~long+ lat
-proj4string(dibackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(dibackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============DIPLOID=============#
-
-
-#===============TETRAPLOID=============#
-
-maskte <- raster(tetraploid)
-res(maskte) <- 0.008333333
-xte <- circles(tetraploid, d=50000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-polte <- gUnaryUnion(xte@polygons)
-sampte <- spsample(polte, 1000, type='random', iter=25)
-extent(maskte) <- extent(polte) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellste <- cellFromXY(maskte, sampte)
-length(cellste)
-cellste <- unique(cellste)
-length(cellste)
-xyte <- xyFromCell(maskte, cellste)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xyte <- as.data.frame(xyte)
-tetraploid <- as.data.frame(tetraploid)
-colnames(xyte) <- colnames(tetraploid)
-tebackcoord <- rbind(xyte, tetraploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(tebackcoord) <- ~long+ lat
-proj4string(tebackcoord) <- crs.geo
-r <- raster(tebackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-tebackcoord_sel <- as.data.frame(gridSample(tebackcoord, r, n=1))
-coordinates(tebackcoord_sel) <- ~long+ lat
-proj4string(tebackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-tebackgroundclim <- extract(variables,tebackcoord_sel)
-tebackgroundsoil <- extract.list(tebackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-tebackgrounddat <- cbind("tebackground",as.data.frame(tebackcoord_sel),tebackgroundclim, tebackgroundsoil)
-tebackgrounddat <- tebackgrounddat[,-42]
-
-coordinates(tebackgrounddat) <- ~long+ lat
-proj4string(tebackgrounddat) <- crs.geo
-trivalues2 <- extract(tri.ext,tebackgrounddat)
-tebackgrounddat <- as.data.frame(tebackgrounddat)
-tebackgrounddat <- tebackgrounddat[,-40]
-tebackgrounddat <- cbind (tebackgrounddat, trivalues2)
-tebackgrounddat <- tebackgrounddat[,c(1:39,50,40:49)]
-colnames(tebackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-tebackgrounddat.c <- na.omit(tebackgrounddat)
-tebackgrounddat.c <- cbind(tebackgrounddat.c,apply(tebackgrounddat.c[,c(42:44)], 1, mean))
-tebackgrounddat.c <- tebackgrounddat.c[,-c(42:44)]
-colnames(tebackgrounddat.c)[48] <- "AWC"
-tebackgrounddat.c <- tebackgrounddat.c[,c(1:41,48,42:47)]
-tebackgrounddat.c <- tebackgrounddat.c[,-48]
-tebackgrounddat.c <- tebackgrounddat.c[,c(2,3,1,4:47)]
-colnames(tebackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(tebackgrounddat.c) <- ~long+ lat
-proj4string(tebackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(tebackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============TETRAPLOID=============#
-
-#===============HEXAPLOID=============#
-
-maskhe <- raster(hexaploid)
-res(maskhe) <- 0.008333333
-xhe <- circles(hexaploid, d=30000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-polhe <- gUnaryUnion(xhe@polygons)
-samphe <- spsample(polhe, 500, type='random', iter=25)
-extent(maskhe) <- extent(polhe) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellshe <- cellFromXY(maskhe, samphe)
-length(cellshe)
-cellshe <- unique(cellshe)
-length(cellshe)
-xyhe <- xyFromCell(maskhe, cellshe)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xyhe <- as.data.frame(xyhe)
-hexaploid <- as.data.frame(hexaploid)
-colnames(xyhe) <- colnames(hexaploid)
-hebackcoord <- rbind(xyhe, hexaploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(hebackcoord) <- ~long+ lat
-proj4string(hebackcoord) <- crs.geo
-r <- raster(hebackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-hebackcoord_sel <- as.data.frame(gridSample(hebackcoord, r, n=1))
-coordinates(hebackcoord_sel) <- ~long+ lat
-proj4string(hebackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-hebackgroundclim <- extract(variables,hebackcoord_sel)
-hebackgroundsoil <- extract.list(hebackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-hebackgrounddat <- cbind("hebackground",as.data.frame(hebackcoord_sel),hebackgroundclim, hebackgroundsoil)
-hebackgrounddat <- hebackgrounddat[,-42]
-
-coordinates(hebackgrounddat) <- ~long+ lat
-proj4string(hebackgrounddat) <- crs.geo
-trivalues3 <- extract(tri.ext,hebackgrounddat)
-hebackgrounddat <- as.data.frame(hebackgrounddat)
-hebackgrounddat <- hebackgrounddat[,-40]
-hebackgrounddat <- cbind (hebackgrounddat, trivalues3)
-hebackgrounddat <- hebackgrounddat[,c(1:39,50,40:49)]
-colnames(hebackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-hebackgrounddat.c <- na.omit(hebackgrounddat)
-hebackgrounddat.c <- cbind(hebackgrounddat.c,apply(hebackgrounddat.c[,c(42:44)], 1, mean))
-hebackgrounddat.c <- hebackgrounddat.c[,-c(42:44)]
-colnames(hebackgrounddat.c)[48] <- "AWC"
-hebackgrounddat.c <- hebackgrounddat.c[,c(1:41,48,42:47)]
-hebackgrounddat.c <- hebackgrounddat.c[,-48]
-hebackgrounddat.c <- hebackgrounddat.c[,c(2,3,1,4:47)]
-colnames(hebackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(hebackgrounddat.c) <- ~long+ lat
-proj4string(hebackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(hebackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============HEXAPLOID=============#
-
-
-#===============DODECAPLOID=============#
-
-maskdo <- raster(dodecaploid)
-res(maskdo) <- 0.008333333
-xdo <- circles(dodecaploid, d=25000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-poldo <- gUnaryUnion(xdo@polygons)
-sampdo <- spsample(poldo, 250, type='random', iter=25)
-extent(maskdo) <- extent(poldo) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellsdo <- cellFromXY(maskdo, sampdo)
-length(cellsdo)
-cellsdo <- unique(cellsdo)
-length(cellsdo)
-xydo <- xyFromCell(maskdo, cellsdo)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xydo <- as.data.frame(xydo)
-dodecaploid <- as.data.frame(dodecaploid)
-colnames(xydo) <- colnames(dodecaploid)
-dobackcoord <- rbind(xydo, dodecaploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(dobackcoord) <- ~long+ lat
-proj4string(dobackcoord) <- crs.geo
-r <- raster(dobackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-dobackcoord_sel <- as.data.frame(gridSample(dobackcoord, r, n=1))
-coordinates(dobackcoord_sel) <- ~long+ lat
-proj4string(dobackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-dobackgroundclim <- extract(variables,dobackcoord_sel)
-dobackgroundsoil <- extract.list(dobackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-dobackgrounddat <- cbind("dobackground",as.data.frame(dobackcoord_sel),dobackgroundclim, dobackgroundsoil)
-dobackgrounddat <- dobackgrounddat[,-42]
-
-coordinates(dobackgrounddat) <- ~long+ lat
-proj4string(dobackgrounddat) <- crs.geo
-trivalues4 <- extract(tri.ext,dobackgrounddat)
-dobackgrounddat <- as.data.frame(dobackgrounddat)
-dobackgrounddat <- dobackgrounddat[,-40]
-dobackgrounddat <- cbind (dobackgrounddat, trivalues4)
-dobackgrounddat <- dobackgrounddat[,c(1:39,50,40:49)]
-colnames(dobackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-dobackgrounddat.c <- na.omit(dobackgrounddat)
-dobackgrounddat.c <- cbind(dobackgrounddat.c,apply(dobackgrounddat.c[,c(42:44)], 1, mean))
-dobackgrounddat.c <- dobackgrounddat.c[,-c(42:44)]
-colnames(dobackgrounddat.c)[48] <- "AWC"
-dobackgrounddat.c <- dobackgrounddat.c[,c(1:41,48,42:47)]
-dobackgrounddat.c <- dobackgrounddat.c[,-48]
-dobackgrounddat.c <- dobackgrounddat.c[,c(2,3,1,4:47)]
-colnames(dobackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(dobackgrounddat.c) <- ~long+ lat
-proj4string(dobackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(dobackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============DODECAPLOID=============#
-
-#===============GENERAL BACKGROUND=============#
-
-backgrounddat.c <- rbind (presvals2, as.data.frame(dibackgrounddat.c), as.data.frame(tebackgrounddat.c), as.data.frame(hebackgrounddat.c), as.data.frame(dobackgrounddat.c))
-backgrounddat.c$ploidy <- "background"
-coordinates(backgrounddat.c) <- ~long+ lat
-proj4string(backgrounddat.c) <- crs.geo
-
-plot(gmap(e, type = "satellite"))
-points(Mercator(backgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============GENERAL BACKGROUND=============#
-
-#===============PCA BACKGROUND=============#
-
-todo <- rbind (presvals2, as.data.frame(backgrounddat.c), as.data.frame(dibackgrounddat.c), as.data.frame(tebackgrounddat.c), as.data.frame(hebackgrounddat.c), as.data.frame(dobackgrounddat.c))
-
-todo.pca <- todo[,-c(1:3)]
-todo.pca <- cbind (todo.pca,1)
-correlations2 <- corSelect (todo.pca, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.75)
-selected2 <- correlations2$selected.var.cols 
-todo.pca.2 <- todo.pca[,c(selected2)]
-
-todoploidy <- todo$ploidy
-w<-c(rep(0,nrow(presvals2)),rep(1,nrow(as.data.frame(backgrounddat.c))),rep(0,nrow(as.data.frame(dibackgrounddat.c))), rep(0,nrow(as.data.frame(tebackgrounddat.c))), rep(0,nrow(as.data.frame(hebackgrounddat.c))), rep(0,nrow(as.data.frame(dobackgrounddat.c))))
-
-pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
-gcol = c("blue", "red", "green", "yellow", "orange", "purple", "cyan", "black", "grey")
-s.label(pcaback$li, clabel = 0.1)
-scatter(pcaback, clab.row = 0, posieig = "none", cex=0.1)
-s.class(pcaback$li, todo[,3], col = gcol, add.plot = TRUE, cstar = 0, clabel = 0, cellipse = 1.5, pch = 16)
-
-
-#===============ECOSPAT=============#
-
-row.di<-which(todo[,3] == "2x")
-row.te<-which(todo[,3] == "4x")
-row.he<-which(todo[,3] == "6x")
-row.do<-which(todo[,3] == "12x")
-row.back<-which(todo[,3] == "background") 
-row.bacdi<-which(todo[,3] == "dibackground") 
-row.bacte<-which(todo[,3] == "tebackground")
-row.bache<-which(todo[,3] == "hebackground")
-row.bacdo<-which(todo[,3] == "dobackground")
-
-scores.clim<- pcaback$li[row.back,] 
-scores.di<- pcaback$li[row.di,]		
-scores.te<- pcaback$li[row.te,]	
-scores.he<- pcaback$li[row.he,]	
-scores.do<- pcaback$li[row.do,]	
-scores.bacdi<- pcaback$li[row.bacdi,]					
-scores.bacte<- pcaback$li[row.bacte,]	
-scores.bache<- pcaback$li[row.bache,]	
-scores.bacdo<- pcaback$li[row.bacdo,]	
-
-zdi<- ecospat.grid.clim.dyn (scores.clim, scores.bacdi, scores.di, R=10)
-zte<- ecospat.grid.clim.dyn (scores.clim, scores.bacte, scores.te, R=10)
-zhe<- ecospat.grid.clim.dyn (scores.clim, scores.bache, scores.he, R=10)
-zdo<- ecospat.grid.clim.dyn (scores.clim, scores.bacdo, scores.do, R=10)
-
-
-#ecospat (tests)
-equivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 1000, alternative = "lower")
-equivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 1000, alternative = "lower")
-equivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 1000, alternative = "lower")
-equivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 1000, alternative = "lower")
-equivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 1000, alternative = "lower")
-equivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 1000, alternative = "lower")
-
-overlap.test.dite<-ecospat.niche.overlap (zdi, zte, cor=FALSE)
-overlap.test.dihe<-ecospat.niche.overlap (zdi, zhe, cor=FALSE)
-overlap.test.dido<-ecospat.niche.overlap (zdi, zdo, cor=FALSE)
-overlap.test.tehe<-ecospat.niche.overlap (zte, zhe, cor=FALSE)
-overlap.test.tedo<-ecospat.niche.overlap (zte, zdo, cor=FALSE)
-overlap.test.hedo<-ecospat.niche.overlap (zhe, zdo, cor=FALSE)
-
-ecospat.plot.niche.dyn (zdi, zdo, quant = 0.75)
-
-
-similarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 1000, alternative = "lower")
-similarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 1000, alternative = "lower")
-similarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 1000, alternative = "lower")
-similarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 1000, alternative = "lower")
-similarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 1000, alternative = "lower")
-similarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 1000, alternative = "lower")
-similarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 1000, alternative = "lower")
-similarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 1000, alternative = "lower")
-similarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 1000, alternative = "lower")
-similarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 1000, alternative = "lower")
-similarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 1000, alternative = "lower")
-similarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 1000, alternative = "lower")
\ No newline at end of file
diff --git a/nichobroteri_newfran.R b/nichobroteri_newfran.R
deleted file mode 100644
index e062535..0000000
--- a/nichobroteri_newfran.R
+++ /dev/null
@@ -1,547 +0,0 @@
- # Diferencias de nicho usando el mismo background (i.e no corrigiendo por l os background de las distintas especies)
-
-
-library (dismo)
-library (raster)
-# library (rJava)
-# library (rgdal)
-library (rgeos)
-# library (rgbif)
-# library (rjson)
-library (gtools)
-# library (maps)
-# library (ggmap)
-library (fuzzySim)
-# library (ade4)
-# library (pcaMethods)
-library (ecospat)
-library (sp)
-library (GSIF)
-# library (caret)
-# library (RCurl)
-# library (gdalUtils)
-# library (plotKML)
-# library (XML)
-# library (lattice)
-# library (aqp)
-# library (soiltexture)
-library (ggbiplot)
-# library (psych)
-library (spatialEco)
-library (spThin)
-
-
-#dataset de poblaciones con coordenadas
-
-dbroteri <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/Poblaciones_Nicho.csv", sep = ";", fileEncoding = "latin1", colClasses = c("factor", "factor", "numeric", "numeric"))
-ploidy <- dbroteri[,2]
-ploidy <- as.data.frame (ploidy)
-rownames (ploidy) <- c("Monte Clerigo","Sao Bras de Alportel","Zafarraya 1","Zafarraya 2","Orgiva","Laroles","Lliria","Troia","Comporta","Donana (Peladillo)","Albufera de Valencia","Alcublas","Azuebar","Sierra de Espadan","Chiclana","Ronda","Calblanque (Cabezo de la Fuente)","Socovos","Cartagena","San Miguel de Salinas","Penon de Ifach","Valverde del Camino","Moguer","Hinojos","Donana (Acebron)","Donana (Puntal)","Huertos del Batan","Isla Cristina (Las Palmeritas)")
-coordinates (dbroteri) <- ~long + lat
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string (dbroteri) <- crs.geo
-
-
-#carga de variables predictoras y union con mismos limites (chelsa, envirem, altitud, SoilGrids)
-#extraccion de datos de las variables predictoras en las poblaciones
-
-e <- extent (-10,3,35,42)
-
-chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chelsafiles)
-che.c <- crop (chelsa,e)
-
-enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (enviremfiles)
-env.c <- crop (envirem, e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-variables <- stack (che.c, env.c, alt.m)
-names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","aridityIndexThornthwaite","climaticMoistureIndex","continentality","embergerQ","growingDegDays0","growingDegDays5","maxTempColdest","minTempWarmest","monthCountByTemp10","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","thermicityIndex","topoWet","tri","elevation")
-
-soilgrids<-extract.list(dbroteri, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-colnames (soilgrids) <- c("ploidy","AWCh1","AWCh2","AWCh3","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT","TEXMHT")
-soilgrids$ploidy<-factor(soilgrids$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="1","clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="2","silty clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="3","sandy clay")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="4","clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="5","silty clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="6","sandy clay loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="7","loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="8","silty loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="9","sandy loam")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="10","silt")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="11","loamy sand")
-soilgrids$TEXMHT<-replace(soilgrids$TEXMHT,soilgrids$TEXMHT=="12","sand")
-soilgrids$TEXMHT<-factor(soilgrids$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
-soilgrids<-cbind(soilgrids,apply(soilgrids[,c(2:4)], 1, mean))
-soilgrids<-soilgrids[,-c(1:4)]
-colnames(soilgrids)[7]<-"AWC"
-soilgrids <- soilgrids[,c(7,1,2,3,4,5,6)]
-
-
-presvals <- extract (variables, dbroteri)
-presvals <- cbind (ploidy, presvals, soilgrids) 
-presvals$PHIHOX <- presvals$PHIHOX/10
-
-#calculo del tri a partir de altitud con libreria spatialEco
-
-tri.ext <- tri(alt.m)
-projection(tri.ext) <- crs.geo 
-trivalues<-extract(tri.ext,dbroteri)
-
-
-#sustitucion de los valores tri por los del dataframe con NAs
-presvals <- presvals[,-38]
-presvals <- cbind (presvals, trivalues)
-presvals <- presvals[,c(1:37,46,38:45)]
-colnames(presvals)[38] <- "tri"
-
-
-#analisis para descartar variables muy correlacionadas
-#PCA de puntos de presencia con variables seleccionadas
-
-colvar <- presvals[c(2:45)]
-x <- cor(colvar, method="pearson")
-ecospat.npred (x, th=0.70)
-
-presvals.pca <- presvals[,-c(1,46)]
-presvals.pca <- cbind (presvals.pca,1)
-correlations <- corSelect (presvals.pca, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.70)
-selected <- correlations$selected.var.cols 
-presvals.pca.2 <- presvals.pca[,c(selected)]
-
-ploidy <- ploidy$ploidy
-ploidy <- factor (ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-
-pca <- prcomp(presvals.pca.2, scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#PCA con todas las variables
-pca <- prcomp(presvals[,-c(1,46)], scale. = TRUE, retx = T)
-ggbiplot(pca, obs.scale = 1,var.scale = 1,
-         groups = ploidy, ellipse = TRUE, circle = FALSE, alpha =  1) +
-  scale_color_discrete(name = '') +
-  geom_point(aes(colour=ploidy), size = 3) +
-  theme(legend.direction = 'vertical', legend.position = 'right')
-
-
-#===============PREPARATION OF DATAFRAMES=============#
-
-dbroteridata <- as.data.frame(dbroteri)
-presvalsdata <- cbind (dbroteridata[,c(3,4)], presvals)
-presvalsdata <- presvalsdata [,-48]
-presvals2 <- presvalsdata
-presvals2 <- presvals2 [,c(2,1,3:47)]
-coordinates(presvalsdata)<- ~long+ lat
-proj4string(presvalsdata) <- crs.geo 
-
-diploid <- as.data.frame(dbroteri[-c(7:28),])
-diploid <- diploid [,-c(1,2)]
-tetraploid <- as.data.frame(dbroteri[c(7:16),])
-tetraploid <- tetraploid [,-c(1,2)]
-hexaploid <- as.data.frame(dbroteri[c(17:21),])
-hexaploid <- hexaploid [,-c(1,2)]
-dodecaploid <- as.data.frame(dbroteri[c(22:28),])
-dodecaploid <- dodecaploid [,-c(1,2)]  
-
-coordinates(diploid)<- ~long+ lat
-proj4string(diploid) <- crs.geo 
-coordinates(tetraploid)<- ~long+ lat
-proj4string(tetraploid) <- crs.geo 
-coordinates(hexaploid)<- ~long+ lat
-proj4string(hexaploid) <- crs.geo 
-coordinates(dodecaploid)<- ~long+ lat
-proj4string(dodecaploid) <- crs.geo 
-
-#===============PREPARATION OF DATAFRAMES=============#
-
-
-#===============DIPLOID=============#
-
-maskdi <- raster(diploid)
-res(maskdi) <- 0.008333333
-xdi <- circles(diploid, d=50000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-poldi <- gUnaryUnion(xdi@polygons)
-sampdi <- spsample(poldi, 1000, type='random', iter=25)
-extent(maskdi) <- extent(poldi) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellsdi <- cellFromXY(maskdi, sampdi)
-length(cellsdi)
-cellsdi <- unique(cellsdi)
-length(cellsdi)
-xydi <- xyFromCell(maskdi, cellsdi)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xydi <- as.data.frame(xydi)
-diploid <- as.data.frame(diploid)
-colnames(xydi) <- colnames(diploid)
-dibackcoord <- rbind(xydi, diploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(dibackcoord) <- ~long+ lat
-proj4string(dibackcoord) <- crs.geo
-r <- raster(dibackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-dibackcoord_sel <- as.data.frame(gridSample(dibackcoord, r, n=1))
-coordinates(dibackcoord_sel) <- ~long+ lat
-proj4string(dibackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-dibackgroundclim <- extract(variables,dibackcoord_sel)
-dibackgroundsoil <- extract.list(dibackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-dibackgrounddat <- cbind("dibackground",as.data.frame(dibackcoord_sel),dibackgroundclim, dibackgroundsoil)
-dibackgrounddat <- dibackgrounddat[,-42]
-
-coordinates(dibackgrounddat) <- ~long+ lat
-proj4string(dibackgrounddat) <- crs.geo
-trivalues1 <- extract(tri.ext,dibackgrounddat)
-dibackgrounddat <- as.data.frame(dibackgrounddat)
-dibackgrounddat <- dibackgrounddat[,-40]
-dibackgrounddat <- cbind (dibackgrounddat, trivalues1)
-dibackgrounddat <- dibackgrounddat[,c(1:39,50,40:49)]
-colnames(dibackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-dibackgrounddat.c <- na.omit(dibackgrounddat)
-dibackgrounddat.c <- cbind(dibackgrounddat.c,apply(dibackgrounddat.c[,c(42:44)], 1, mean))
-dibackgrounddat.c <- dibackgrounddat.c[,-c(42:44)]
-colnames(dibackgrounddat.c)[48] <- "AWC"
-dibackgrounddat.c <- dibackgrounddat.c[,c(1:41,48,42:47)]
-dibackgrounddat.c <- dibackgrounddat.c[,-48]
-dibackgrounddat.c <- dibackgrounddat.c[,c(2,3,1,4:47)]
-colnames(dibackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(dibackgrounddat.c) <- ~long+ lat
-proj4string(dibackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(dibackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============DIPLOID=============#
-
-
-#===============TETRAPLOID=============#
-
-maskte <- raster(tetraploid)
-res(maskte) <- 0.008333333
-xte <- circles(tetraploid, d=50000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-polte <- gUnaryUnion(xte@polygons)
-sampte <- spsample(polte, 1000, type='random', iter=25)
-extent(maskte) <- extent(polte) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellste <- cellFromXY(maskte, sampte)
-length(cellste)
-cellste <- unique(cellste)
-length(cellste)
-xyte <- xyFromCell(maskte, cellste)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xyte <- as.data.frame(xyte)
-tetraploid <- as.data.frame(tetraploid)
-colnames(xyte) <- colnames(tetraploid)
-tebackcoord <- rbind(xyte, tetraploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(tebackcoord) <- ~long+ lat
-proj4string(tebackcoord) <- crs.geo
-r <- raster(tebackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-tebackcoord_sel <- as.data.frame(gridSample(tebackcoord, r, n=1))
-coordinates(tebackcoord_sel) <- ~long+ lat
-proj4string(tebackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-tebackgroundclim <- extract(variables,tebackcoord_sel)
-tebackgroundsoil <- extract.list(tebackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-tebackgrounddat <- cbind("tebackground",as.data.frame(tebackcoord_sel),tebackgroundclim, tebackgroundsoil)
-tebackgrounddat <- tebackgrounddat[,-42]
-
-coordinates(tebackgrounddat) <- ~long+ lat
-proj4string(tebackgrounddat) <- crs.geo
-trivalues2 <- extract(tri.ext,tebackgrounddat)
-tebackgrounddat <- as.data.frame(tebackgrounddat)
-tebackgrounddat <- tebackgrounddat[,-40]
-tebackgrounddat <- cbind (tebackgrounddat, trivalues2)
-tebackgrounddat <- tebackgrounddat[,c(1:39,50,40:49)]
-colnames(tebackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-tebackgrounddat.c <- na.omit(tebackgrounddat)
-tebackgrounddat.c <- cbind(tebackgrounddat.c,apply(tebackgrounddat.c[,c(42:44)], 1, mean))
-tebackgrounddat.c <- tebackgrounddat.c[,-c(42:44)]
-colnames(tebackgrounddat.c)[48] <- "AWC"
-tebackgrounddat.c <- tebackgrounddat.c[,c(1:41,48,42:47)]
-tebackgrounddat.c <- tebackgrounddat.c[,-48]
-tebackgrounddat.c <- tebackgrounddat.c[,c(2,3,1,4:47)]
-colnames(tebackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(tebackgrounddat.c) <- ~long+ lat
-proj4string(tebackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(tebackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============TETRAPLOID=============#
-
-#===============HEXAPLOID=============#
-
-maskhe <- raster(hexaploid)
-res(maskhe) <- 0.008333333
-xhe <- circles(hexaploid, d=30000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-polhe <- gUnaryUnion(xhe@polygons)
-samphe <- spsample(polhe, 500, type='random', iter=25)
-extent(maskhe) <- extent(polhe) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellshe <- cellFromXY(maskhe, samphe)
-length(cellshe)
-cellshe <- unique(cellshe)
-length(cellshe)
-xyhe <- xyFromCell(maskhe, cellshe)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xyhe <- as.data.frame(xyhe)
-hexaploid <- as.data.frame(hexaploid)
-colnames(xyhe) <- colnames(hexaploid)
-hebackcoord <- rbind(xyhe, hexaploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(hebackcoord) <- ~long+ lat
-proj4string(hebackcoord) <- crs.geo
-r <- raster(hebackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-hebackcoord_sel <- as.data.frame(gridSample(hebackcoord, r, n=1))
-coordinates(hebackcoord_sel) <- ~long+ lat
-proj4string(hebackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-hebackgroundclim <- extract(variables,hebackcoord_sel)
-hebackgroundsoil <- extract.list(hebackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-hebackgrounddat <- cbind("hebackground",as.data.frame(hebackcoord_sel),hebackgroundclim, hebackgroundsoil)
-hebackgrounddat <- hebackgrounddat[,-42]
-
-coordinates(hebackgrounddat) <- ~long+ lat
-proj4string(hebackgrounddat) <- crs.geo
-trivalues3 <- extract(tri.ext,hebackgrounddat)
-hebackgrounddat <- as.data.frame(hebackgrounddat)
-hebackgrounddat <- hebackgrounddat[,-40]
-hebackgrounddat <- cbind (hebackgrounddat, trivalues3)
-hebackgrounddat <- hebackgrounddat[,c(1:39,50,40:49)]
-colnames(hebackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-hebackgrounddat.c <- na.omit(hebackgrounddat)
-hebackgrounddat.c <- cbind(hebackgrounddat.c,apply(hebackgrounddat.c[,c(42:44)], 1, mean))
-hebackgrounddat.c <- hebackgrounddat.c[,-c(42:44)]
-colnames(hebackgrounddat.c)[48] <- "AWC"
-hebackgrounddat.c <- hebackgrounddat.c[,c(1:41,48,42:47)]
-hebackgrounddat.c <- hebackgrounddat.c[,-48]
-hebackgrounddat.c <- hebackgrounddat.c[,c(2,3,1,4:47)]
-colnames(hebackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(hebackgrounddat.c) <- ~long+ lat
-proj4string(hebackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(hebackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============HEXAPLOID=============#
-
-
-#===============DODECAPLOID=============#
-
-maskdo <- raster(dodecaploid)
-res(maskdo) <- 0.008333333
-xdo <- circles(dodecaploid, d=25000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-poldo <- gUnaryUnion(xdo@polygons)
-sampdo <- spsample(poldo, 250, type='random', iter=25)
-extent(maskdo) <- extent(poldo) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cellsdo <- cellFromXY(maskdo, sampdo)
-length(cellsdo)
-cellsdo <- unique(cellsdo)
-length(cellsdo)
-xydo <- xyFromCell(maskdo, cellsdo)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xydo <- as.data.frame(xydo)
-dodecaploid <- as.data.frame(dodecaploid)
-colnames(xydo) <- colnames(dodecaploid)
-dobackcoord <- rbind(xydo, dodecaploid)
-
-# Quitamos los puntos en el mismo km^2
-coordinates(dobackcoord) <- ~long+ lat
-proj4string(dobackcoord) <- crs.geo
-r <- raster(dobackcoord)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-dobackcoord_sel <- as.data.frame(gridSample(dobackcoord, r, n=1))
-coordinates(dobackcoord_sel) <- ~long+ lat
-proj4string(dobackcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-dobackgroundclim <- extract(variables,dobackcoord_sel)
-dobackgroundsoil <- extract.list(dobackcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-dobackgrounddat <- cbind("dobackground",as.data.frame(dobackcoord_sel),dobackgroundclim, dobackgroundsoil)
-dobackgrounddat <- dobackgrounddat[,-42]
-
-coordinates(dobackgrounddat) <- ~long+ lat
-proj4string(dobackgrounddat) <- crs.geo
-trivalues4 <- extract(tri.ext,dobackgrounddat)
-dobackgrounddat <- as.data.frame(dobackgrounddat)
-dobackgrounddat <- dobackgrounddat[,-40]
-dobackgrounddat <- cbind (dobackgrounddat, trivalues4)
-dobackgrounddat <- dobackgrounddat[,c(1:39,50,40:49)]
-colnames(dobackgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-dobackgrounddat.c <- na.omit(dobackgrounddat)
-dobackgrounddat.c <- cbind(dobackgrounddat.c,apply(dobackgrounddat.c[,c(42:44)], 1, mean))
-dobackgrounddat.c <- dobackgrounddat.c[,-c(42:44)]
-colnames(dobackgrounddat.c)[48] <- "AWC"
-dobackgrounddat.c <- dobackgrounddat.c[,c(1:41,48,42:47)]
-dobackgrounddat.c <- dobackgrounddat.c[,-48]
-dobackgrounddat.c <- dobackgrounddat.c[,c(2,3,1,4:47)]
-colnames(dobackgrounddat.c) <- colnames(presvals2)
-
-# Representacion de los puntos en el mapa alrededor de los de presencia
-coordinates(dobackgrounddat.c) <- ~long+ lat
-proj4string(dobackgrounddat.c) <- crs.geo
-plot(gmap(e, type = "satellite"))
-points(Mercator(dobackgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============DODECAPLOID=============#
-
-#===============GENERAL BACKGROUND=============#
-
-backgrounddat.c <- rbind (presvals2, as.data.frame(dibackgrounddat.c), as.data.frame(tebackgrounddat.c), as.data.frame(hebackgrounddat.c), as.data.frame(dobackgrounddat.c))
-backgrounddat.c$ploidy <- "background"
-coordinates(backgrounddat.c) <- ~long+ lat
-proj4string(backgrounddat.c) <- crs.geo
-
-plot(gmap(e, type = "satellite"))
-points(Mercator(backgrounddat.c), col = 'orange', pch=20, cex=1)
-points(Mercator(presvalsdata), col=presvalsdata$ploidy, pch=20, cex=1)
-
-#===============GENERAL BACKGROUND=============#
-
-#===============PCA BACKGROUND=============#
-
-todo <- rbind (presvals2, as.data.frame(backgrounddat.c))
-todo.pca <- todo[,-c(1:3)]
-source("vif_fun.R")
-selected2<-vif_func(todo.pca)
-todo.pca.2 <- todo.pca[,c(selected2)]
-
-todoploidy <- todo$ploidy
-w<-c(rep(0,nrow(presvals2)),rep(1,nrow(as.data.frame(backgrounddat.c))))
-
-pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
-gcol = c("blue", "red", "green", "yellow", "grey")
-s.label(pcaback$li, clabel = 0.1)
-scatter(pcaback, clab.row = 0, posieig = "none", cex=0.1)
-s.class(pcaback$li, todo[,3], col = gcol, add.plot = TRUE, cstar = 0, clabel = 0, cellipse = 1.5, pch = 16)
-
-
-#===============ECOSPAT=============#
-
-row.di<-which(todo[,3] == "2x")
-row.te<-which(todo[,3] == "4x")
-row.he<-which(todo[,3] == "6x")
-row.do<-which(todo[,3] == "12x")
-row.back<-which(todo[,3] == "background") 
-
-
-scores.clim<- pcaback$li 
-scores.di<- pcaback$li[row.di,]		
-scores.te<- pcaback$li[row.te,]	
-scores.he<- pcaback$li[row.he,]	
-scores.do<- pcaback$li[row.do,]	
-
-
-zdi<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.di, R=100)
-zte<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.te, R=100)
-zhe<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.he, R=100)
-zdo<- ecospat.grid.clim.dyn (scores.clim, scores.clim, scores.do, R=100)
-
-ecospat.plot.niche(zdi)
-ecospat.plot.niche(zte)
-ecospat.plot.niche(zhe)
-ecospat.plot.niche(zdo)
-
-#ecospat (tests)
-equivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 100, alternative = "lower", ncores = detectCores() - 1)
-equivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 100, alternative = "lower", ncores = detectCores() - 1)
-equivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 100, alternative = "lower", ncores = detectCores() - 1)
-equivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 100, alternative = "lower", ncores = detectCores() - 1)
-equivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 100, alternative = "lower", ncores = detectCores() - 1)
-equivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 100, alternative = "lower", ncores = detectCores() - 1)
-
-gequivalency.test.dite<-ecospat.niche.equivalency.test (zdi, zte, 100, alternative = "greater", ncores = detectCores() - 1)
-gequivalency.test.dihe<-ecospat.niche.equivalency.test (zdi, zhe, 100, alternative = "greater", ncores = detectCores() - 1)
-gequivalency.test.dido<-ecospat.niche.equivalency.test (zdi, zdo, 100, alternative = "greater", ncores = detectCores() - 1)
-gequivalency.test.tehe<-ecospat.niche.equivalency.test (zte, zhe, 100, alternative = "greater", ncores = detectCores() - 1)
-gequivalency.test.tedo<-ecospat.niche.equivalency.test (zte, zdo, 100, alternative = "greater", ncores = detectCores() - 1)
-gequivalency.test.hedo<-ecospat.niche.equivalency.test (zhe, zdo, 100, alternative = "greater", ncores = detectCores() - 1)
-
-overlap.test.dite<-ecospat.niche.overlap (zdi, zte, cor=FALSE)
-overlap.test.dihe<-ecospat.niche.overlap (zdi, zhe, cor=FALSE)
-overlap.test.dido<-ecospat.niche.overlap (zdi, zdo, cor=FALSE)
-overlap.test.tehe<-ecospat.niche.overlap (zte, zhe, cor=FALSE)
-overlap.test.tedo<-ecospat.niche.overlap (zte, zdo, cor=FALSE)
-overlap.test.hedo<-ecospat.niche.overlap (zhe, zdo, cor=FALSE)
-
-ecospat.plot.niche.dyn (zte, zhe, quant = 0.75)
-
-
-similarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 100, alternative = "lower", ncores = detectCores() - 1)
-similarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 100, alternative = "lower", ncores = detectCores() - 1)
-similarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 100, alternative = "lower", ncores = detectCores() - 1)
-similarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 100, alternative = "lower", ncores = detectCores() - 1)
-similarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 100, alternative = "lower", ncores = detectCores() - 1)
-similarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 1000, alternative = "lower", ncores = detectCores() - 1)
-similarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 1000, alternative = "lower", ncores = detectCores() - 1)
-
-gsimilarity.testdite<-ecospat.niche.similarity.test (zdi, zte, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testtedi<-ecospat.niche.similarity.test (zte, zdi, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testdihe<-ecospat.niche.similarity.test (zdi, zhe, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testhedi<-ecospat.niche.similarity.test (zhe, zdi, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testdido<-ecospat.niche.similarity.test (zdi, zdo, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testdodi<-ecospat.niche.similarity.test (zdo, zdi, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testtehe<-ecospat.niche.similarity.test (zte, zhe, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testhete<-ecospat.niche.similarity.test (zhe, zte, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testtedo<-ecospat.niche.similarity.test (zte, zdo, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testdote<-ecospat.niche.similarity.test (zdo, zte, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testhedo<-ecospat.niche.similarity.test (zhe, zdo, 1000, alternative = "greater", ncores = detectCores() - 1)
-gsimilarity.testdohe<-ecospat.niche.similarity.test (zdo, zhe, 1000, alternative = "greater", ncores = detectCores() - 1)
-
-
-# Tabla
-
-overlap<-rbind(overlap.test.dite$D,overlap.test.dihe$D, overlap.test.dido$D,overlap.test.tehe$D, overlap.test.tedo$D, overlap.test.hedo$D)
-
-lsimilarityab<-rbind(similarity.testdite$p.D, similarity.testdihe$p.D, similarity.testdido$p.D, similarity.testtehe$p.D, similarity.testtedo$p.D, similarity.testhedo$p.D)
-lsimilarityba<-rbind(similarity.testtedi$p.D, similarity.testhedi$p.D, similarity.testdodi$p.D, similarity.testhete$p.D, similarity.testdote$p.D, similarity.testdohe$p.D)
-
-gsimilarityab<-rbind(gsimilarity.testdite$p.D, gsimilarity.testdihe$p.D, gsimilarity.testdido$p.D, gsimilarity.testtehe$p.D, gsimilarity.testtedo$p.D, gsimilarity.testhedo$p.D)
-gsimilarityba<-rbind(gsimilarity.testtedi$p.D, gsimilarity.testhedi$p.D, gsimilarity.testdodi$p.D, gsimilarity.testhete$p.D, gsimilarity.testdote$p.D, gsimilarity.testdohe$p.D)
-
-equivalency<-rbind(equivalency.test.dite$p.D, equivalency.test.dihe$p.D, equivalency.test.dido$p.D, equivalency.test.tehe$p.D, equivalency.test.tedo$p.D, equivalency.test.hedo$p.D)
-gequivalency<-rbind(gequivalency.test.dite$p.D, gequivalency.test.dihe$p.D, gequivalency.test.dido$p.D, gequivalency.test.tehe$p.D, gequivalency.test.tedo$p.D, gequivalency.test.hedo$p.D)
-
-tablaresul<-data.frame(overlap,lsimilarityab,gsimilarityab, lsimilarityba,gsimilarityba,equivalency, gequivalency)
\ No newline at end of file
diff --git a/phy_prepdataset.R b/phy_prepdataset.R
deleted file mode 100644
index c3e654b..0000000
--- a/phy_prepdataset.R
+++ /dev/null
@@ -1,163 +0,0 @@
-library (sp)
-library (GSIF)
-library (gtools)
-library (raster)
-library (spatialEco)
-library (dismo)
-library (ade4)
-
-
-dbroteri_arbol <- read.delim2 (file="D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/Poblaciones_Nicho_arbol.csv", sep = ";", fileEncoding = "latin1", colClasses = c("factor", "factor", "numeric", "numeric"))
-ploidy_arbol <- dbroteri_arbol[,2]
-ploidy_arbol <- as.data.frame (ploidy_arbol)
-coordinates(dbroteri_arbol) <- ~long+ lat
-crs.geo <- CRS ("+proj=longlat +ellps=WGS84 +datum=WGS84")
-proj4string(dbroteri_arbol) <- crs.geo
-
-soilgrids_arbol<-extract.list(dbroteri_arbol, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-colnames (soilgrids_arbol) <- c("ploidy","AWCh1","AWCh2","AWCh3","BLDFIE","CECSOL","ORCDRC","PHIHOX","SNDPPT","TEXMHT")
-soilgrids_arbol$ploidy<-factor(soilgrids_arbol$ploidy, levels = c("2x", "4x", "6x", "12x"), ordered = TRUE)
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="1","clay")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="2","silty clay")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="3","sandy clay")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="4","clay loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="5","silty clay loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="6","sandy clay loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="7","loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="8","silty loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="9","sandy loam")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="10","silt")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="11","loamy sand")
-soilgrids_arbol$TEXMHT<-replace(soilgrids_arbol$TEXMHT,soilgrids_arbol$TEXMHT=="12","sand")
-soilgrids_arbol$TEXMHT<-factor(soilgrids_arbol$TEXMHT,levels = c("clay", "silty clay", "sandy clay", "clay loam","silty clay loam","sandy clay loam","loam","silty loam","sandy loam","silt","loamy sand","sand"))
-soilgrids_arbol<-cbind(soilgrids_arbol,apply(soilgrids_arbol[,c(2:4)], 1, mean))
-soilgrids_arbol<-soilgrids_arbol[,-c(1:4)]
-colnames(soilgrids_arbol)[7]<-"AWC"
-soilgrids_arbol <- soilgrids_arbol[,c(7,1,2,3,4,5,6)]
-
-#carga de variables predictoras y union con mismos limites (chelsa, envirem, altitud, SoilGrids)
-#extraccion de datos de las variables predictoras en las poblaciones
-
-e <- extent (-10,3,35,42)
-
-chelsafiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chelsafiles)
-che.c <- crop (chelsa,e)
-
-enviremfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (enviremfiles)
-env.c <- crop (envirem, e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-variables <- stack (che.c, env.c, alt.m)
-names (variables) <- c("bio1","bio2","bio3","bio4","bio5","bio6","bio7","bio8","bio9","bio10","bio11","bio12","bio13","bio14","bio15","bio16","bio17","bio18","bio19","annualPET","aridityIndexThornthwaite","climaticMoistureIndex","continentality","embergerQ","growingDegDays0","growingDegDays5","maxTempColdest","minTempWarmest","monthCountByTemp10","PETColdestQuarter","PETDriestQuarter","PETseasonality","PETWarmestQuarter","PETWettestQuarter","thermicityIndex","topoWet","tri","elevation")
-
-presvals_arbol <- extract (variables, dbroteri_arbol)
-presvals_arbol <- cbind (ploidy_arbol, presvals_arbol, soilgrids_arbol) 
-presvals_arbol$PHIHOX <- presvals_arbol$PHIHOX/10
-
-#calculo del tri a partir de altitud con libreria spatialEco
-
-tri.ext <- tri(alt.m)
-projection(tri.ext) <- crs.geo 
-trivalues_arbol<-extract(tri.ext,dbroteri_arbol)
-
-presvals_arbol <- presvals_arbol[,-38]
-presvals_arbol <- cbind (presvals_arbol, trivalues_arbol)
-presvals_arbol <- presvals_arbol[,c(1:37,46,38:45)]
-colnames(presvals_arbol)[38] <- "tri"
-
-#BACKGROUND
-dbroteridata_arbol <- as.data.frame(dbroteri_arbol)
-presvalsdata_arbol <- cbind (dbroteridata_arbol[,c(3,4)], presvals_arbol)
-presvalsdata_arbol <- presvalsdata_arbol [,-48]
-presvalsdata_arbol <- presvalsdata_arbol [,c(2,1,3:47)]
-presvals2_arbol <- presvalsdata_arbol
-coordinates(presvals2_arbol) <- ~long+ lat
-proj4string(presvals2_arbol) <- crs.geo
-
-backgroundcoord_arbol <- presvalsdata_arbol [,c(1,2)]
-coordinates(backgroundcoord_arbol) <- ~long+ lat
-proj4string(backgroundcoord_arbol) <- crs.geo
-
-
-mask <- raster(backgroundcoord_arbol)
-res(mask) <- 0.008333333
-x <- circles(backgroundcoord_arbol, d=50000, lonlat=TRUE)
-#Se podria hacer un clip de los poligonos y el continente para que no salgan puntos en el mar, solucion provisional aumentar el N
-pol <- gUnaryUnion(x@polygons)
-samp <- spsample(pol, 1000, type='random', iter=25)
-extent(mask) <- extent(pol) # Sirve para que las submuestras de los poligonos salgan en el extent de la muestra
-cells <- cellFromXY(mask, samp)
-length(cells)
-cells <- unique(cells)
-length(cells)
-xy <- xyFromCell(mask, cells)
-# A los puntos generados para el background unimos los de presencia para este citotipo (util para que no de error la funcion posterior ecospat.grid.clim.dyn)
-xy <- as.data.frame(xy)
-colnames(xy) <- c("long", "lat")
-
-# Quitamos los puntos en el mismo km^2
-coordinates(xy) <- ~long+ lat
-proj4string(xy) <- crs.geo
-r <- raster(xy)
-res(r) <- 0.008333333
-r <- extend(r, extent(r)+1)
-backcoord_sel <- as.data.frame(gridSample(xy, r, n=1))
-coordinates(backcoord_sel) <- ~long+ lat
-proj4string(backcoord_sel) <- crs.geo
-
-# Extraccion de variables para el background y modificacion de la tabla
-backgroundclim <- extract(variables,backcoord_sel)
-backgroundsoil <- extract.list(backcoord_sel, list.files("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas"),path = "D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2017_nicho/soilgrids/capas", ID = "ploidy")
-backgrounddat <- cbind("background",as.data.frame(backcoord_sel),backgroundclim, backgroundsoil)
-backgrounddat <- backgrounddat[,-42]
-
-coordinates(backgrounddat) <- ~long+ lat
-proj4string(backgrounddat) <- crs.geo
-trivalues1 <- extract(tri.ext,backgrounddat)
-backgrounddat <- as.data.frame(backgrounddat)
-backgrounddat <- backgrounddat[,-40]
-backgrounddat <- cbind (backgrounddat, trivalues1)
-backgrounddat <- backgrounddat[,c(1:39,50,40:49)]
-colnames(backgrounddat)[40] <- "tri" # Sustitucion de los valores tri por los del dataframe con NAs
-
-backgrounddat.c <- na.omit(backgrounddat)
-backgrounddat.c <- cbind(backgrounddat.c,apply(backgrounddat.c[,c(42:44)], 1, mean))
-backgrounddat.c <- backgrounddat.c[,-c(42:44)]
-colnames(backgrounddat.c)[48] <- "AWC"
-backgrounddat.c <- backgrounddat.c[,c(1:41,48,42:47)]
-backgrounddat.c <- backgrounddat.c[,-48]
-backgrounddat.c <- backgrounddat.c[,c(2,3,1,4:47)]
-colnames(backgrounddat.c) <- colnames(presvalsdata_arbol)
-
-todo <- rbind (presvalsdata_arbol, as.data.frame(backgrounddat.c))
-
-todo.pca <- todo[,-c(1:3)]
-selected2 <- vif_func(todo.pca)
-todo.pca.2 <- todo.pca[,c(selected2)]
-
-# presvals.pca.corselect2 <- cbind (todo.pca,1)
-# correlations2 <- corSelect (presvals.pca.corselect2, var.cols = 1:44, sp.cols = 45, cor.thresh = 0.75)
-# presvals.pca.2.corselect2 <- presvals.pca.corselect2[,correlations2$selected.var.cols]
-
-todoploidy <- factor (todo$ploidy, levels = c("2x", "4x", "6x", "12x", "background"), ordered = TRUE)
-w<-c(rep(0,nrow(presvalsdata_arbol)),rep(1,nrow(as.data.frame(backgrounddat.c))))
-
-pcaback <-dudi.pca(todo.pca.2, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
-gcol = c("blue", "red", "green", "purple", "black")
-s.label(pcaback$li, clabel = 0.1)
-scatter(pcaback, clab.row = 0, posieig = "none", cex=0.1, clab.col = 0.5)
-s.class(pcaback$li, todoploidy, col = gcol, add.plot = TRUE, cstar = 0, clabel = 0, cellipse = 1.5, pch = 16)
-legend (7.5,-1.8,c("2x", "4x", "6x", "12x","Background"), col = gcol, pch =19, text.width = 1.8, y.intersp = 0.5, cex = 0.8)
-
-muestras<-todo$ploidy_arbol!="background"
-PCAphylo<-pcaback$li[muestras,]
-saveRDS(PCAphylo,"PCAphylo.RDS")
-
-write.table (pcaback$co, "PCAvars_phy.txt", dec = ",", sep = "\t")
diff --git a/rotation.tsv b/rotation.tsv
deleted file mode 100644
index 00f48b8..0000000
--- a/rotation.tsv
+++ /dev/null
@@ -1,45 +0,0 @@
-"PC1"/t"PC2"/t"PC3"/t"PC4"/t"PC5"/t"PC6"/t"PC7"/t"PC8"/t"PC9"/t"PC10"/t"PC11"/t"PC12"/t"PC13"/t"PC14"/t"PC15"/t"PC16"/t"PC17"/t"PC18"/t"PC19"/t"PC20"/t"PC21"/t"PC22"/t"PC23"/t"PC24"/t"PC25"/t"PC26"/t"PC27"/t"PC28"
-"bio1"/t-0.224328935738146/t-0.108135933667616/t-0.0192166317771092/t0.0811585504650087/t0.00213208856911033/t0.0317203266434497/t-0.0075186051395045/t-0.0563904169627873/t0.0438870448071573/t-0.028122829487266/t0.0204946375787728/t0.00868629198882989/t0.0770010491110456/t0.067920896667818/t0.0759087690549226/t0.00855492744460884/t-0.0595180690790423/t0.0460642799418226/t-0.0560790458936688/t0.136969505125186/t-0.0115267919190899/t0.180050955391694/t0.163811455893673/t-0.144437599855246/t0.196803789701143/t0.0498448780599188/t-0.167061767580596/t0.170149064844788
-"bio2"/t0.139964068004947/t0.0738837073124654/t-0.227352637807745/t0.236826341981264/t-0.0286607072966779/t0.0290275338256475/t0.0746676391866181/t0.0257631352996279/t-0.0834931967133748/t-0.10801405483137/t0.0679800402147643/t0.0663068811515626/t0.126420222846291/t0.0591452869274151/t-0.098827656628493/t0.0226477263264403/t-0.000199326021780646/t0.0739210172717185/t-0.0495917925858231/t-0.0618478600156289/t0.0257171221295562/t-0.0361177425912496/t0.036018934641065/t0.144604517473474/t0.221713627083418/t-0.273102576095056/t0.0796916748208247/t0.187546339312862
-"bio3"/t0.0528886098576754/t0.20130071676867/t-0.163231368325188/t0.265821553880965/t0.059720105306165/t-0.133494310263219/t0.240719824256872/t-0.0472710083905515/t-0.124861757641513/t-0.234178196332056/t0.0271123029264932/t0.0467242577878455/t0.284765695445983/t0.027932830211907/t-0.23528018775289/t-0.023525577113862/t0.316713899841813/t0.0496718604105513/t-0.247590770213207/t0.158100676075622/t0.0470686652163026/t0.183049460540892/t-0.247831942724228/t-0.13095089335859/t-0.184444604756347/t-0.123671616285522/t-0.269165660626522/t0.0370077939384998
-"bio4"/t0.175575901517895/t-0.116337710679734/t-0.151505337928262/t0.16315996779917/t-0.0888258497951207/t0.195943654165896/t-0.113553673132863/t0.0406707898029753/t-0.0250409178515241/t-0.0132894237586358/t0.031787573143816/t-0.013109791399167/t-0.0367696949633614/t-0.0343503467506036/t-0.0347292699874004/t-0.00440130428746535/t0.00913817033241563/t-0.00727785539383569/t0.102355233244006/t-0.181071607175102/t-0.00578862624381892/t-0.119307534329242/t-0.0227810132864881/t0.086033469668373/t-0.0559263165452369/t0.334489798632357/t-0.231919739144149/t0.0787335973554748
-"bio5"/t-0.0228524841740174/t-0.104201969550059/t-0.245362901080693/t0.328751420028439/t-0.0769391375702561/t0.14781941293015/t-0.0126439693932101/t-0.0188097251510023/t0.000800677024635526/t-0.0404541427133232/t0.0479542562830434/t0.0708299178324132/t0.138249088474119/t0.124804593812749/t0.0276164442782185/t0.00643357175851452/t-0.100582156582138/t0.0701564071564559/t-0.00106898443403391/t-0.0616224618227148/t0.0132024074643105/t-0.068225194109358/t0.209089317190267/t0.0389962413492879/t0.124114793653623/t-0.113873901853235/t0.371864099897156/t-0.340461174519786
-"bio6"/t-0.224297312773061/t-0.0438370005091041/t0.113536961439594/t-0.0494898241972046/t0.0236155206340591/t-0.0542981957602048/t0.0235156261902036/t-0.0677996148266875/t0.0786391486789488/t0.0396413077707623/t-0.0310701247058464/t-0.00345146923423589/t0.0387352508466192/t0.0592835178729229/t0.0975965663108616/t-0.0149888487531879/t-0.0341260047123138/t-0.00931828337185415/t-0.0330662348292573/t0.129027844158091/t-0.0590854605855679/t0.0692975525415152/t0.201708353144527/t-0.100419924508101/t0.00828932586489826/t-0.0892479396135318/t0.137627441426391/t0.23422249513986
-"bio7"/t0.165189456058806/t-0.0177718933702543/t-0.213534124233689/t0.205715985408034/t-0.0578314011473431/t0.116273533903968/t-0.0239428618981012/t0.0394155507771865/t-0.0618612307504168/t-0.0546369468619572/t0.0479946391093488/t0.0384095763666996/t0.0404029498913832/t0.029859496054149/t-0.0649024407648289/t0.0100528722106082/t-0.0183462256480998/t0.0431954434706819/t0.0300310495047816/t-0.120479001008693/t0.0441019971704605/t-0.11219918598983/t-0.0292788782813024/t0.0807963642409813/t0.0973637672633268/t0.0245888101190404/t-0.0431558044038454/t-0.0989580598256977
-"bio8"/t-0.0807275692132175/t-0.222888950633431/t0.221638680071944/t-0.0197212353533235/t-0.0414756871392428/t0.0437162315154033/t0.00639918524815399/t0.00920738372450987/t0.116385884884644/t0.0478302483013544/t0.0359287581606769/t0.0417487269718488/t0.0385598731670806/t-0.317475410387318/t-0.575948058097652/t-0.0703846124803775/t-0.00295344901208148/t-0.0187030065295181/t-0.329517755333647/t-0.0210142925143645/t-0.412058338899018/t-0.117167185000066/t0.085415756367112/t0.187442406441073/t0.254878855549771/t0.0508058090716903/t-0.0557696965039789/t-0.0251309751016619
-"bio9"/t-0.136008826607553/t-0.175901439979415/t-0.178915575965478/t0.190238104767456/t-0.00706001255162697/t0.125358600895453/t-0.055377428887669/t-0.095437550569745/t0.0114399492100258/t-0.0191743858211837/t-0.00659723991817387/t0.0647366439252747/t0.0152360171092358/t0.0758553980382982/t0.040579528032953/t-0.0446746584053742/t-0.0740150132614736/t-0.0167412703781161/t-0.165689696823577/t0.0457475915180974/t-0.0666017959847031/t0.0776183388938223/t-0.25692999671899/t-0.0322651092469687/t-0.0353577960880398/t0.240424294430528/t0.254400330197027/t-0.00364302305903496
-"bio10"/t-0.122087555388051/t-0.211815558576371/t-0.121312559502315/t0.213249241752894/t-0.0588588618105384/t0.192633509130845/t-0.100771565745341/t-0.0319984745434518/t0.037064954754127/t-0.0302729107357906/t0.0401568532696747/t-0.00288884572487647/t0.0605299010048209/t0.0420859851465468/t0.0513352567327501/t-0.0144823119357991/t-0.0502190569279166/t0.0438175387587519/t0.0108858510248177/t0.0252028668273325/t0.00779003077807668/t0.0381276685263598/t0.0724160512125461/t-0.0890638823732715/t0.0299073320935973/t0.265389769278108/t-0.0231088541321337/t0.543341459923896
-"bio11"/t-0.235322198031178/t-0.0445599663853468/t0.055468473970033/t-0.00458649097812645/t0.0354135579688872/t-0.0435963503832855/t0.03081883985737/t-0.0579359635966208/t0.0449708028225693/t-1.29298234103972e-05/t0.0080137961283914/t0.0103564370736636/t0.0673401590064133/t0.0572568304152086/t0.0660805476770162/t0.00044084163168696/t-0.00774797585315138/t0.03369133817178/t-0.0738393866619167/t0.154294847157002/t-0.0283079775944016/t0.11059735040465/t0.0289256617090446/t-0.0941380672050338/t-0.0782510828917371/t-0.104206629565338/t0.121392621636845/t-0.0028062391566825
-"bio12"/t-0.0725307438171991/t0.28593680697532/t0.0110582912296724/t0.0920821548378869/t-0.203962471570015/t0.0850557650757668/t0.159493248746905/t-0.0850020140985716/t0.149222537807831/t0.0391491969700159/t0.121433816704082/t-0.098741653587019/t-0.132252987499037/t0.0657749555689423/t0.0274596910598766/t0.0436614200954679/t-0.458923930913251/t-0.15737426180952/t-0.119133615677169/t0.151501625948218/t0.168002665082231/t-0.0985333871528145/t-0.152747095781552/t0.0402140502181371/t0.180007299821801/t-0.00724334827510461/t-0.0653543493628536/t-0.0938882647860404
-"bio13"/t-0.0920806780726117/t0.269790657025488/t0.0294750434746603/t0.100334213195544/t-0.151807296111834/t0.200596124892898/t0.0629013903725691/t-0.0269597085764402/t0.207524351227219/t0.111168403493979/t-0.0298142310138828/t-0.301560221315984/t-0.0588757535587367/t-0.0413795265068081/t-0.208283835457058/t0.0064348239227133/t0.0753959528197465/t0.0451094081924796/t0.128185575494008/t-0.130811067633469/t-0.181948126156292/t0.118838045669719/t-0.20046065985114/t-0.0760300393963822/t-0.332734881862752/t0.0270318359966738/t0.196611117606505/t0.10252407129687
-"bio14"/t0.164034541994682/t-0.0936365433865695/t0.218008007318956/t0.11074489277833/t-0.0941510016185772/t0.0184930493325922/t0.122372304074818/t-0.105400733993421/t0.132545654413691/t-0.0982089611269703/t-0.0428535139196642/t0.0488432138178762/t0.134633488791768/t-0.0377810516400632/t-0.0570949918691031/t-0.0808016675338867/t-0.0263276150725949/t0.00460550011931616/t0.212279884581085/t-0.0297857839273159/t-0.0423398334941024/t-0.163423203291482/t-0.0339714809488422/t-0.0642408363553065/t-0.33943264799207/t0.153031641970418/t0.366863963911568/t-0.0380902006502719
-"bio15"/t-0.160330127507147/t0.203024514747095/t-0.113852570009908/t-0.0211339965773361/t0.0963801117901862/t0.0982627057153258/t-0.125237832981925/t0.0458809028708733/t0.012242988887706/t0.112423910709578/t-0.106010837806472/t-0.266870336073056/t-0.00544210183540184/t0.0122190558711746/t0.0112672153725288/t0.239340923490242/t0.436307755474952/t0.371052454822368/t-0.0429688126276107/t0.104718401251493/t0.0820895499226338/t-0.0668735129076415/t0.317134036758452/t0.229061267533267/t0.0061657695910496/t0.219079104569069/t0.0656841679544706/t-0.0951794870982389
-"bio16"/t-0.0912274025497038/t0.283755849605586/t0.000571140233341816/t0.103415362338651/t-0.128964439378172/t0.143482833111823/t0.0400184191024118/t-0.0102635246027168/t0.180731766901472/t0.0870483278500592/t0.01181154990115/t-0.231724484755862/t-0.0839148988006927/t0.0145339750145121/t-0.146355367812479/t0.0328781627183407/t-0.05606815640917/t0.00940224546172848/t-0.0451753848313089/t0.0481999416189897/t-0.0265672389812509/t0.194970401154086/t-0.065862936034582/t-0.0453509052441023/t0.103862154640621/t0.0691563601475962/t-0.0658145320651485/t-0.0879458163701042
-"bio17"/t0.165835300326787/t-0.0902046008635743/t0.217636455987803/t0.116038313813937/t-0.113899328371588/t0.064613152381725/t0.0822092227164777/t0.0133634180554367/t0.0861111636876059/t-0.0520577836274791/t0.0263813184519956/t-0.00355379510690291/t0.05579060321815/t-0.117182565198934/t-0.134710648265116/t-0.0267028544166423/t-0.0798219327649053/t0.0270553677474192/t0.185620812464909/t0.0801651324278929/t0.305777317551849/t0.155796177523281/t0.293596918189224/t-0.0143018458205006/t-0.0176127366529004/t-0.111678680553097/t-0.0757341038185484/t0.064349694392251
-"bio18"/t0.162048672915887/t-0.083989111551446/t0.224223537646862/t0.114283883741873/t-0.108456206652718/t0.090430169927621/t0.06379708362473/t0.00159585049685929/t0.116911636273258/t-0.0935242862367341/t-0.0128558077086383/t-0.0509228434053479/t0.0624127942859179/t-0.168545395462505/t-0.119798959653117/t-0.0251710044191471/t0.00436183897219739/t0.0408354459298768/t0.0823270602009548/t0.171459247441223/t0.373030110206253/t0.216640985467439/t0.207642018029765/t0.0167168018841127/t-0.0265103827881532/t0.0428946865601112/t-0.0467157498024832/t-0.0906470142612253
-"bio19"/t-0.117310550522094/t0.268703870385105/t-0.0989714403599621/t0.0610266958006039/t-0.0587132238305808/t0.0127516529893893/t0.070705046143366/t-0.0992581972529779/t0.0610088964632306/t0.0877464965439966/t0.0523709101618525/t-0.0975149612000122/t-0.0530272635319355/t0.0809325168331301/t-0.0134953559975252/t0.0999795849635975/t-0.157680332041878/t-0.130666660610488/t-0.00790227405656937/t-0.177928834550906/t-0.0318961245975375/t-0.278737543544228/t0.339813924220854/t-0.0252267267929255/t-0.0670441505394695/t-0.1439289532194/t-0.170472559743587/t0.151645029301563
-"annualPET"/t-0.0377985121773322/t-0.129320370640655/t-0.276604707266857/t-0.192829228806995/t-0.158154617564522/t-0.0954046273653742/t0.119762153447739/t0.0510726402583764/t-0.0142318688592219/t0.255266776223223/t-0.0254236394343699/t-0.0481125254538953/t0.131416454704796/t-0.0467261079384709/t-0.1209264632335/t-0.149317393742887/t-0.0997073641074179/t0.0554312381664394/t-0.140611750707065/t-0.126944386968083/t0.091448194411322/t0.013544230713913/t0.0927409909157001/t0.0413680701294197/t-0.312960531777297/t-0.285143638228478/t0.106695660284358/t0.206603196464808
-"aridityIndexThornthwaite"/t-0.16225755299065/t0.0410143014467005/t-0.213797932519977/t-0.0783771809010973/t0.100465253925717/t0.0802898640941684/t0.0251172403399018/t0.183506166434245/t0.0065803131318867/t-0.182444490359371/t0.240307893699651/t0.352303363570871/t-0.424248978471339/t0.202886828016964/t-0.293323695832166/t-0.0970547807534873/t-0.0408791558553349/t-0.0885301654218743/t0.144766387577896/t0.238242770924547/t-0.0536046511868054/t0.161260762030618/t0.215726769408404/t0.234631885085431/t-0.307672556134853/t0.0426006732665687/t0.00296880982220487/t-0.0418351960698094
-"climaticMoistureIndex"/t0.00580234338257277/t0.280923295108522/t0.0955856648850828/t0.0897223092546453/t-0.00138131289756106/t0.100657112374963/t-0.0846038339268504/t0.487416341075031/t-0.0792271812107403/t0.0888640163917425/t-0.111113692426787/t0.172096397017067/t0.0898794249045964/t-0.222001708081931/t0.0453432455009901/t-0.106704569888866/t-0.0466529108078469/t0.137707106034411/t-0.0614321472108247/t-0.0262116787389872/t0.0196381038280695/t-0.14993602636701/t0.0428125469754819/t0.0421469794002272/t0.0225871788029835/t-0.194864132575898/t0.193525152394791/t0.145823893935445
-"continentality"/t0.176183153438376/t-0.110519274944742/t-0.153505935943669/t0.070884726120517/t-0.128628947751001/t0.118411433913068/t-0.243516735096778/t0.00593024596500284/t0.0404503042380969/t0.323310203318892/t-0.254436660592253/t0.0203690286639744/t-0.111199972144615/t-0.0414493098383781/t0.0224812634692863/t-0.100138766394176/t0.0496473429300255/t-0.10334783365737/t-0.00536405449167938/t0.0726712779854823/t-0.0845030158974389/t0.20525470324539/t0.0375017198501386/t-0.0863103942965513/t-0.090113459489324/t-0.210609928263733/t-0.0826791842103903/t-0.0664892687312569
-"embergerQ"/t-0.0775280110263514/t0.267845441916043/t0.103963242024708/t0.0593442063842523/t0.0268582317874724/t0.0387488290214487/t0.0166632862494337/t0.455077573252199/t-0.0895838960695103/t-0.0269775500832536/t-0.0339798264553154/t0.224779804192738/t0.110417679962971/t-0.121052194285279/t0.143078090453554/t-0.129053039898774/t-0.202565749316273/t0.125302357761192/t0.0157302225167563/t-0.0502864904524643/t-0.134130884613037/t0.213635300777435/t-0.00524885016073674/t-0.18863660559636/t0.0263254919161838/t0.268852132535534/t-0.045889683796431/t-0.0242437453229883
-"growingDegDays0"/t-0.227676045848928/t-0.099180640451683/t-0.0337052346763953/t0.0205699307809611/t-0.0523315627947985/t0.0402849812001279/t-0.0460009425437976/t0.0595065475383176/t0.0130133070870755/t-0.00117548521571254/t-0.0123939094999675/t-0.000571301003641484/t0.113698051696388/t-0.0587610526663474/t0.0565324723965495/t0.00745040906408579/t-0.0421529709442394/t-0.0271910275863506/t0.0411378532541797/t0.130641977714854/t-0.0462611310373157/t-0.0504210901041321/t-0.131210801348564/t0.111101889345834/t-0.0847519546702425/t-0.127403880933037/t-0.0347417369058467/t-0.316013934360502
-"growingDegDays5"/t-0.225697128863733/t-0.0929741289309785/t-0.0508194665932009/t0.00140499937891402/t-0.0665599155353057/t0.0658240559898157/t-0.0162787524590744/t0.0875236328162942/t0.0662562293430908/t-0.106899418902238/t-0.0688031577735141/t-0.0217456496637579/t0.0606470019070226/t-0.0860206412382681/t0.0821324946000522/t-0.0213538830990788/t0.0201217394830918/t-0.0834213026771358/t0.049623563146613/t0.101326357613049/t-0.033351767120729/t-0.047260348501764/t-0.133403893664755/t0.110402547266662/t-0.0268623864371646/t-0.130090442875561/t-0.0155819714542014/t0.195853238729979
-"maxTempColdest"/t-0.221316317199991/t-0.102499116959946/t0.0229633843323997/t-0.0899764238510021/t-0.009703956673697/t0.0229352532358427/t0.115987002600082/t0.113173063895536/t-0.0264859150878108/t-0.107546857503829/t0.145442032356096/t-0.0356134910185745/t0.0938922551902227/t-0.0729181298010529/t-0.0780541198660653/t0.0370312187960912/t-0.0188814875528736/t0.0900586124386688/t0.0995015963528143/t-0.179278411736647/t0.232297760688978/t-0.366341727520588/t-0.0762868504550143/t-0.0421421407107268/t-0.176230549651899/t0.0867616430821745/t-0.206796297274622/t0.0620801166443522
-"minTempWarmest"/t-0.159373419961284/t-0.179064523899487/t0.0360477670564315/t0.159063135922672/t-0.0562232111761212/t0.257156894462404/t-0.295381977917926/t0.0992938647226184/t0.055257728950723/t0.120050367490029/t-0.111889489764691/t0.00560290885765775/t-0.0204154065124398/t-0.0850309714910451/t0.0839242236171106/t0.023640689847799/t0.00905272323630068/t0.00212869960800026/t0.0428862671396123/t0.0874862647916598/t0.0617503506152553/t-0.036587664817085/t-0.137461861773598/t0.00428836319426407/t-0.164424142859599/t-0.224453832148598/t-0.366575628231777/t-0.13595192113554
-"monthCountByTemp10"/t-0.233670755100617/t-0.0557028942577398/t-0.0183044560057545/t-0.0121205842388315/t-0.020160872190532/t-0.0109588063301384/t0.106946478847715/t-0.0642967940333929/t0.0501100811301048/t-0.0515339542472987/t0.0553355896686774/t0.0578885223721612/t0.184439445942949/t0.0583146049922012/t-0.0586168994442327/t0.109593040413892/t-0.00214321830307859/t0.144798809844614/t0.43982424411085/t-0.43834133764723/t-0.249595627683499/t0.320341082546053/t0.0788089722649358/t-0.0984230894256262/t0.150169410918621/t-0.16363855336832/t-0.12082680570163/t-0.148204450084268
-"PETColdestQuarter"/t-0.150860718839305/t-0.164111131346627/t-0.136991191136438/t-0.218074762095281/t0.0737072537491323/t0.0670123720653381/t0.0529869984629142/t0.0942239440331043/t-0.00205303653131963/t0.179474974929775/t0.091211892672495/t-0.023995130171473/t-0.161359296283338/t-0.0225687494263269/t-0.204475659795453/t-0.110556064450031/t-0.0466254210784407/t0.233270445111606/t-0.200011513811157/t-0.250876112762453/t0.469300867568603/t0.15781541945001/t-0.155167712618213/t-0.220336098070738/t0.137100242864328/t0.0700060973079136/t0.124910962151422/t-0.106562363269222
-"PETDriestQuarter"/t-0.0637717423579856/t-0.0332658448129387/t-0.25312737066886/t-0.177494974234045/t-0.188208691444502/t0.0727676279989646/t0.209627449379399/t0.17780523393079/t0.276934752384958/t-0.468379170735816/t-0.331456109830669/t-0.133873221178143/t-0.133182292470815/t-0.185373889827358/t0.160606527580051/t-0.194566586832141/t0.237923891618433/t-0.292331194385711/t-0.0237090301118913/t-0.108724962930673/t0.0362111384919226/t0.00985140587179486/t0.0651198958389308/t-0.010804246975264/t0.174727649310794/t0.000379994594272451/t0.0457445893788213/t-0.032462324937663
-"PETseasonality"/t0.0916022715118237/t-0.00838577592316629/t-0.279681387438799/t-0.128060397906005/t-0.230227923179474/t-0.193919550449045/t0.140393553174271/t0.0106557928609329/t-0.0345479225461165/t0.155106514142379/t-0.0796337720737982/t-0.0263605560008455/t0.2808935580476/t-0.0574169208083356/t0.0697940036440614/t-0.0882025400347951/t-0.183986397788584/t-0.00859674416450782/t-0.143775516426568/t0.102670142484407/t-0.147362197815107/t0.0127876405168721/t0.29557972055394/t-0.108921314781534/t-0.187374682101734/t0.333340805528522/t-0.184642539095961/t-0.227207637025508
-"PETWarmestQuarter"/t0.0582814935856609/t-0.0108941891335254/t-0.299306903234295/t-0.184886344922517/t-0.181537004487034/t-0.143013323940105/t0.0905805814127839/t0.0579584686439406/t-0.0238367024586243/t0.171653543755891/t-0.0259518671167875/t0.0413012260087688/t0.0311085081340637/t-0.197079283013722/t-0.170749836964206/t0.202184984669121/t-0.0306166487646765/t0.130585353579856/t0.493564906115165/t0.415463289342805/t-0.0269208729936465/t-0.0828362825858246/t-0.234273654186078/t-0.0188484579851025/t0.260711244139359/t0.0353798733572134/t0.0659193062471091/t0.124184765885302
-"PETWettestQuarter"/t0.0799026812628079/t-0.219730012459824/t0.120764918266885/t0.129269565501614/t-0.0795756735746461/t0.0433900191762119/t0.535595753957845/t0.0727343659889786/t0.0475916264660216/t0.0549581959176639/t-0.346028995580876/t0.135828733630064/t-0.349181224238744/t0.185430037791358/t0.132481130529805/t0.21159902834898/t-0.0545779533084998/t0.416829171223472/t-0.124369903819554/t0.0169319753674072/t-0.141074951765974/t-0.0778229727474098/t-0.043057332761674/t0.0371532320727768/t-0.0484741070080367/t-0.0551275620412566/t-0.121289313272932/t0.0144569354097411
-"thermicityIndex"/t-0.236127463593907/t-0.053912879490521/t0.0212593269067196/t-0.0112121839643976/t-0.0131059511751412/t0.0122713162814375/t0.0420413166049116/t0.0674254929567633/t-0.0116730624222692/t-0.108627991918002/t0.0828975256571634/t-0.00950677494377056/t0.129807301855702/t-0.0381372306045542/t0.0410002543561623/t0.0528993227528436/t-0.0548368686693341/t0.0401906095822468/t0.0423898583326168/t0.0062288566231805/t0.0273704936711147/t-0.177280215522779/t-0.0489039730573515/t0.0568504585783085/t-0.120509200054357/t-0.0234274920021169/t-0.0984095520104995/t-0.166573472059069
-"topoWet"/t-0.201715188550227/t-0.0261796724184499/t0.0603701475994947/t0.180224880668954/t-0.121580058662333/t-0.0851745165991217/t0.0749077735863486/t0.102067970256619/t-0.275142808786501/t0.154128709508179/t-0.206149694163368/t0.204994016284862/t-0.0643296921167888/t-0.0629472986617315/t-0.175180054571128/t0.508819545366736/t0.154282116589396/t-0.490482977539483/t-0.0323564209714221/t-0.0964777334518348/t0.154052131933524/t-0.00065136317605818/t0.0745504084471151/t-0.208829547537192/t-0.0195674321770304/t0.10538102314685/t0.127780451109565/t-0.0290691090548688
-"tri"/t0.155833040940001/t0.0185788813967378/t-0.148457491787441/t-0.221189610164462/t0.132123896541521/t0.187499211729266/t-0.217385427133626/t-0.108512274085223/t-0.00894989674905656/t-0.405177166795009/t-0.0615576331498494/t-0.0195399084353021/t-0.0235818951250484/t-0.299221179367911/t-0.0380611925465029/t0.507451552265232/t-0.327606780554933/t0.1673306909143/t-0.178474642717886/t0.0137875469059488/t-0.0819433867344176/t0.0825393878623676/t0.0447921167332966/t-0.198992262384866/t-0.174289332066324/t-0.111822427516273/t0.0627856400344508/t0.0123735682239494
-"elevation"/t0.227397694909475/t0.0680974277243243/t-0.0727638914266795/t-0.06066428749251/t0.0335641253865587/t0.0152070694182943/t-0.0519391462099498/t0.0158393427936218/t-0.0648895257177605/t0.0563948806098457/t-0.0158969811811115/t0.0443966904076983/t-0.222251441754658/t0.0275915595799931/t-0.0417504628872487/t-0.0239109272949297/t0.00561397942694227/t-0.0763329370403259/t0.0687659299421901/t-0.271236832346286/t-0.0727465892957934/t0.0786151806676599/t-0.0341115445463991/t-0.0504061041590598/t-0.00659289275550947/t-0.00125095589394334/t-0.128304040303232/t0.0670954409776237
-"AWC"/t0.127418131471096/t-0.162349945024699/t0.110154210563711/t-0.0763953365361789/t-0.280257587124581/t-0.0594120415302956/t-0.192601558038028/t0.422195678261961/t-0.0703389742459672/t-0.178904981134533/t0.136110861511861/t-0.331847790031113/t-0.0139028455807081/t0.482281505448491/t-0.171573804768705/t0.053213387527198/t0.0460208283715072/t0.0246068807592731/t-0.0375408786803816/t0.119358143187406/t-0.132263170966367/t-0.122244585361775/t0.0311640727085706/t-0.359624608047802/t0.0467092039820658/t-0.0916761800545616/t0.028023891512895/t-0.0327049415094918
-"BLDFIE"/t-0.111153749600858/t-0.0482771200225892/t-0.101123111150358/t0.335561494755849/t0.26429795304239/t-0.183498703778612/t0.0497200021667333/t-0.0863750233720208/t-0.267061532621259/t-0.0489415261073795/t-0.0130502221639374/t-0.282420541858149/t-0.356114883358267/t-0.332897482622503/t-0.0431945099077874/t-0.290186683887773/t-0.0562580681486159/t0.0810248124969172/t0.146441848039126/t0.102347812191302/t-0.0922064780446633/t-0.175685612794837/t0.126574636433/t-0.361349909166683/t0.00369743632070352/t-0.0631634517190948/t0.0111023791312974/t-0.0831380902376236
-"CECSOL"/t0.0347564706302964/t0.0230788210533776/t-0.000405255622623181/t-0.203422630301272/t0.428314276312698/t0.504828292531205/t0.205139899894017/t-0.0263526816306639/t-0.192255440437848/t0.0634040363412119/t-0.357295024599633/t-0.0909078622964542/t0.255988894931676/t0.254528597289147/t-0.212586942290355/t-0.135253249698069/t-0.115295516891011/t-0.170117952930023/t0.127705592172835/t0.130445674232739/t0.00120709257492445/t-0.105697208913595/t0.0298010433561572/t-0.085489237848571/t0.0664797553796919/t0.00817721270981004/t-0.0287806188145003/t-0.0201492661504882
-"ORCDRC"/t-0.00714747454363387/t0.0188189872531611/t0.102400203611736/t-0.191131632043013/t-0.381513286668628/t0.439530218108651/t0.196599433053599/t-0.170519065657104/t-0.532058290533785/t0.00180268026991911/t0.329689816203732/t-0.00132155276458602/t-0.0776961685957538/t-0.164165654184379/t0.205373880266806/t-0.0731576098778127/t0.193719255120277/t0.0138787635464308/t-0.0530028905684151/t0.0587902685259578/t-0.112038324542497/t0.0614405723425727/t0.0417782909671735/t-0.020417718485655/t0.0159999970223876/t-0.0623193156839113/t0.0422487070614047/t-0.00988804594776375
-"PHIHOX"/t0.0937805011232968/t-0.230682971198385/t-0.0224130285744476/t0.105313272066117/t0.26707510622708/t-0.0498528696065704/t0.237177337258313/t0.340654511127417/t-0.125316034362824/t0.0900501276420434/t0.206031791303369/t-0.46680731284736/t0.0377500403320047/t-0.103224513598306/t0.148387202598381/t0.199568234116294/t-0.146759189446147/t-0.213830776169051/t-0.00333487538128612/t-0.0252354190062977/t-0.0706607677237541/t0.275277341901093/t-0.0330721409185081/t0.338807582414083/t-0.0658500678836616/t0.00773267162599152/t0.0476806436279749/t0.0397325218935467
-"SNDPPT"/t-0.167832039660242/t0.0497283303668456/t0.100934013125017/t0.0047906141872552/t-0.25528420945601/t-0.205504953699394/t-0.19722549176782/t-0.118832611597444/t-0.455405286960333/t-0.224977765382666/t-0.421993682338782/t-0.168898702111153/t-0.0352764922397226/t0.114545168251548/t-0.160543013030712/t-0.122862502160769/t-0.257602797334541/t0.115959635275861/t0.0480152604383469/t-0.0905111124294882/t0.0895188506772944/t0.150697445136912/t-0.0479619358583037/t0.352353048292371/t-0.0187155253755434/t0.0639360820174991/t0.0118710398647973/t0.0523069269724753
diff --git a/spThin.R b/spThin.R
deleted file mode 100644
index 8cbe5f0..0000000
--- a/spThin.R
+++ /dev/null
@@ -1,50 +0,0 @@
-#' @name spThin
-#'
-#' @title Spatial thinning of Presence-only Data
-#'
-#' @description Apply the stochastic spatial thinning algorithm implemented in the spThin package to presence data in a presence-background dataset
-#'
-#' @details Full details of the algorithm are available in the open-access article by Aiello-Lammens et al. (2015): dx.doi.org/10.1111/ecog.01132
-#'
-#' @param .data \strong{Internal parameter, do not use in the workflow function}. \code{.data} is a list of a data frame and a raster object returned from occurrence modules and covariate modules respectively. \code{.data} is passed automatically in workflow from the occurrence and covariate modules to the process module(s) and should not be passed by the user.
-#'
-#' @param thin Thinning parameter - the required minimum distance (in kilometres) between points after applying the thinning procedure
-#'
-#' @family process
-#'
-#' @author zoon Developers, \email{zoonproject@@gmail.com}
-#'
-#' @section Data type: presence-only
-#'
-#' @section Version: 0.1
-#'
-#' @section Date submitted:  2017-11-22
-spThin <- function (.data, thin = 1) {
-  
-  # check these are presence-background data
-  stopifnot(all(.data$df$type %in% c('presence', 'background')))
-  
-  # install & load the package
-  zoon::GetPackage('spThin')
-  
-  # get dataframe & index to presence data
-  df <- na.omit(.data$df)
-  pres_idx <- which(df$type == 'presence')
-  
-  # prepare presence data subset and apply thinning
-  sub_df <- data.frame(LAT = df$latitude[pres_idx],
-                       LONG = df$longitude[pres_idx],
-                       SPEC = NA)
-  th <- thin(loc.data = sub_df,
-             thin.par = thin,
-             reps = 1,
-             locs.thinned.list.return = TRUE,
-             write.files = FALSE,
-             write.log.file = FALSE)
-  
-  # get index to rows in sub_df, update the full dataset and return
-  pres_keep_idx <- as.numeric(rownames(th[[1]]))
-  .data$df <- rbind(df[pres_idx,][pres_keep_idx, ],
-                    df[-pres_idx, ])
-  return (.data)
-}
diff --git a/spatial_thin_log.txt b/spatial_thin_log.txt
deleted file mode 100644
index 5e6de3e..0000000
--- a/spatial_thin_log.txt
+++ /dev/null
@@ -1,26 +0,0 @@
-********************************************** 
- Beginning Spatial Thinning.
- Script Started at: Sun Oct 08 19:54:37 2017
-********************************************** 
- Beginning Spatial Thinning.
- Script Started at: Sun Oct 08 19:56:40 2017
-
-Thinning Parameter Used (in km): 1
-********************************************** 
- Beginning Spatial Thinning.
- Script Started at: Sun Oct 08 19:59:13 2017
-
-Thinning Parameter Used (in km): 1
-Number of replicates of thinning script: 100
-********************************************** 
- Beginning Spatial Thinning.
- Script Started at: Sun Oct 08 20:32:24 2017
-
-Thinning Parameter Used (in km): 1
-Number of replicates of thinning script: 10
-********************************************** 
- Beginning Spatial Thinning.
- Script Started at: Mon Oct 09 10:29:38 2017
-
-Thinning Parameter Used (in km): 1
-Number of replicates of thinning script: 10
diff --git a/tree.nex b/tree.nex
deleted file mode 100755
index 1d58dca..0000000
--- a/tree.nex
+++ /dev/null
@@ -1,32 +0,0 @@
-#NEXUS
-begin taxa;
-	dimensions ntax=22;
-	taxlabels
-	Zafarraya2
-	Zafarraya1
-	Orgiva
-	Laroles
-	SMiguel
-	Ifach
-	Cartagena
-	Alcublas
-	Azuebar
-	Lliria
-	Albufera1
-	Albufera2
-	Chiclana
-	Valverde
-	LagMadres
-	Acebron
-	Peladillo
-	Ronda
-	Comporta
-	Troia
-	BArportel
-	MClerigo
-;
-end;
-
-begin trees;
-	tree tree_1 = [&R] (((((Zafarraya2:0.022285533,(Zafarraya1:0.027696343999999998,(Orgiva:0.026052324,Laroles:0.020145377):8.400930000000001E-4):0.0023136480000000015):0.0018481340000000013,(SMiguel:0.022707511999999996,(Ifach:0.025276318,Cartagena:0.020187512000000005):0.0020802579999999946):0.0044702980000000045):6.267999999999968E-4,((Alcublas:0.021412639000000004,(Azuebar:0.01708129,Lliria:0.018697491):9.236420000000023E-4):0.0014286300000000002,(Albufera1:0.014228144000000002,Albufera2:0.01359442):0.006433400000000002):0.005309706999999997):6.213470000000013E-4,((Chiclana:0.023106505,(Valverde:0.026610595,(LagMadres:0.021779520000000004,Acebron:0.019812812):0.0018987879999999985):0.01089772):9.436229999999976E-4,((Peladillo:0.022100548999999997,Ronda:0.019762446999999995):4.444430000000027E-4,(Comporta:0.017539237,Troia:0.019412256000000003):0.004525745999999997):0.0014704499999999981):0.0014811200000000024):2.68034E-4,(BArportel:0.034570733,MClerigo:0.033750249):2.6803400000000065E-4);
-end;
diff --git a/zoon_server_javi.R b/zoon_server_javi.R
deleted file mode 100644
index aa87c56..0000000
--- a/zoon_server_javi.R
+++ /dev/null
@@ -1,99 +0,0 @@
-library (zoon)
-library (future)
-library (ulimit)
-library (gtools)
-library (spatialEco)
-
-#==================POPULATIONS=====================#
-
-# dbrot <- LocalOccurrenceData (filename="dbroteri.csv")
-# 
-e <- extent (-10,4.5,35.5,44)
-chefiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/vars_selected_gbif/chelsa", pattern = ".tif", full.names = TRUE))
-chelsa <- stack (chefiles)
-che.c <- crop (chelsa,e)
-
-alt15files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_15", pattern = ".bil", full.names = TRUE))
-alt16files <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/altitud_16", pattern = ".bil", full.names = TRUE))
-alt15 <- stack (alt15files)
-alt16 <- stack (alt16files)
-alt.m <- merge (alt15, alt16, ext=e)
-
-envfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/vars_selected_gbif/envirem", pattern = ".bil", full.names = TRUE))
-envirem <- stack (envfiles)
-env.c <- crop (envirem,e)
-
-soilfiles <- mixedsort (list.files ("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/vars_selected_gbif/soil", pattern = ".tif", full.names = TRUE))
-soilgrids <- stack (soilfiles)
-soil.c <- crop (soilgrids,e)
-
-tri.ext <- tri(alt.m)
-
-projection (che.c) <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
-projection (alt.m) <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
-projection (env.c) <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
-projection (soil.c) <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
-projection (tri.ext) <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
-
-combras <- CombineRasters(c(che.c, env.c, soil.c, tri.ext))
-vars.stack <- stack (combras [[1]], combras [[2]], combras [[3]], combras[[4]])
-writeRaster(vars.stack,"stack_zoon_alllayers.grd", format="raster")
-
-
-# Chain(PrintMap, PerformanceMeasures),
-
-LoadModule('PerformanceMeasures')
-
-vars.stack <- stack("D:/Copia de seguridad JAVI/UNIVERSIDAD DE SEVILLA/Experimentos Dianthus/Lopez_Juradoetal2018_nicho/stack_zoon_alllayers.grd")
-plan (multicore)
-
-work <- function(){ workflow (occurrence = LocalOccurrenceData (filename="dbroteri.csv",
-                                                                occurrenceType='presence/absence',
-                                                                columns=c(long = 'longitude', lat = 'latitude', value = 'value', type = 'type')),
-                              covariate  = LocalRaster(vars.stack),
-                              process = Chain(StandardiseCov,
-                                              Crossvalidate),
-                              model = MaxNet,
-                              output = PrintMap (dir = "/home/javlopez", size = c (600,600)),
-                              forceReproducible = TRUE)}
-
-workres<-work()
-save(workres, file = 'workflow_javi.RData')
-
-#==================POPULATIONS=====================#
-
-
-#==================GBIF=====================#
-
-LoadModule('PerformanceMeasures')
-LoadModule('ROCcurve')
-
-vars.stack.gbif <- stack("stack_zoon_alllayers.grd")
-plan (multicore)
-ulimit::memory_limit(100000)
-
-workgbif <- function(){ workflow (occurrence = LocalOccurrenceData (filename = "dbroterigbif.csv"),
-                              covariate  = LocalRaster(vars.stack.gbif),
-                              process = Chain (Background (n = 10000, bias = 100, seed = 999999), 
-                                               StandardiseCov),
-                              model = MaxNet,
-                              output = Chain (PrintMap (dir = "/home/javlopez", 
-                                                        size = c (600,600)),
-                                              PerformanceMeasures, ROCcurve (newwin = FALSE)),
-                              forceReproducible = TRUE)}
-
-workresgbif <- workgbif()
-save (output, file = 'workflow_gbif_nooutput.RData')
-
-#==================GBIF=====================#
-
-LoadModule ('ChangeWorkflow')
-
-ChangeWorkflow (workresgbif, occurrence = NULL,
-                covariate = NULL,
-                process = NULL,
-                model = NULL,
-                output = DataSummary)
-
-LocalOccurrenceData (filename = "dbroterigbif.csv", occurrenceType='presence', columns=c(long = 'longitude', lat = 'latitude', 
-            value = 'value', type = 'type', fold = 'fold'), externalValidation = TRUE)
\ No newline at end of file