Last active
May 27, 2025 18:12
-
-
Save pablogarciadiaz/aa2c9b0eea7933c739f1eb4b22d04353 to your computer and use it in GitHub Desktop.
R script for creating, analysing, and visualising a conceptual system map of natural capital, ecosystem services, and ecosystem disservices
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ##### R SCRIPT FOR CREATING, ANALYSING, AND VISUALISING A CONCEPTUAL SYSTEM MAP OF NATURAL CAPITAL, | |
| ##### ECOSYSTEM SERVICES AND DISSERVICES USING NETWORK ANALYSES. | |
| ### CASE STUDY: RIO LULES WATERSHED, TUCUMAN, ARGENTINA | |
| #### This script should be used alongside the CSV files 'Variables-CompleteWithDescriptions.csv' and | |
| #### 'LinksbtwDrivers-ForiGraph.csv' | |
| #### Pablo Garcia Diaz, Instituto de Ecologia Regional (UNT-CONICET, Tucuman, Argentina), 2025 | |
| #### E-mail: pablo.garciadiaz.eco@gmail.com | |
| set.seed(2021) | |
| #### Load packages | |
| library(igraph) | |
| library(reshape2) | |
| library(ggplot2) | |
| library(raster) | |
| library(network) | |
| library(RColorBrewer) | |
| library(cowplot) | |
| library(ggrepel) | |
| #### Read data | |
| ### Read the variables | |
| variables<-read.csv("./Variables-CompleteWithDescriptions.csv", header=TRUE, sep=",") | |
| ### Select the columns with the name and type of the drivers | |
| variables<-variables[ , c(1, 2)] | |
| ### Change column names | |
| colnames(variables)<-c("id", "Tipo") | |
| ### Assign membership to each driver (type) | |
| variables$membership<-as.numeric(as.factor((variables$Tipo))) | |
| head(variables) | |
| summary(variables) | |
| ### Read the csv with the links between drivers | |
| links<-read.csv("./LinksbtwDrivers-ForiGraph.csv", header=TRUE, sep=",") | |
| head(links) | |
| summary(links) | |
| #### Convert to a network using iGraph | |
| net.tot<-graph_from_data_frame(d=links, vertices=variables, directed=T) | |
| #### Calculate betweenness (remember that this is a directed graph) | |
| btw<-betweenness(net.tot, directed=T, weights=NA) | |
| ### Calculate out degree (number of outgoing connections from each driver) | |
| degree.out<-degree(net.tot, mode="out") | |
| ##### Simple paths from each driver to each of the EDS and total number of paths to all EDS | |
| ### Number of drivers | |
| n.vars<-length(V(net.tot)) | |
| ### Empty objects to store the results | |
| paths<-paths.leish<-paths.dengue<-paths.carbon<-paths.quality<-paths.quantity<-rep(NA, n.vars) | |
| ### Loop to calculate the simple paths | |
| for(i in 1:n.vars){ | |
| #### Total number of simple paths from each driver to each EDS | |
| paths[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Leishmaniasis"))+length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Aboveground_Carbon"))+length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Water_quality"))+length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Water_quantity"))+length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Dengue")) | |
| #### Simple paths from each driver to leishmaniasis transmission | |
| paths.leish[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Leishmaniasis")) | |
| #### Simple paths from each driver to aboveground carbon | |
| paths.carbon[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Aboveground_Carbon")) | |
| #### Simple paths from each driver to water quality | |
| paths.quality[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Water_quality")) | |
| #### Simple paths from each driver to water quantity | |
| paths.quantity[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Water_quantity")) | |
| #### Simple paths from each driver to dengue transmission | |
| paths.dengue[i]<-length(all_simple_paths(net.tot, from=V(net.tot)[i], to="Dengue")) | |
| } | |
| paths | |
| #### Join all results in a single dataframe | |
| summ.graph<-data.frame(Variable=names(btw), betweenness=btw, degree_out=degree.out, tot.paths=paths, | |
| carbon.paths=paths.carbon, leish.paths=paths.leish, dengue.paths=paths.dengue, quality.paths=paths.quality, | |
| quantity.paths=paths.quantity) | |
| summ.graph | |
| ### Bivariate scatterplots showing the relationships between different network metrics | |
| variables2<-variables | |
| colnames(variables2)<-c("Variable", "Tipo", "membership") | |
| head(variables2) | |
| #### Get the data together | |
| scatter.for.plot<-merge(summ.graph, variables2, by=c("Variable"), all=TRUE) | |
| scatter.for.plot$Tipo<-as.factor(scatter.for.plot$Tipo) | |
| summary(scatter.for.plot) | |
| #### Spearman correlations | |
| #### Betweenness and out-degree | |
| cor(scatter.for.plot$betweenness, scatter.for.plot$degree_out, method="spearman") | |
| #### Betweenness and log(total simple paths) | |
| cor(scatter.for.plot$betweenness, log(scatter.for.plot$tot.paths+1), method="spearman") | |
| ### Out degree and log(total simple paths) | |
| cor(scatter.for.plot$degree_out, log(scatter.for.plot$tot.paths+1), method="spearman") | |
| #### Scatterplots | |
| ### Define colours for each type of driver | |
| scatter.for.plot$col.points<-rainbow(5, alpha=1)[scatter.for.plot$membership] | |
| leg.col<-data.frame(Tipo=as.factor(scatter.for.plot$Tipo), Colour=scatter.for.plot$col.points) | |
| col.leg<-unique(leg.col) | |
| col.leg1<-data.frame(rbind(col.leg[1, ], leg.col[4,], leg.col[5, ], leg.col[2, ], leg.col[3, ])) | |
| colnames(col.leg1)<-c("Tipo", "Colour") | |
| col.leg1 | |
| ### Plotting including exporting in png format | |
| png('./CorrelationsBtw2.png', | |
| width=2500*3+100, height=2500, res=300, pointsize=60) | |
| ### Betweenness and out degree | |
| p1<-ggplot(scatter.for.plot, aes(x=degree_out, y=betweenness, color=factor(Tipo))) + | |
| geom_point(size=5, shape=21, aes(fill=factor(Tipo))) + | |
| scale_fill_manual(values=scatter.for.plot$col.points, breaks=c('ESD', 'Abiotic', 'Biotic', 'Policies', 'Socio-economic'))+ | |
| xlab("Out degree") + | |
| ylab("Betweenness") + | |
| theme_bw() + | |
| theme(text = element_text(size=15)) + | |
| theme(legend.position="none") + | |
| #### Add labels to outliers: Human population; deforestation and fragmentation; elevation; vegetation cover; temperature | |
| geom_label_repel( | |
| aes(label = Variable), | |
| data = scatter.for.plot[c(8, 10, 14, 35, 41), ], | |
| color = "black", fill= NA, box.padding = 1.5) | |
| #### Betweenness and log(total simple paths) | |
| p2<-ggplot(scatter.for.plot, aes(x=log(tot.paths+1), y=betweenness, color=factor(Tipo))) + | |
| geom_point(size=5, shape=21, aes(fill=factor(Tipo))) + | |
| scale_fill_manual(values=scatter.for.plot$col.points, breaks=c('ESD', 'Abiotic', 'Biotic', 'Policies', 'Socio-economic'))+ | |
| xlab("Simple paths (log[x+1] transformed)") + | |
| ylab("Betweenness") + | |
| theme_bw() + | |
| theme(text = element_text(size=15)) + | |
| theme(legend.position="none") + | |
| #### Add labels to outliers: Human population; deforestation and fragmentation; elevation; vegetation cover; temperature | |
| geom_label_repel( | |
| aes(label = Variable), | |
| data = scatter.for.plot[c(8, 10, 14, 35, 41), ], | |
| color = "black", fill= NA, box.padding = 1.5) | |
| ### Out degree and log(total simple paths) | |
| p3<-ggplot(scatter.for.plot, aes(x=log(tot.paths+1), y=degree_out, color=factor(Tipo))) + | |
| geom_point(size=5, shape=21, aes(fill=factor(Tipo))) + | |
| scale_fill_manual(values=scatter.for.plot$col.points, breaks=c('ESD', 'Abiotic', 'Biotic', 'Policies', 'Socio-economic'))+ | |
| xlab("Simple paths (log[x+1] transformed)") + | |
| ylab("Out degree") + | |
| theme_bw() + | |
| theme(text = element_text(size=15)) + | |
| #### Add legend | |
| theme(legend.position="right", legend.title=element_blank(), legend.text = element_text(size = 20)) + | |
| #### Add labels to outliers: Human population; deforestation and fragmentation; elevation; vegetation cover; temperature | |
| geom_label_repel( | |
| aes(label = Variable), | |
| data = scatter.for.plot[c(8, 10, 14, 35, 41), ], | |
| color = "black", fill= NA, box.padding = 1.5) | |
| ### Remove multiple legends | |
| p4<-p3+guides(color = FALSE, size = FALSE) | |
| ##### Plot the three graphs in a grid | |
| plot_grid(p1, p2,p4, labels = c('A', 'B', 'C'), ncol=3, rel_widths = c(1, 1, 1.3), label_size = 15) | |
| dev.off() | |
| ##### Plotting the networks and exporting them in png format | |
| #### Define colours | |
| V(net.tot)$group<-variables$membership | |
| vertex.col<-rainbow(5, alpha=0.4)[variables$membership] | |
| uni.cols<-unique(cbind(vertex.col, variables$membership )) | |
| uni.cols | |
| uni.cols<-uni.cols[order(uni.cols[, 2]), ] | |
| ### Define layout for the networks | |
| LO1 = layout_with_fr(net.tot) | |
| #### Network with all vertices of the same size | |
| png('./iGraphTot-NoSize.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| size.node<-rep(8, n.vars) | |
| #### EDS assigned larger sizes | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=15 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=15 | |
| size.node[which(names(V(net.tot))=="Dengue")]=15 | |
| ### Change vertex shape of EDS to square | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| ### Plotting | |
| plot(net.tot,layout=LO1, margin=0, | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| #### Adding the legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ##### Network in which the size of the vertices is proportional to the total number of simple paths | |
| png('./iGraphTot-SizePath.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| ### Define the scaling according to the number of simple paths | |
| size.node<-log10(paths.carbon+1)*5 | |
| #### EDS vertices with the same size | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=15 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=15 | |
| size.node[which(names(V(net.tot))=="Dengue")]=15 | |
| ### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| #### Plotting | |
| plot(net.tot,layout=LO1, margin=0, | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| #### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ################ Visualising the networks based on the number of single paths to each EDS | |
| ##### Size of the vertices proportional to the number of single paths to aboveground carbon | |
| png('./iGraphTot-SizePathCarbono.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| #### Scaling | |
| size.node<-log10(paths.carbon+1)*5 | |
| #### The target EDS is assigned a larger size than the other four EDS | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=2 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=2 | |
| size.node[which(names(V(net.tot))=="Dengue")]=2 | |
| ### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| #### Plotting | |
| plot(net.tot,layout=LO1, margin=0, | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| title("Aboveground Carbon", cex.main=3) | |
| #### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ##### Size of the vertices proportional to the number of single paths to leishmaniasis transmision | |
| png('./iGraphTot-SizePathLeishmaniasis.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| ### Scaling | |
| size.node<-log10(paths.leish+1)*5 | |
| #### The target EDS is assigned a larger size than the other four EDS | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=15 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=2 | |
| size.node[which(names(V(net.tot))=="Dengue")]=2 | |
| #### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| ### Plotting | |
| plot(net.tot,layout=LO1, margin=0, main="", | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| title("Leishmaniasis transmission", cex.main=3) | |
| #### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ##### Size of the vertices proportional to the number of single paths to dengue transmision | |
| png('./iGraphTot-SizePathDengue.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| ### Scaling | |
| size.node<-log10(paths.dengue+1)*5 | |
| #### The target EDS is assigned a larger size than the other four EDS | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=2 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=2 | |
| size.node[which(names(V(net.tot))=="Dengue")]=15 | |
| #### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| ### Plotting | |
| plot(net.tot,layout=LO1, margin=0, | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| title("Dengue transmission", cex.main=3) | |
| ### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ##### Size of the vertices proportional to the number of single paths to water quality | |
| png('./iGraphTot-SizePathWaterQuality.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| ### Scaling | |
| size.node<-log10(paths.quality+1)*5 | |
| #### The target EDS is assigned a larger size than the other four EDS | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=2 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=15 | |
| size.node[which(names(V(net.tot))=="Dengue")]=2 | |
| #### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| #### Plotting | |
| plot(net.tot,layout=LO1, margin=0, | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| title("Water Quality", cex.main=3) | |
| ### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| ##### Size of the vertices proportional to the number of single paths to water quantity | |
| png('./iGraphTot-SizePathWaterQuantity.png', | |
| width=6000, height=6000, res=300, pointsize=15) | |
| #### Scaling | |
| size.node<-log10(paths.quantity+1)*5 | |
| #### The target EDS is assigned a larger size than the other four EDS | |
| size.node[which(names(V(net.tot))=="Leishmaniasis")]=2 | |
| size.node[which(names(V(net.tot))=="Aboveground_Carbon")]=2 | |
| size.node[which(names(V(net.tot))=="Water_quantity")]=15 | |
| size.node[which(names(V(net.tot))=="Water_quality")]=2 | |
| size.node[which(names(V(net.tot))=="Dengue")]=2 | |
| #### Change shape of EDS vertices | |
| vertex.shape<-rep("circle", n.vars) | |
| vertex.shape[which(names(V(net.tot))=="Leishmaniasis")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Aboveground_Carbon")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quantity")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Water_quality")]="square" | |
| vertex.shape[which(names(V(net.tot))=="Dengue")]="square" | |
| V(net.tot)$label.cex<-2 | |
| #### Plotting | |
| plot(net.tot,layout=LO1, margin=0, main="", | |
| vertex.color=vertex.col, vertex.label=variables$id, vertex.shape=vertex.shape, vertex.size=size.node) | |
| title("Water Quantity", cex.main=3) | |
| #### Add legend | |
| data.legend1<-cbind(levels(as.factor(variables$Tipo)), uni.cols[,1]) | |
| data.legend<-rbind(data.legend1[3,], data.legend1[1, ], data.legend1[2, ], data.legend1[4, ], data.legend1[5, ]) | |
| legend(x=-1.1, y=1.1,bty = "n", cex=2, | |
| legend=data.legend[,1], | |
| fill=data.legend[, 2], border=NA) | |
| dev.off() | |
| #### END OF SCRIPT #### |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment