# Analýza mnohostranných vyjednávání ve Světové obchodní organizaci
# Author   Michal Parízek
# URL 	http://mv.iir.cz/article/view/910
# Volume 	49
# Issue 	2
# Pages 	5-31
# Publication 	Mezinárodní vztahy

#This file gives the source code for the construction of the graphs used in the article.
#The file has two parts, one showing the construction of the radar charts (approximately lines 30-400), one showing the contour plots (approximately lines 400-600)


##################################################
#Install the necessary packages

#install.packages("fmsb")
library(fmsb)








############################################################
############################################################
#Radar charts




##################################################
#load dataset


###This file is NOT yet available for publication.
interests_dataset<-read.delim("interests_dataset.csv")


###############################################
#Preparations

#extract countries ideal points (columns 3 to 12 of the dataset)
positions<-interests_dataset[3:12]
#articulation_c_130827_fornash


#define the range of the radar charts for each of the 10 negotiation dimensions "direct..."
maxmin<-data.frame(
  directagr=c(10,0),
  directnam=c(10,0),
  directser=c(10,0),
  directsin=c(10,0),
  directdev=c(10,0),
  directldc=c(10,0),
  directipr=c(10,0),
  directenv=c(10,0),
  directsoc=c(10,0),
  directrul=c(10,0)
  )


#indicate the values for the Nash bargaining solutions (see table 1 on p. 22 in the article for the values)
#note that this value is a result of previous computation, here taken as an input for graph construction
nashes<-data.frame(
  directagr=c(5.59,6.66),
  directnam=c(7.21,6.97),
  directser=c(6.73,6.37),
  directsin=c(5.61,4.72),
  directdev=c(4.77,5.09),
  directldc=c(4.39,4.42),
  directipr=c(6.48,6.27),
  directenv=c(5.35,4.24),
  directsoc=c(4.57,2.81),
  directrul=c(7.03,6.95))



#extract from the main dataset a smaller dataframe with the columns for the salience each state attaches to all the areas
#note that the value "20" is arbitrary - it was chosen to make the data visible in the charts, it has no substantive meaning on its own
saliences<-data.frame(
  directagr=20*interests_dataset[13],
  directnam=20*interests_dataset[14],
  directser=20*interests_dataset[15],
  directsin=20*interests_dataset[16],
  directdev=20*interests_dataset[17],
  directldc=20*interests_dataset[18],
  directipr=20*interests_dataset[19],
  directenv=20*interests_dataset[20],
  directsoc=20*interests_dataset[21],
  directrul=20*interests_dataset[22])


#give the columns in the data frame "saliences" the right names
names(saliences) <- c("directagr","directnam","directser", "directsin","directdev","directldc","directipr","directenv","directsoc","directrul") 

#create a separate list with countries names to label the graphs later
names<-interests_dataset[1]



###############################################
#connect the individual dataframes into coherent dataframes with all the necessary information

#main descriptive graphs in two steps
radarchart_data<-rbind(maxmin,positions)
radarchart_data<-rbind(radarchart_data,saliences)

#chart with the Nash bargaining solution
radarchart_nash<-rbind(maxmin,nashes)




###############################################
#Radar chart with the positions of the Czech Republic and the EU (figure 1 on p. 17)


radarchart(radarchart_data[c(1,2,13,17), ], axistype=0,seg=5, plwd=3, 
           caxislabels=c("0", "", "", "","","10"),
           vlabels=c("Zemědělství (agr)","Průmysl (nam)","Služby (ser)","Singapurské otázky (sin)","Rozvoj (dev)","LDCs (ldc)"," Práve duševního vlastnictví (ipr)","Ochrana prostředí (env)","Sociální otázy (soc)","Pravidla WTO (rul)"),
           pcol=rgb(102,0,0,max = 255), plty=c(1,3), centerzero=TRUE)
           



###############################################
#Radar chart with the Nash bargaining solution (figure 5 on p. 23)


