################################# ### This script takes a list of predicted occupancy maps (coming from modelling procedures) ### calculates the niche overlapping among every two species ### and then applies a weighing according to centrality indexes ### obtained from a network (data for the network provided in other ### directories of this dataset) ################################# library(raster) library(sp) library(maptools) library(ggplot2) library(scales) setwd("~/Documents/Cosas/Networks/PathoMaps/HostandTicksModels_asc") ## point to the directory to write the results (replace as necessary) especies<-list.files(pattern=".asc") ################################# ### This is an example for Borrelia ### (users must to change list of files for other patheogens) ################################# Borrelia_ticks <- read.csv("~/Documents/Cosas/Networks/PathoMaps/Borrelia ticks-hosts.csv") ### the centrality indexes of Borrelia and ticks Borrelia_other <- read.csv("~/Documents/Cosas/Networks/PathoMaps/Borrelia paths-ticks-hosts.csv") ### the centrality indexes of Borrelia and hosts Borrelia_ticks[,1] <- paste(Borrelia_ticks[,1],".asc",sep="") Borrelia_ticks[,2] <- paste(Borrelia_ticks[,2],".asc",sep="") Borrelia_other[,2] <- paste(Borrelia_other[,2],".asc",sep="") # nombres <- as.vector(1,mode="any") acumula_ticks <- stack() acumula_hosts <- stack() #### We will plot reslts in the dimensions of LSTD and NDVI a <- raster("~/Documents/Cosas/MODIS_Fourier/Nuevas0.05/Europe/ParamTday_1.asc"); NAvalue(a) <- -9999 ### point to the directory for LSTD b <- raster("~/Documents/Cosas/MODIS_Fourier/Nuevas0.05/Europe/ParamNDVI_1.asc"); NAvalue(b) <- -9999 ### point to the directory for NDVI LSTD <- as.vector(a) NDVI <- as.vector(b) for (i in 1:nrow(Borrelia_ticks)) { if (file.exists(Borrelia_ticks[i,2])) { ######## Computes interaction of the vertebrate with the vector imagen_host <- raster(Borrelia_ticks[i,2]); NAvalue(imagen_host) <- -9999 imagen_tick <- raster(Borrelia_ticks[i,1]); NAvalue(imagen_tick) <- -9999 imagen_final <- (1+exp(2.297565*(sqrt(imagen_host/100)+sqrt(imagen_tick/100))))*Borrelia_ticks[i,3]*Borrelia_ticks[i,9] if (file.exists(Borrelia_other[i,2])) { ######## Computes interaction of the pathogen with the host temporal_nombres <- unlist(strsplit(Borrelia_ticks[i,2], split=".",fixed=TRUE))[1] imagen_final <- imagen_final*Borrelia_other[i,3]*Borrelia_other[i,5] Values_hosts <- as.vector(imagen_final) latabla <- as.data.frame(cbind(LSTD,NDVI,Values_hosts)) Values_hosts_2 <- rescale(Values_hosts, to=c(0,100)) latabla_2 <- as.data.frame(cbind(LSTD,NDVI,Values_hosts_2)) latabla <- na.omit(latabla) latabla_2 <- na.omit(latabla_2) latabla_plot <- subset(latabla_2, latabla_2[,3]>20) titulo <- paste(c(temporal_nombres),".png",sep="") titulo_legend <- c(temporal_nombres) binplot <- ggplot(latabla_plot, aes(x=LSTD, y=NDVI, z=latabla_plot[,3], colour=latabla_plot[,3])) + geom_point(aes(size=latabla_plot[,3],alpha=latabla_plot[,3]/100)) + scale_colour_gradientn(colours=rainbow(4),trans="reverse") + theme(axis.text=element_text(size=16), axis.title=element_text(size=18))+ ggtitle(titulo_legend) +theme(plot.title = element_text(size=18, face="italic"))+ theme(legend.title=element_blank()) ggsave(filename=titulo, plot=binplot) #### Saves the plot in the disk } } print(i*100/nrow(Borrelia_ticks)) }