radarchart(radarchart_nash[c(1,2,3), ], axistype=0,seg=5, plwd=3, na.itp=FALSE,
           caxislabels=c("0", "", "", "","","10"),
           vlabels=c("Zemědělství (agr)","Průmysl (nam)","Služby (ser)","Singapurské otázky (sin)","Rozvoj (dev)","LDCs (ldc)"," Práva duševního vlastnictví (ipr)","Ochrana prostředí (env)","Sociální otázky (soc)","Pravidla WTO (rul)"),
           pcol=rgb(102,0,0,max = 255), plty=c(1,3), centerzero=TRUE)
           

###############################################
#Radar charts with positions and saliences of individual countries: ARE - IRL (figure on p. 30)

op <- par(mar=c(1, 1, 1, 1),mfrow=c(6, 4))

#ARE
radarchart(radarchart_data[c(1,2,3,50), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[1,])


radarchart(radarchart_data[c(1,2,4,51), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[2,])

radarchart(radarchart_data[c(1,2,5,52), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[3,])

radarchart(radarchart_data[c(1,2,6,53), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[4,])

radarchart(radarchart_data[c(1,2,7,54), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[5,])

radarchart(radarchart_data[c(1,2,8,55), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[6,])

radarchart(radarchart_data[c(1,2,9, 56), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[7,])

radarchart(radarchart_data[c(1,2,10, 57), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[8,])

radarchart(radarchart_data[c(1,2,11,58), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[9,])

radarchart(radarchart_data[c(1,2,12, 59), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[10,])

radarchart(radarchart_data[c(1,2,13, 60), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[11,])

radarchart(radarchart_data[c(1,2,14,61), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[12,])

radarchart(radarchart_data[c(1,2,15,62), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[13,])

radarchart(radarchart_data[c(1,2,16,63), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[14,])

radarchart(radarchart_data[c(1,2,17, 64), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[15,])

radarchart(radarchart_data[c(1,2,18,65), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[16,])

radarchart(radarchart_data[c(1,2,19,66), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[17,])

radarchart(radarchart_data[c(1,2,20,67), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[18,])

radarchart(radarchart_data[c(1,2,21,68), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[19,])

radarchart(radarchart_data[c(1,2,22,69), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[20,])
  
radarchart(radarchart_data[c(1,2,23,70), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[21,])

radarchart(radarchart_data[c(1,2,24,71), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[22,])

radarchart(radarchart_data[c(1,2,25,72), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[23,])

radarchart(radarchart_data[c(1,2,26,73), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"),pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[24,])



###############################################
#Radar charts with positions and saliences of individual countries: ISR - ZAF (figure on p. 31)

op <- par(mar=c(1, 1, 1, 1),mfrow=c(6, 4))



radarchart(radarchart_data[c(1,2,27,74), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[25,])

radarchart(radarchart_data[c(1,2,28,75), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[26,])

radarchart(radarchart_data[c(1,2,29,76), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[27,])

radarchart(radarchart_data[c(1,2,30,77), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[28,])

radarchart(radarchart_data[c(1,2,31,78), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[29,])

radarchart(radarchart_data[c(1,2,32,79), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[30,])

radarchart(radarchart_data[c(1,2,33,80), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[31,])

radarchart(radarchart_data[c(1,2,34,81), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[32,])
  
radarchart(radarchart_data[c(1,2,35,82), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[33,])

radarchart(radarchart_data[c(1,2,36,83), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[34,])

radarchart(radarchart_data[c(1,2,37,84), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[35,])

radarchart(radarchart_data[c(1,2,38,85), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[36,])

radarchart(radarchart_data[c(1,2,39,86), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[37,])

radarchart(radarchart_data[c(1,2,40,87), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[38,])

radarchart(radarchart_data[c(1,2,41,88), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[39,])

radarchart(radarchart_data[c(1,2,42,89), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[40,])

radarchart(radarchart_data[c(1,2,43,90), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[41,])

radarchart(radarchart_data[c(1,2,44,91), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[42,])

radarchart(radarchart_data[c(1,2,45,92), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[43,])

radarchart(radarchart_data[c(1,2,46,93), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[44,])

radarchart(radarchart_data[c(1,2,47,94), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[45,])

radarchart(radarchart_data[c(1,2,48,95), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[46,])

radarchart(radarchart_data[c(1,2,49,96), ], axistype=0,seg=5, plwd=3, na.itp=FALSE, 
           caxislabels=c("0", "", "", "","","10"), pcol=c(rgb(102,0,0,max = 255),rgb(204,0,0,max = 255)), pfcol=c(0,"darkred"), pdensity=c(0,100), plty=c(1,1), pty=32, centerzero=TRUE,
           vlabels=c("agr","nam","ser","sin","dev","ldc","ipr","env","soc","rul"),
           title=names[47,])








############################################################
############################################################
#Contour plots


#set to your working directory in the form dir="X:/folder/subfolder/"
dir="F:/BTS/mv_paper/"
setwd(dir)



###############################################
#Define the contour plot funciton
#Taken from:
#https://r-forge.r-project.org/scm/viewvc.php/pkg/R/ouch%20contour/Filled.contour2.R?view=annotate&root=phylooptim&sortby=rev


filled.contour2 <-
  function (x = seq(0, 1, length.out = nrow(z)),
            y = seq(0, 1, length.out = ncol(z)), z, xlim = range(x, finite = TRUE), 
            ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE), 
            levels = pretty(zlim, nlevels), nlevels = 20, color.palette = cm.colors, 
            col = color.palette(length(levels) - 1), plot.title, plot.axes, 
            key.title, key.axes, asp = NA, xaxs = "i", yaxs = "i", las = 1, 
            axes = TRUE, frame.plot = axes,mar, ...) 
  {
    # modification by Ian Taylor of the filled.contour function
    # to remove the key and facilitate overplotting with contour()
    if (missing(z)) {
      if (!missing(x)) {
        if (is.list(x)) {
          z <- x$z
          y <- x$y
          x <- x$x
        }
        else {
          z <- x
          x <- seq.int(0, 1, length.out = nrow(z))
        }
      }
      else stop("no 'z' matrix specified")
    }
    else if (is.list(x)) {
      y <- x$y
      x <- x$x
    }
    if (any(diff(x) <= 0) || any(diff(y) <= 0)) 
      stop("increasing 'x' and 'y' values expected")
    mar.orig <- (par.orig <- par(c("mar", "las", "mfrow")))$mar
    on.exit(par(par.orig))
    w <- (3 + mar.orig[2]) * par("csi") * 2.54
    par(las = las)
    mar <- mar.orig
    plot.new()
    par(mar=mar)
    plot.window(xlim, ylim, "", xaxs = xaxs, yaxs = yaxs, asp = asp)
    if (!is.matrix(z) || nrow(z) <= 1 || ncol(z) <= 1) 
      stop("no proper 'z' matrix specified")
    if (!is.double(z)) 
      storage.mode(z) <- "double"
    .Internal(filledcontour(as.double(x), as.double(y), z, as.double(levels), 
                            col = col))
    if (missing(plot.axes)) {
      if (axes) {
        title(main = "", xlab = "", ylab = "")
        Axis(x, side = 1)
        Axis(y, side = 2)
      }
    }
    else plot.axes
    if (frame.plot) 
      box()
    if (missing(plot.title)) 
      title(...)
    else plot.title
    invisible()
  }



###############################################
#Load all the datasets with the result of the calculations of states' support
#Read them directly from the file, or download them into your working directory
#Note that each of these datasets is itself a result of a simulation (done in Python)


#Figure 2, p. 19
#environmental and social issues
#reversion point at status quo (O) on all dimensions
zone_agreement_envsoc_rev0<-read.table("zone_agreement_envsoc_rev0.txt", header = FALSE, sep = ",", dec = ".")


#Figure 3, p. 20
#environmental and social issues
#reversion point at the least common denominator (O) on all dimensions
zone_agreement_envsoc_revmin<-read.table("zone_agreement_envsoc_revmin.txt", header = FALSE, sep = ",", dec = ".")

#Figure 4, p. 21
#agriculture and NAMA
#reversion point at value 3 on all dimensions
zone_agreement_agrnam_rev3<-read.table("zone_agreement_agrnam_rev3.txt", header = FALSE, sep = ",", dec = ".")



###############################################
#Some common preparations

#define colors with 14grades, high resolution on black and white
mycolors<-c(rgb(255,255,255,max = 255),rgb(255,229,255,max = 255),rgb(255,204,204,max = 255),rgb(255,153,153,max = 255),
            rgb(255,102,102,max = 255),rgb(255,51,51,max = 255),rgb(255,0,0,max = 255),
            rgb(204,0,0,max = 255),rgb(170,0,0,max = 255),rgb(140,0,0,max = 255),
            rgb(110,0,0,max = 255),rgb(80,0,0,max = 255),rgb(50,0,0,max = 255),rgb(0,0,0,max = 255))

#set the frame for each graph (margins, etc.)
op <- par(mar=c(2, 2, 2, 2),mfrow=c(1, 1))




###############################################
###Figure 2, p. 19

#preparation of the data matrix
mymatrixdata<-zone_agreement_envsoc_rev0
mymatrix<-data.matrix(mymatrixdata, rownames.force = NA)
z<-mymatrix
x <- 0:20
y <- 0:20


# define the levels of the graph (this has impact on the clarity of the picture and the size of the steps used)
mylevels <- seq(95,100,0.5) 

#the contour plot
filled.contour(x = seq(0, 10, length.out = nrow(z)),
               y = seq(0, 10, length.out = ncol(z)),
               z,
               ylab="Sociální otázky (soc)",
               xlab="Ochrana životního prostředí (env)", 
               cex.lab=1.5,
               color.palette=function(x)mycolors,
               #color.palette=function(x)rev(heat.colors(x)), #this is a sunny alternative
               levels=mylevels)




###############################################
###Figure 3, p. 20

#preparation of the data matrix
mymatrixdata<-zone_agreement_envsoc_revmin
mymatrix<-data.matrix(mymatrixdata, rownames.force = NA)
z<-mymatrix
x <- 0:20
y <- 0:20


# define the levels of the graph (this has impact on the clarity of the picture and the size of the steps used)
mylevels <- seq(90,100,1) 


#the contour plot
filled.contour(x = seq(0, 10, length.out = nrow(z)),
               y = seq(0, 10, length.out = ncol(z)),
               z,
               ylab="Sociální otázky (soc)",
               xlab="Ochrana životního prostředí (env)", 
               cex.lab=1.5,
               color.palette=function(x)mycolors,
               #color.palette=function(x)rev(heat.colors(x)), #this is a sunny alternative
               levels=mylevels)




###############################################
###Figure 4, p. 21

#preparation of the data matrix
mymatrixdata<-zone_agreement_agrnam_rev3
mymatrix<-data.matrix(mymatrixdata, rownames.force = NA)
z<-mymatrix
x <- 0:20
y <- 0:20

# define the levels of the graph (this has impact on the clarity of the picture and the size of the steps used)
mylevels <- seq(30,100,5) 

#the contour plot
filled.contour(x = seq(0, 10, length.out = nrow(z)),
               y = seq(0, 10, length.out = ncol(z)),
               z,
               color.palette=function(x)mycolors, #high resolution
               #color.palette=function(x)rev(heat.colors(x)), #this is a sunny alternative
               ylab="Průmyslové zboží (nam)",
               xlab="Zemědělství (agr)",cex.lab=1.5,
               levels=mylevels)