#libraries
library(plyr)
library(dplyr)
library(pbapply)
library(igraph)
library(openxlsx)
library(writexl)
#####Word cloud############
library(ggplot2)
library(ggwordcloud)
####point diagram##############
library(ggplot2)
library(ggrepel)
#####Maps############
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggrepel) 

#xyz t?bla el??ll?t:
#xyz- to word cloud
#xyz.kw_keywords
#xyz.c:clusters number; freq:size, szavak freq; kw:a kulcsszavak, hex:hexadecimalis szink?dok szavank?nt klaszterenkent
#xyz.hex:hexadecimal colours

# kw frec by clusters #w.cloud1.clus1$c=1 
#kw.p.clus.1 cl6 to cl2 sust, tech clust;  in kw.p.clus.1 cl1: sustainability=cl1, cl2-6:technological oriented clusters=cl2
kw.p.cl2.p1=kw.p.clus.1
kw.p.cl2.p1$clus=ifelse(kw.p.clus.1$clus>1,2,1)

kw.p.cl2.p2=kw.p.clus.2
kw.p.cl2.p2$clus=ifelse(kw.p.clus.2$clus==2,1,2) #cl1 sust, cl2:tech

#w.cloud1.cl1p1 from kw.countKL
w.cl1p1=subset(kw.countKL,clus==1)
kw.countKL.to2clusp1=kw.countKL
kw.countKL.to2clusp1$clus=ifelse(kw.countKL$clus==1,1,2)
p=subset(kw.countKL.to2clusp1,clus==2)
#aggregate by kw
p2=plyr::count(p, vars=c("kw"))
w.cl2p1=p2
w.cl2p1$clus=2
rm(p,p2)
#Period2
w.cl1p2=subset(kw.countKL2,clus==2)
w.cl1p2$clus=1
kw.countKL.to2clus=kw.countKL2
kw.countKL.to2clus$clus=ifelse(kw.countKL2$clus==2,1,2)
p=subset(kw.countKL.to2clus,clus==2)
#aggregate by kw
p2=plyr::count(p, vars=c("kw"))
w.cl2p2=p2
w.cl2p2$clus=2

w.cl1p1=w.cl1p1[,c(2,3,1)]
w.cl1p2=w.cl1p2[,c(2,3,1)]

c=colnames(w.cloud1.clus1)
colnames(w.cl1p1)=c
colnames(w.cl2p1)=c
colnames(w.cl1p2)=c
colnames(w.cl2p2)=c

#FIGUREA3A, FigureA3b-----------------------
#word clouds for 2 aggregated cluster for 2 periods
#word cloud smart city 2 periods

#colors and tables to word clouds
xyz.p1=do.call(rbind,list(w.cl1p1,w.cl2p1))

c.minmeanmax=lapply(c(1:2),function(x){
  a=subset(xyz.p1,clus==x)
  min=min(a$freq)
  mean=mean(a$freq)
  max=max(a$freq)
  b=as.data.frame(min,optional=TRUE,stringAsFactots=F)
  colnames(b)="min"
  b$mean=mean    
  b$max=max
  b$c=x
  return(b)})
p1.minmeanmax=do.call(rbind,c.minmeanmax)
# hexa color codes add
t1=c("#C9C7C3","#16CB97","#038D67",1)
t=c("#C9C7C3","#FFA07A","#CD8162",2)

h.col.p1=as.data.frame(do.call(rbind,list(t1,t)))
colnames(h.col.p1)=c("hmin" , "hmean", "hmax" , "c" )
library(pbapply)

xyz.h=pblapply(c(1:2),function(x){
  a=subset(xyz.p1,clus==x)
  #COLOR scale with pal function
  pal <- scales::gradient_n_pal(colours = c(h.col.p1[x,1],h.col.p1[x,2],h.col.p1[x,3]),
                                
                                values= c(p1.minmeanmax[x,1],p1.minmeanmax[x,2],p1.minmeanmax[x,3]))
  
  a$hex <- pal(a$freq)
  return(a)})
xyz.p1=do.call(rbind,xyz.h)

library(plyr)
xyz.p1=arrange(xyz.p1,desc(freq))

wc.p1 <- ggplot(transform(xyz.p1, c=factor(clus,levels=c("1","2"))),
               aes(label = kw, size = freq, color=hex)) +
  geom_text_wordcloud_area(shape="circle", eccentricity = 1) +
  scale_size_area(max_size = 4) +
  scale_size(range = c(2,12))+
  theme_minimal() +
  theme(strip.text.x = element_blank())+
  scale_color_identity()+
  #scale_fill_manual(values = xyz$hex)+
  #scale_color_gradientn(colors=brewer.pal(9,"Reds"))+
  facet_wrap(~c, ncol=3)

ggsave("wc.cl2p1.png", wc.p1, height = 8, width = 12, dpi = 240)

#2-es period
xyz.p2=do.call(rbind,list(w.cl1p2,w.cl2p2))

c.minmeanmax=lapply(c(1:2),function(x){
  a=subset(xyz.p2,clus==x)
  min=min(a$freq)
  mean=mean(a$freq)
  max=max(a$freq)
  b=as.data.frame(min,optional=TRUE,stringAsFactots=F)
  colnames(b)="min"
  b$mean=mean    
  b$max=max
  b$c=x
  return(b)})

p2.minmeanmax=do.call(rbind,c.minmeanmax)

library(pbapply)

xyz.h=pblapply(c(1:2),function(x){
  a=subset(xyz.p2,clus==x)
  #fv def az ?tmenetes hexadec. szinekhez  #hexadec szinek pal fv-el
  pal <- scales::gradient_n_pal(colours = c(h.col.p1[x,1],h.col.p1[x,2],h.col.p1[x,3]),
                                
                                values= c(p2.minmeanmax[x,1],p2.minmeanmax[x,2],p2.minmeanmax[x,3]))
  
  a$hex <- pal(a$freq)
  return(a)})
xyz.p2=do.call(rbind,xyz.h)

library(plyr)
xyz.p2=arrange(xyz.p2,desc(freq))

wc.p2 <- ggplot(transform(xyz.p2, c=factor(clus,levels=c("1","2"))),
                aes(label = kw, size = freq, color=hex)) +
  geom_text_wordcloud_area(shape="circle", eccentricity = 1) +
  scale_size_area(max_size = 4) +
  scale_size(range = c(2,12))+
  theme_minimal() +
  theme(strip.text.x = element_blank())+
  scale_color_identity()+
  #scale_fill_manual(values = xyz$hex)+
  #scale_color_gradientn(colors=brewer.pal(9,"Reds"))+
  facet_wrap(~c, ncol=3)

ggsave("wc.cl2p2.png", wc.p2, height = 8, width = 12, dpi = 240)


#-------------------------------------------
#inter, intra relations Macro regio level
#macro reg codes: Af:Africa, NAm: U.S., Canada, LAm: Latin American countries, EU:Europe, As: Asia, Pac:Pacific and Australia
#c1.cityid.mregio: UT, City coordinates, macroreg kod
#kw.p.cl2.p1: UT, Clus id; macroreg categories. in period1
#kw.p.cl2.p2: UT, Clus id; macroreg categories. in period2
#macro.regio add to UT
macro.reg=c1.cityid.mregio[,c("UT","macro.reg")]

p=merge(kw.p.cl2.p1,macro.reg,by.x="UT",by.y="UT",all.x=T,all.y=F,sort=F) #pubs which have inter macro regional collaborations appear with all macro regio name, so there are more records in this table
kw.p.cl2.p1.macroreg=p
rm(p)

p=merge(kw.p.cl2.p2,macro.reg,by.x="UT",by.y="UT",all.x=T,all.y=F,sort=F) #pubs which have inter macro regional collaborations appear with all macro regio name, so there are more records in this table
kw.p.cl2.p2.macroreg=p
rm(p)

p=merge(kw.p.clus.1,macro.reg,by.x="name",by.y="UT",all.x=T,all.y=F,sort=F) #pubs which have inter macro regional collaborations appear with all macro regio name, so there are more records in this table
kw.p.clus.1.macroreg=p
rm(p)

p=merge(kw.p.clus.2,macro.reg,by.x="name",by.y="UT",all.x=T,all.y=F,sort=F) #pubs wich have inter macro regional collaborations appear with all macro regio name, so there are more records in this table
kw.p.clus.2.macroreg=p
rm(p)
kw.p.cl2.p1.macroreg=unique(kw.p.cl2.p1.macroreg)
kw.p.cl2.p2.macroreg=unique(kw.p.cl2.p2.macroreg)
kw.p.clus.1.macroreg=unique(kw.p.clus.1.macroreg)
kw.p.clus.2.macroreg=unique(kw.p.clus.2.macroreg)
unique(kw.p.cl2.p1.macroreg$macro.reg)

#FOR Figure02; Figure04; Figure05; Figure06; based on Mathiensen et al., 2010, Figure6

#CNCI from incites: https://access.clarivate.com/login?app=incites
  #1WoS pubs find by UT-s
  #2 Export from WoS to Incites;
  #3. sign in Incites and dowload CNCI, top10, %top5% by papers
  #4. read the xlsx
#p1.imp, p2.imp, instead of UTs with randomized ids of the publication and std CNCI

#------------------------------------------------------------------------------------

 ##city number in both periods all SC fields 3-5 collab cat.
 p1.cl2=subset(kw.p.cl2.p1,regional.collab.category>2)
 p2.cl2=subset(kw.p.cl2.p2,regional.collab.category>2)
 #cluster and reg.collab.categories added impact values

 #p1.imp, p2.imp, instead of UTs with randomized ids of the publication and std CNCI
 p1.cl2=merge(p1.cl2,p1.imp,by.x="UT",by.y="UT",all=F)
 p2.cl2=merge(p2.cl2,p2.imp,by.x="UT",by.y="UT",all=F) #in Incites 4 pubs not fiend with CNCI, these are erised
 
 #AFFIL geography info added
 
p1.cl.aff=merge(p1.cl2,c1.cityid.mregio,by.x="UT",by.y="UT",all.x=T,all.y=F) #most csak a intercity collab. 
          
rm(p2.cl.aff)                   
p2.cl.aff=merge(p2.cl2,c1.cityid.mregio,by.x="UT",by.y="UT",all.x=T,all.y=F) #most csak a intercity collab. 
       


nc1=unique(p1.cl.aff$city.id) #1282
nc2=unique(p2.cl.aff$city.id) #2635
nc=intersect(nc1,nc2) #1058 city)

# 
p1.city=subset(p1.cl.aff, city.id %in% nc)
p2.city=subset(p2.cl.aff, city.id %in% nc)


#output for 1058 cities
p1.output=plyr::count(p1.city, vars=c("city.id"))
colnames(p1.output)=c("city.id","output")

cnci=aggregate(p1.city$CNCI, by=list(p1.city$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI")
p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
p1.output$mean.CNCI=p1.output$szum.CNCI/p1.output$output
colnames(p1.output)=c("city.id",   "output.p1",    "szum.CNCI.p1", "mean.CNCI.p1")

p2.output=plyr::count(p2.city, vars=c("city.id"))
colnames(p2.output)=c("city.id","output")

cnci=aggregate(p2.city$CNCI, by=list(p2.city$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI")
p2.output=merge(p2.output,cnci,by.x="city.id",by.y="city.id",all=T)
p2.output$mean.CNCI=p2.output$szum.CNCI/p2.output$output
colnames(p2.output)=c("city.id",   "output.p2",    "szum.CNCI.p2", "mean.CNCI.p2")

imp.output=merge(p1.output,p2.output,by.x="city.id",by.y="city.id",all=T)
#cl1 and cl2 calculate output, szum.CNCI,mean.CNCI by two periods 2010-2017, 2018-2022


#output for CL1 (Sustainability cluster) 1058 cities
p1.city.cl1=subset(p1.city,clus==1)
p1.city.cl2=subset(p1.city,clus==2)
p2.city.cl1=subset(p2.city,clus==1)
p2.city.cl2=subset(p2.city,clus==2)

#period1, period2 cl1,cl2
#p1 cl1
p1.output=plyr::count(p1.city.cl1, vars=c("city.id"))
colnames(p1.output)=c("city.id","output.cl1")

cnci=aggregate(p1.city.cl1$CNCI, by=list(p1.city.cl1$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI.cl1")
p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
p1.output$mean.CNCI.cl1=p1.output$szum.CNCI.cl1/p1.output$output.cl1
colnames(p1.output)=c("city.id",   "output.p1.cl1",    "szum.CNCI.p1.cl1", "mean.CNCI.p1.cl1")
p1.output.cl1=p1.output
 #p1 cl2
p1.output=plyr::count(p1.city.cl2, vars=c("city.id"))
colnames(p1.output)=c("city.id","output.cl2")

cnci=aggregate(p1.city.cl2$CNCI, by=list(p1.city.cl2$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI.cl2")
p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
p1.output$mean.CNCI.cl2=p1.output$szum.CNCI.cl2/p1.output$output.cl2
colnames(p1.output)=c("city.id",   "output.p1.cl2",    "szum.CNCI.p1.cl2", "mean.CNCI.p1.cl2")
p1.output.cl2=p1.output

#p2 cl1
p1.output=plyr::count(p2.city.cl1, vars=c("city.id"))
colnames(p1.output)=c("city.id","output.cl1")

cnci=aggregate(p2.city.cl1$CNCI, by=list(p2.city.cl1$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI.cl1")
p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
p1.output$mean.CNCI.cl1=p1.output$szum.CNCI.cl1/p1.output$output.cl1
colnames(p1.output)=c("city.id",   "output.p2.cl1",    "szum.CNCI.p2.cl1", "mean.CNCI.p2.cl1")
p2.output.cl1=p1.output
#p2 cl2
p1.output=plyr::count(p2.city.cl2, vars=c("city.id"))
colnames(p1.output)=c("city.id","output.cl2")

cnci=aggregate(p2.city.cl2$CNCI, by=list(p2.city.cl2$city.id), FUN=sum)
colnames(cnci)=c("city.id","szum.CNCI.cl2")
p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
p1.output$mean.CNCI.cl2=p1.output$szum.CNCI.cl2/p1.output$output.cl2
colnames(p1.output)=c("city.id",   "output.p2.cl2",    "szum.CNCI.p2.cl2", "mean.CNCI.p2.cl2")
p2.output.cl2=p1.output
imp.output=imp.output[,c(1:7)]
imp.output=merge(imp.output,p1.output.cl1,by.x="city.id",by.y="city.id",all=T)
imp.output=merge(imp.output,p2.output.cl1,by.x="city.id",by.y="city.id",all=T)
imp.output=merge(imp.output,p1.output.cl2,by.x="city.id",by.y="city.id",all=T)
imp.output=merge(imp.output,p2.output.cl2,by.x="city.id",by.y="city.id",all=T)

nc.cls=imp.output[,c("city.id","output.p1.cl1","output.p1.cl2","output.p2.cl1","output.p2.cl2")]
nc.cls$all=paste(nc.cls$output.p1.cl1,nc.cls$output.p1.cl2,nc.cls$output.p2.cl1,nc.cls$output.p2.cl2,sep="")
nc.cls.all=grep("\\.*NA",nc.cls$all)
nc.cls2=imp.output
nc.cls2$all[nc.cls.all]=0
nc.cls2$all=ifelse(is.na(nc.cls2$all)==T,1,0)
#subset imp.out remain those cities if the city appear both periods and both cluster, so the city have a diverse reserch output
imp.output.r=subset(nc.cls2,all==1) #346 city 

#categorise 346 cities based on Matthiensen, 2010 Figure 6; unless for citation change we use normalize citation score namely WoS CNCI value, and calculate the mean CNCI for each city in each periods
#subset based on output 5 papers/cluster/periods
imp.output.49=subset(imp.output.r,output.p1.cl1>4)
imp.output.49=subset(imp.output.49,output.p2.cl1>4)
imp.output.49=subset(imp.output.49,output.p1.cl2>4)
imp.output.49=subset(imp.output.49,output.p2.cl2>4)

#Categorised 49 cities by Mathiensen methods
#Change of output the whole SC field
#percent growth table
c49.growth=as.data.frame(imp.output.49[,1])
colnames(c49.growth)="city.id"
c49.growth$ch.output=(imp.output.49$output.p2/imp.output.49$output.p1)*100
c49.growth$ch.meanCNCI=(imp.output.49$mean.CNCI.p2/imp.output.49$mean.CNCI.p1)*100

#Mean output, Mean meanCNCI
M49.output=mean(c49.growth$ch.output)
M49.output.k25=M49.output*0.75  #253.24
M49.output.n25=M49.output*1.25  #422.06

M49.CNCI=mean(c49.growth$ch.meanCNCI)
M49.CNCI.k25=M49.CNCI*0.75   #76.54
M49.CNCI.n25=M49.CNCI*1.25   #127.56

#categorise
c49.growth$cat.output=ifelse(c49.growth$ch.output>M49.output.n25,"i","n")
c49.growth$cat.output=ifelse(c49.growth$ch.output<M49.output.k25,"d",c49.growth$cat.output)
c49.growth$cat.CNCI=ifelse(c49.growth$ch.meanCNCI>M49.CNCI.n25,"i","n")
c49.growth$cat.CNCI=ifelse(c49.growth$ch.meanCNCI<M49.CNCI.k25,"d",c49.growth$cat.CNCI)

c49.growth$cat=paste(c49.growth$cat.output,c49.growth$cat.CNCI,sep=";")
c49.growth$cat=ifelse(c49.growth$cat=="n;n","n","")

c49.growth$cat2.output=ifelse(c49.growth$ch.output>M49.output,"i","d")
c49.growth$cat2.CNCI=ifelse(c49.growth$ch.meanCNCI>M49.CNCI,"i","d")
c49.growth$cat2=paste(c49.growth$cat2.output,c49.growth$cat2.CNCI,sep=";")
c49.growth$cat3=ifelse(c49.growth$cat=="n","n",c49.growth$cat2)
c49.growth$cat=c49.growth$cat3
c49.growth=c49.growth[,-c(4,5,7:10)]

#cl1 and cl2-re ugyanez a 49 cities-re
c49.growth$ch.output.cl1=(imp.output.49$output.p2.cl1/imp.output.49$output.p1.cl1)*100
c49.growth$ch.meanCNCI.cl1=(imp.output.49$mean.CNCI.p2.cl1/imp.output.49$mean.CNCI.p1.cl1)*100

#Mean output, Mean meanCNCI
M49.outputcl1=mean(c49.growth$ch.output.cl1)
M49.outputcl1.k25=M49.outputcl1*0.75  #315.36
M49.outputcl1.n25=M49.outputcl1*1.25  #525.60

M49.CNCIcl1=mean(c49.growth$ch.meanCNCI.cl1)
M49.CNCIcl1.k25=M49.CNCIcl1*0.75   #100.72
M49.CNCIcl1.n25=M49.CNCIcl1*1.25   #167.87

#categorise
c49.growth$cat.output.cl1=ifelse(c49.growth$ch.output.cl1>M49.outputcl1.n25,"i","n")
c49.growth$cat.output.cl1=ifelse(c49.growth$ch.output.cl1<M49.outputcl1.k25,"d",c49.growth$cat.output.cl1)
c49.growth$cat.CNCI.cl1=ifelse(c49.growth$ch.meanCNCI.cl1>M49.CNCIcl1.n25,"i","n")
c49.growth$cat.CNCI.cl1=ifelse(c49.growth$ch.meanCNCI.cl1<M49.CNCIcl1.k25,"d",c49.growth$cat.CNCI.cl1)

c49.growth$cat.cl1=paste(c49.growth$cat.output.cl1,c49.growth$cat.CNCI.cl1,sep=";")
c49.growth$cat.cl1=ifelse(c49.growth$cat.cl1=="n;n","n","")

c49.growth$cat2.output.cl1=ifelse(c49.growth$ch.output.cl1>M49.outputcl1,"i","d")
c49.growth$cat2.CNCI.cl1=ifelse(c49.growth$ch.meanCNCI.cl1>M49.CNCIcl1,"i","d")
c49.growth$cat2.cl1=paste(c49.growth$cat2.output.cl1,c49.growth$cat2.CNCI.cl1,sep=";")
c49.growth$cat3.cl1=ifelse(c49.growth$cat.cl1=="n","n",c49.growth$cat2.cl1)
c49.growth$cat.cl1=c49.growth$cat3.cl1
c49.growth=c49.growth[,-c(7,8,10:13)]
 
#--------------
#cl2
c49.growth$ch.output.cl2=(imp.output.49$output.p2.cl2/imp.output.49$output.p1.cl2)*100
c49.growth$ch.meanCNCI.cl2=(imp.output.49$mean.CNCI.p2.cl2/imp.output.49$mean.CNCI.p1.cl2)*100

#Mean output, Mean meanCNCI
M49.outputcl2=mean(c49.growth$ch.output.cl2)
M49.outputcl2.k25=M49.outputcl2*0.75  #243.23
M49.outputcl2.n25=M49.outputcl2*1.25  #405.38

M49.CNCIcl2=mean(c49.growth$ch.meanCNCI.cl2)
M49.CNCIcl2.k25=M49.CNCIcl2*0.75   #84.52
M49.CNCIcl2.n25=M49.CNCIcl2*1.25   #140.86

#categorise
c49.growth$cat.output.cl2=ifelse(c49.growth$ch.output.cl2>M49.outputcl2.n25,"i","n")
c49.growth$cat.output.cl2=ifelse(c49.growth$ch.output.cl2<M49.outputcl2.k25,"d",c49.growth$cat.output.cl2)
c49.growth$cat.CNCI.cl2=ifelse(c49.growth$ch.meanCNCI.cl2>M49.CNCIcl2.n25,"i","n")
c49.growth$cat.CNCI.cl2=ifelse(c49.growth$ch.meanCNCI.cl2<M49.CNCIcl2.k25,"d",c49.growth$cat.CNCI.cl2)

c49.growth$cat.cl2=paste(c49.growth$cat.output.cl2,c49.growth$cat.CNCI.cl2,sep=";")
c49.growth$cat.cl2=ifelse(c49.growth$cat.cl2=="n;n","n","")

c49.growth$cat2.output.cl2=ifelse(c49.growth$ch.output.cl2>M49.outputcl2,"i","d")
c49.growth$cat2.CNCI.cl2=ifelse(c49.growth$ch.meanCNCI.cl2>M49.CNCIcl2,"i","d")
c49.growth$cat2.cl2=paste(c49.growth$cat2.output.cl2,c49.growth$cat2.CNCI.cl2,sep=";")
c49.growth$cat3.cl2=ifelse(c49.growth$cat.cl2=="n","n",c49.growth$cat2.cl2)
c49.growth$cat.cl2=c49.growth$cat3.cl2
c49.growth=c49.growth[,-c(10,11,13:16)]

#city name add
cities=c1.cityid.mregio[,-2]
cities=unique(cities) #3234
c=unique(c1.cityid.mregio$city.id)
rm(c)
c49.growth=merge(c49.growth,cities,by.x="city.id",by.y="city.id",all.x=T,all.y=F,sort=F)
View(c49.growth)

#Figure02; Figure04; Figure05; Figure06; based on Mathiensen et al., 2010, Figure6
#point graf; 
#https://r-charts.com/color-palette-generator/
c49.growth=arrange(c49.growth,(cat))
unique(c49.growth$cat)
# x axe: c49.growth$ch.meanCNCI  
# y axe: c49.growth$ch.output
#axes cross point M49.CNCI, M.49.output
# colors by c49.growth$cat (i;i  )
#"n"   c(col "grey")      "#C9C7C3"  
#"i;i" c(col "red")       "#EE2C2C"
#"i;d" c(col "darkgreen") "#43CD80"
#"d;i" c(col "gold")      "#FFD700"    "#FFB90F"
#"d;d" c(col "darkblue")  "#27408B"         "#4682B4" "#27408B"
 

 cols <- c("#27408B","#FFB90F","#43CD80","#EE2C2C","#C9C7C3")
 # labels of point c49.growth$city_ascii2
 library(ggplot2)
 library(ggrepel)
 # Define Colours
 cols <- c(
   "n"   = "#A9A9A9",
   "i;i" = "#EE2C2C",
   "i;d" = "#43CD80",
   "d;i" = "#FFD700",
   "d;d" = "#27408B"
 )
 
 #point graf SC
 
 g <- ggplot(c49.growth, aes(x = ch.output, y = ch.meanCNCI, color = cat, label = city_ascii2)) +  
   geom_point(size = 3, alpha = 0.8) +
   
   # Feliratok GGREPEL segítségével
   geom_text_repel(
     size = 3,
     force = 1,            # Mennyire tolja el a feliratokat egymástól
     max.overlaps = 30,    # Növelheted, ha sok feliratot akar mindenképp kitenni
     box.padding = 0.2,    # Távolság a pont és a szöveg között
     point.padding = 0.1,  # Üres hely a pont körül
     segment.color = "grey70", # Kis vonal, ami összeköti a szöveget a ponttal
     segment.linewidth = 0.2,
     show.legend = FALSE
   ) +
   
   scale_color_manual(
     values = cols,
     labels = c(
       "d;d" = "Black holes",
       "i;i" = "Hot spots",
       "i;d" = "Volume growth but loss in reputation",
       "d;i" = "Focus on succes",
       "n"   = "Neutral"
     )
   ) +
   
   # Átlagoknál a tengelyek metszéspontja
   geom_vline(xintercept = M49.output, linetype = "solid", color = "black", linewidth = 0.6) +
   geom_hline(yintercept = M49.CNCI, linetype = "solid", color = "black", linewidth = 0.6) +
   
   scale_x_continuous(breaks = seq(-10000, 10000, 100)) +
   scale_y_continuous(breaks = seq(-10000, 10000, 100)) +
   
   coord_fixed(ratio = 1) +
   
   theme_minimal() +
   labs(
     x = "Scientific Output Change (%), 2010-2017 to 2018-2022",
     y = "Scientific Impact Change (%), 2010-2017 to 2018-2022"
   ) +
   theme(
     panel.grid.major = element_line(color = "grey92", linewidth = 0.2), 
     panel.grid.minor = element_line(color = "grey95", linewidth = 0.1),
     
     legend.position = "bottom",
     legend.box = "horizontal",
     legend.title = element_blank(),
     axis.line = element_line(color = "black", linewidth = 0.5),
     axis.ticks = element_line(color = "black"),
     axis.ticks.length = unit(0.2, "cm"),
     plot.title = element_text(hjust = 0.5)
   )
 
 # Megjelenítés és mentés
 print(g)
 ggsave("gv2.pdf", plot = g, width = 25, height = 25, units = "cm")
 
 #############
 #SC cluster 1: output, impact change point graf 
 
 
 g.cl1 <- ggplot(c49.growth, aes(x = ch.output.cl1, y = ch.meanCNCI.cl1, color = cat.cl1, label = city_ascii2)) +  
   geom_point(size = 3, alpha = 0.8) +
   
   # Feliratok GGREPEL segítségével
   geom_text_repel(
     size = 3,
     force = 1,            # Mennyire tolja el a feliratokat egymástól
     max.overlaps = 30,    # Növelheted, ha sok feliratot akar mindenképp kitenni
     box.padding = 0.2,    # Távolság a pont és a szöveg között
     point.padding = 0.1,  # Üres hely a pont körül
     segment.color = "grey70", # Kis vonal, ami összeköti a szöveget a ponttal
     segment.linewidth = 0.2,
     show.legend = FALSE
   ) +
   
   scale_color_manual(
     values = cols,
     labels = c(
       "d;d" = "Black holes",
       "i;i" = "Hot spots",
       "i;d" = "Volume growth but loss in reputation",
       "d;i" = "Focus on succes",
       "n"   = "Neutral"
     )
   ) +
   
   # Átlagoknál a tengelyek metszéspontja
   geom_vline(xintercept = M49.outputcl1, linetype = "solid", color = "black", linewidth = 0.6) +
   geom_hline(yintercept = M49.CNCIcl1, linetype = "solid", color = "black", linewidth = 0.6) +
   
   scale_x_continuous(breaks = seq(-10000, 10000, 100)) +
   scale_y_continuous(breaks = seq(-10000, 10000, 100)) +
   
   coord_fixed(ratio = 1) +
   
   theme_minimal() +
   labs(
     x = "Scientific Output Change (%), 2010-2017 to 2018-2022",
     y = "Scientific Impact Change (%), 2010-2017 to 2018-2022"
   ) +
   theme(
     panel.grid.major = element_line(color = "grey92", linewidth = 0.2), 
     panel.grid.minor = element_line(color = "grey95", linewidth = 0.1),
     
     legend.position = "bottom",
     legend.box = "horizontal",
     legend.title = element_blank(),
     axis.line = element_line(color = "black", linewidth = 0.5),
     axis.ticks = element_line(color = "black"),
     axis.ticks.length = unit(0.2, "cm"),
     plot.title = element_text(hjust = 0.5)
   )
 
 # Megjelenítés és mentés
 print(g.cl1)
 ggsave("g.cl1.pdf", plot = g.cl1, width = 25, height = 25, units = "cm")
 
 #########-------------
 #SC cluster 2: output, impact change point graf
 
 g.cl2 <- ggplot(c49.growth, aes(x = ch.output.cl2, y = ch.meanCNCI.cl2, color = cat.cl2, label = city_ascii2)) +  
   geom_point(size = 3, alpha = 0.8) +
   
   # Feliratok GGREPEL segítségével
   geom_text_repel(
     size = 3,
     force = 1,            # Mennyire tolja el a feliratokat egymástól
     max.overlaps = 30,    # Növelheted, ha sok feliratot akar mindenképp kitenni
     box.padding = 0.2,    # Távolság a pont és a szöveg között
     point.padding = 0.1,  # Üres hely a pont körül
     segment.color = "grey70", # Kis vonal, ami összeköti a szöveget a ponttal
     segment.linewidth = 0.2,
     show.legend = FALSE
   ) +
   
   scale_color_manual(
     values = cols,
     labels = c(
       "d;d" = "Black holes",
       "i;i" = "Hot spots",
       "i;d" = "Volume growth but loss in reputation",
       "d;i" = "Focus on succes",
       "n"   = "Neutral"
     )
   ) +
   
   # Átlagoknál a tengelyek metszéspontja
   geom_vline(xintercept = M49.outputcl2, linetype = "solid", color = "black", linewidth = 0.6) +
   geom_hline(yintercept = M49.CNCIcl2, linetype = "solid", color = "black", linewidth = 0.6) +
   
   scale_x_continuous(breaks = seq(-10000, 10000, 100)) +
   scale_y_continuous(breaks = seq(-10000, 10000, 100)) +
   
   coord_fixed(ratio = 1) +
   
   theme_minimal() +
   labs(
     x = "Scientific Output Change (%), 2010-2017 to 2018-2022",
     y = "Scientific Impact Change (%), 2010-2017 to 2018-2022"
   ) +
   theme(
     panel.grid.major = element_line(color = "grey92", linewidth = 0.2), 
     panel.grid.minor = element_line(color = "grey95", linewidth = 0.1),
     
     legend.position = "bottom",
     legend.box = "horizontal",
     legend.title = element_blank(),
     axis.line = element_line(color = "black", linewidth = 0.5),
     axis.ticks = element_line(color = "black"),
     axis.ticks.length = unit(0.2, "cm"),
     plot.title = element_text(hjust = 0.5)
   )
 
 # Megjelenítés és mentés
 print(g.cl2)
 ggsave("g.cl2.pdf", plot = g.cl2, width = 25, height = 25, units = "cm")
 
 ##########cl1
 #Prague (Focus on succes category), Sydney (Hot spot category) and Guanzhou (Volume growth but loss in reputation category) are extreme change
 #Without them the same figure
 ol=c("Prague","Sydney","Guangzhou")
 c49.growth.woutliers=subset(c49.growth, !(city_ascii2 %in% ol))
 ############-
 #g.cl1 the same just 46 cities without outliers
 
 g.cl1.v2 <- ggplot(c49.growth.woutliers, aes(x = ch.output.cl1, y = ch.meanCNCI.cl1, color = cat.cl1, label = city_ascii2)) +  
   geom_point(size = 3, alpha = 0.8) +
   
   # Feliratok GGREPEL segítségével
   geom_text_repel(
     size = 3,
     force = 1,            # Mennyire tolja el a feliratokat egymástól
     max.overlaps = 30,    # Növelheted, ha sok feliratot akar mindenképp kitenni
     box.padding = 0.2,    # Távolság a pont és a szöveg között
     point.padding = 0.1,  # Üres hely a pont körül
     segment.color = "grey70", # Kis vonal, ami összeköti a szöveget a ponttal
     segment.linewidth = 0.2,
     show.legend = FALSE
   ) +
   
   scale_color_manual(
     values = cols,
     labels = c(
       "d;d" = "Black holes",
       "i;i" = "Hot spots",
       "i;d" = "Volume growth but loss in reputation",
       "d;i" = "Focus on succes",
       "n"   = "Neutral"
     )
   ) +
   
   # Átlagoknál a tengelyek metszéspontja
   geom_vline(xintercept = M49.outputcl1, linetype = "solid", color = "black", linewidth = 0.6) +
   geom_hline(yintercept = M49.CNCIcl1, linetype = "solid", color = "black", linewidth = 0.6) +
   
   scale_x_continuous(breaks = seq(-10000, 10000, 100)) +
   scale_y_continuous(breaks = seq(-10000, 10000, 100)) +
   
   coord_fixed(ratio = 1) +
   
   theme_minimal() +
   labs(
     x = "Scientific Output Change (%), 2010-2017 to 2018-2022",
     y = "Scientific Impact Change (%), 2010-2017 to 2018-2022"
   ) +
   theme(
     panel.grid.major = element_line(color = "grey92", linewidth = 0.2), 
     panel.grid.minor = element_line(color = "grey95", linewidth = 0.1),
     
     legend.position = "bottom",
     legend.box = "horizontal",
     legend.title = element_blank(),
     axis.line = element_line(color = "black", linewidth = 0.5),
     axis.ticks = element_line(color = "black"),
     axis.ticks.length = unit(0.2, "cm"),
     plot.title = element_text(hjust = 0.5)
   )
 
 # Megjelenítés és mentés
 print(g.cl1.v2)
 ggsave("g.cl1.v2.pdf", plot = g.cl1.v2, width = 25, height = 25, units = "cm")
########################
 
 #download WOS, export, import Incites (Analyse/RA, all years Schemes/Gipp overall; top5%; top10%)
 ######TABLE1, TABLE2, TABLE3; incite: https://access.clarivate.com/login?app=incites
 #p1.imp.all, p2.imp.all, instead of UTs with randomized ids of the publication and std CNCI 
 #------------------------------------------------------------------------------------
 
 #cluster and reg.collab.categories added impact values
 rm(p1.cl2,p2.cl2)
 p1.cl2=merge(kw.p.cl2.p1,p1.imp.all,by.x="UT",by.y="UT",all.x=F,all.y=T)
 p2.cl2=merge(kw.p.cl2.p2,p2.imp.all,by.x="UT",by.y="UT",all.x=F,all.y=T) #in Incites 4 pubs not fiend with CNCI, these are erised
 
 #AFFIL geography info added
 
 p1.cl.aff=merge(p1.cl2,c1.cityid.mregio,by.x="UT",by.y="UT",all.x=T,all.y=F) #most csak a intercity collab. 
 p2.cl.aff=merge(p2.cl2,c1.cityid.mregio,by.x="UT",by.y="UT",all.x=T,all.y=F) #most csak a intercity collab. 
 #output for all SC pubs cities
 
 p1.output.all=plyr::count(p1.cl.aff, vars=c("city.id"))
 colnames(p1.output.all)=c("city.id","output")
 
 #top5, top10 pubs number
 top5=aggregate(p1.cl.aff$top5, by=list(p1.cl.aff$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p1.output.all=merge(p1.output.all,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p1.cl.aff$top10, by=list(p1.cl.aff$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p1.output.all=merge(p1.output.all,top10,by.x="city.id",by.y="city.id",all=T)
 
 
 cnci=aggregate(p1.cl.aff$CNCI, by=list(p1.cl.aff$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI")
 p1.output.all=merge(p1.output.all,cnci,by.x="city.id",by.y="city.id",all=T)
 p1.output.all$mean.CNCI=p1.output.all$szum.CNCI/p1.output.all$output
 colnames(p1.output.all)=c("city.id",   "output.p1",    "output.top5.p1","output.top10.p1","szum.CNCI.p1", "mean.CNCI.p1")
 
 p2.output.all=plyr::count(p2.cl.aff, vars=c("city.id"))
 colnames(p2.output.all)=c("city.id","output")
 #top5, top10 pubs number
 top5=aggregate(p2.cl.aff$top5, by=list(p2.cl.aff$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p2.output.all=merge(p2.output.all,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p2.cl.aff$top10, by=list(p2.cl.aff$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p2.output.all=merge(p2.output.all,top10,by.x="city.id",by.y="city.id",all=T)
 
 cnci=aggregate(p2.cl.aff$CNCI, by=list(p2.cl.aff$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI")
 p2.output.all=merge(p2.output.all,cnci,by.x="city.id",by.y="city.id",all=T)
 p2.output.all$mean.CNCI=p2.output.all$szum.CNCI/p2.output.all$output
 colnames(p2.output.all)=c("city.id","output.p2", "output.top5.p2","output.top10.p2","szum.CNCI.p2", "mean.CNCI.p2")
 
 
 imp.output.all=merge(p1.output.all,p2.output.all,by.x="city.id",by.y="city.id",all=T)

 #2010-2022 together
 imp.output.all$output.p1=ifelse(is.na(imp.output.all$output.p1)==T,0,imp.output.all$output.p1)
 imp.output.all$output.p2=ifelse(is.na(imp.output.all$output.p2)==T,0,imp.output.all$output.p2)
 imp.output.all$output.top5.p1=ifelse(is.na(imp.output.all$output.top5.p1)==T,0,imp.output.all$output.top5.p1)
 imp.output.all$output.top5.p2=ifelse(is.na(imp.output.all$output.top5.p2)==T,0,imp.output.all$output.top5.p2)
 imp.output.all$output.top10.p1=ifelse(is.na(imp.output.all$output.top10.p1)==T,0,imp.output.all$output.top10.p1)
 imp.output.all$output.top10.p2=ifelse(is.na(imp.output.all$output.top10.p2)==T,0,imp.output.all$output.top10.p2)
 imp.output.all$szum.CNCI.p1=ifelse(is.na(imp.output.all$szum.CNCI.p1)==T,0,imp.output.all$szum.CNCI.p1)
 imp.output.all$szum.CNCI.p2=ifelse(is.na(imp.output.all$szum.CNCI.p2)==T,0,imp.output.all$szum.CNCI.p2)
 imp.output.all$mean.CNCI.p1=ifelse(is.na(imp.output.all$mean.CNCI.p1)==T,0,imp.output.all$mean.CNCI.p1)
 imp.output.all$mean.CNCI.p2=ifelse(is.na(imp.output.all$mean.CNCI.p2)==T,0,imp.output.all$mean.CNCI.p2)
 
 
 imp.output.all$output=imp.output.all$output.p1+imp.output.all$output.p2
 imp.output.all$top5=imp.output.all$output.top5.p1+imp.output.all$output.top5.p2
 imp.output.all$top10=imp.output.all$output.top10.p1+imp.output.all$output.top10.p2
 imp.output.all$szumCNCI=imp.output.all$szum.CNCI.p1+imp.output.all$szum.CNCI.p2
 imp.output.all$meanCNCI=imp.output.all$szumCNCI/imp.output.all$output
 View(imp.output.all)
 
 #cl1 and cl2 calculate output, szum.CNCI,mean.CNCI by two periods 2010-2017, 2018-2022
 
 
 #output for CL1 (Sustainability cluster) 
 p1.city.cl1=subset(p1.cl.aff,clus==1)
 p1.city.cl2=subset(p1.cl.aff,clus==2)
 p2.city.cl1=subset(p2.cl.aff,clus==1)
 p2.city.cl2=subset(p2.cl.aff,clus==2)
 
 #period1, period2 cl1,cl2
 #p1 cl1
 p1.output=plyr::count(p1.city.cl1, vars=c("city.id"))
 colnames(p1.output)=c("city.id","output.cl1")
 
 #top5, top10 pubs number
 top5=aggregate(p1.city.cl1$top5, by=list(p1.city.cl1$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p1.output=merge(p1.output,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p1.city.cl1$top10, by=list(p1.city.cl1$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p1.output=merge(p1.output,top10,by.x="city.id",by.y="city.id",all=T)
 
 cnci=aggregate(p1.city.cl1$CNCI, by=list(p1.city.cl1$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI.cl1")
 p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
 p1.output$mean.CNCI.cl1=p1.output$szum.CNCI.cl1/p1.output$output.cl1
 colnames(p1.output)=c("city.id",   "output.p1.cl1", "output.top5.p1.cl1", "output.top10.p1.cl1",   "szum.CNCI.p1.cl1", "mean.CNCI.p1.cl1")
 p1.output.cl1=p1.output
 #p1 cl2
 p1.output=plyr::count(p1.city.cl2, vars=c("city.id"))
 colnames(p1.output)=c("city.id","output.cl2")
 
 #top5, top10 pubs number
 top5=aggregate(p1.city.cl2$top5, by=list(p1.city.cl2$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p1.output=merge(p1.output,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p1.city.cl2$top10, by=list(p1.city.cl2$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p1.output=merge(p1.output,top10,by.x="city.id",by.y="city.id",all=T)
 
 cnci=aggregate(p1.city.cl2$CNCI, by=list(p1.city.cl2$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI.cl2")
 p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
 p1.output$mean.CNCI.cl2=p1.output$szum.CNCI.cl2/p1.output$output.cl2
 colnames(p1.output)=c("city.id",   "output.p1.cl2", "output.top5.p1.cl2", "output.top10.p1.cl2",   "szum.CNCI.p1.cl2", "mean.CNCI.p1.cl2")
 p1.output.cl2=p1.output
 
 #p2 cl1
 p1.output=plyr::count(p2.city.cl1, vars=c("city.id"))
 colnames(p1.output)=c("city.id","output.cl1")
 
 #top5, top10 pubs number
 top5=aggregate(p2.city.cl1$top5, by=list(p2.city.cl1$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p1.output=merge(p1.output,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p2.city.cl1$top10, by=list(p2.city.cl1$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p1.output=merge(p1.output,top10,by.x="city.id",by.y="city.id",all=T)
 
 cnci=aggregate(p2.city.cl1$CNCI, by=list(p2.city.cl1$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI.cl1")
 p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
 p1.output$mean.CNCI.cl1=p1.output$szum.CNCI.cl1/p1.output$output.cl1
 colnames(p1.output)=c("city.id",   "output.p2.cl1", "output.top5.p2.cl1", "output.top10.p2.cl1",   "szum.CNCI.p2.cl1", "mean.CNCI.p2.cl1")
 p2.output.cl1=p1.output
 #p2 cl2
 p1.output=plyr::count(p2.city.cl2, vars=c("city.id"))
 colnames(p1.output)=c("city.id","output.cl2")
 
 #top5, top10 pubs number
 top5=aggregate(p2.city.cl2$top5, by=list(p2.city.cl2$city.id), FUN=sum)
 colnames(top5)=c("city.id","output.top5")
 p1.output=merge(p1.output,top5,by.x="city.id",by.y="city.id",all=T)
 
 top10=aggregate(p2.city.cl2$top10, by=list(p2.city.cl2$city.id), FUN=sum)
 colnames(top10)=c("city.id","output.top10")
 p1.output=merge(p1.output,top10,by.x="city.id",by.y="city.id",all=T)
 
 cnci=aggregate(p2.city.cl2$CNCI, by=list(p2.city.cl2$city.id), FUN=sum)
 colnames(cnci)=c("city.id","szum.CNCI.cl2")
 p1.output=merge(p1.output,cnci,by.x="city.id",by.y="city.id",all=T)
 p1.output$mean.CNCI.cl2=p1.output$szum.CNCI.cl2/p1.output$output.cl2
 colnames(p1.output)=c("city.id",   "output.p2.cl2", "output.top5.p2.cl2", "output.top10.p2.cl2",   "szum.CNCI.p2.cl2", "mean.CNCI.p2.cl2")
 p2.output.cl2=p1.output
 
 imp.output.allcl=merge(imp.output.all,p1.output.cl1,by.x="city.id",by.y="city.id",all=T)
 imp.output.allcl=merge(imp.output.allcl,p2.output.cl1,by.x="city.id",by.y="city.id",all=T)
 imp.output.allcl=merge(imp.output.allcl,p1.output.cl2,by.x="city.id",by.y="city.id",all=T)
 imp.output.allcl=merge(imp.output.allcl,p2.output.cl2,by.x="city.id",by.y="city.id",all=T)
 
 #
 #-------------------------------------------
 #Figure01
 
 #no min output:
 map2=merge(imp.output.all,cities,by.x="city.id",by.y="city.id",all.x=T,all.y=F)
 
 #create map bars output 2010-2022
 map2.s=map2[,c("city.id","city_ascii2","lat","lng","output","szumCNCI","meanCNCI")]
 
 map2.s$lat=as.numeric(map2.s$lat)
 map2.s$lng=as.numeric(map2.s$lng)
 map2$lat=as.numeric(map2$lat)
 map2$lng=as.numeric(map2$lng)
 
 #the same period1
 #period 1 ahol output0 az ne legyen
 map1.p1=subset(map2,output.p1>0)
 #FIGURE1-------------------------------------------------------
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.12         # Magassági szorzó rate to 2010-2022 max value *3
 bar_fill <- "blue4"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p1, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (output.p1 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "output2010-2017.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 
 #mean CNCI
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.6         # Magassági szorzó rate to 2010-2022 max value /3
 bar_fill <- "#d63031"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p1, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (mean.CNCI.p1 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "meanCNCI2010-2017.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #max value have approx. the same high
 #period 2
 
 #period 2 ahol output0 az ne legyen
 map1.p2=subset(map2,output.p2>0)
 
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.056         # Magassági szorzó rate to 2010-2022 max value *1.4
 bar_fill <- "blue4"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p2, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (output.p2 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "output2018-2022.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #mean CNCI
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 2.8         # Magassági szorzó rate to 2010-2022 max value *1.4
 bar_fill <- "#d63031"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p2, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (mean.CNCI.p2 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "meanCNCI2018-2022.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 
 #----------------------------------------------------------------------------
 #Cluster1 Socio-oriented SC cluster period1 and period2
 
 #the same period1
 #period 1 ahol output0 az ne legyen
 map1.p1=subset(map2,output.p1>0)
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.12         # Magassági szorzó rate to 2010-2022 max value *3
 bar_fill <- "blue4"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p1, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (output.p1 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "output2010-2017.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #mean CNCI
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.6         # Magassági szorzó rate to 2010-2022 max value /3
 bar_fill <- "#d63031"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p1, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (mean.CNCI.p1 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "meanCNCI2010-2017.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #max value have approx. the same high
 #period 2
 
 #period 2 ahol output0 az ne legyen
 map1.p2=subset(map2,output.p2>0)
 
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 0.056         # Magassági szorzó rate to 2010-2022 max value *1.4
 bar_fill <- "blue4"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p2, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (output.p2 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "output2018-2022.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #mean CNCI
 # 1. Világtérkép lekérése
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Beállítások
 bar_w <- 0.35          # Vékony oszlopok
 h_scale <- 2.8         # Magassági szorzó rate to 2010-2022 max value *1.4
 bar_fill <- "#d63031"  # Oszlop színe
 bar_border <- "black"  # Oszlop kerete
 country_border <- "gray40" # Sötétebb, kontrasztos országhatárok
 
 # 3. Grafika elkészítése
 p <- ggplot() +
   # Kontinensek: Világos kitöltés, sötétebb és határozottabb országhatárok
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + # Vastagabb és sötétebb vonal az országhatároknak
   
   # Oszlopok vékony kerettel
   geom_rect(data = map1.p2, aes(
     xmin = lng - bar_w / 2,
     xmax = lng + bar_w / 2,
     ymin = lat,
     ymax = lat + (mean.CNCI.p2 * h_scale)
   ), 
   fill = bar_fill, 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # Téma beállításai (Fehér óceán)
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA)
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Exportálás 300dpi TIFF-be
 ggsave(
   filename = "meanCNCI2018-2022.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #-------------------------
 #FIGURE03
 #maps for RQ2 output, impact (mean CNCI) for all cities, cl1, cl2, P1, P2 FIGURE3
 
 #no min output:
 #MAP imp.output.all at least 10 pubs in a city during 2010-2022
 map3=merge(imp.output.allcl,cities,by.x="city.id",by.y="city.id",all.x=T,all.y=F)
 
 
 #create map bars output Period1 cl1 and cl2 2010-2017
 map3.s=map3[,c("city.id","city_ascii2","lat","lng","output.p1.cl1","output.p1.cl2","output.p2.cl1","output.p2.cl2","mean.CNCI.p1.cl1","mean.CNCI.p1.cl2","mean.CNCI.p2.cl1","mean.CNCI.p2.cl2")]
 #NA's csere 0-ra
 map3.s$output.p1.cl1=ifelse(is.na(map3.s$output.p1.cl1)==T,0,map3.s$output.p1.cl1)
 map3.s$output.p1.cl2=ifelse(is.na(map3.s$output.p1.cl2)==T,0,map3.s$output.p1.cl2)
 map3.s$output.p2.cl1=ifelse(is.na(map3.s$output.p2.cl1)==T,0,map3.s$output.p2.cl1)
 map3.s$output.p2.cl2=ifelse(is.na(map3.s$output.p2.cl2)==T,0,map3.s$output.p2.cl2)
 map3.s$mean.CNCI.p1.cl1=ifelse(is.na(map3.s$mean.CNCI.p1.cl1)==T,0,map3.s$mean.CNCI.p1.cl1)
 map3.s$mean.CNCI.p1.cl2=ifelse(is.na(map3.s$mean.CNCI.p1.cl2)==T,0,map3.s$mean.CNCI.p1.cl2)
 map3.s$mean.CNCI.p2.cl1=ifelse(is.na(map3.s$mean.CNCI.p2.cl1)==T,0,map3.s$mean.CNCI.p2.cl1)
 map3.s$mean.CNCI.p2.cl2=ifelse(is.na(map3.s$mean.CNCI.p2.cl2)==T,0,map3.s$mean.CNCI.p2.cl2)
 
 map3.s$lat=as.numeric(map3.s$lat)
 map3.s$lng=as.numeric(map3.s$lng)
 summary(map3.s)
 
 #period1 cites
 map3.p1=subset(map3.s, city.id %in% nc1) #1566 cities
 map3.p2=subset(map3.s, city.id %in% nc2) #2836 cities
 
 library(ggplot2)
 library(sf)
 library(rnaturalearth)
 library(rnaturalearthdata)
 library(ggrepel) 
 
 #OUTPUT PERIOD 1 CL1 AND CL2
 # 1. Load world map
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Settings
 bar_w <- 0.30          
 h_scale <- 0.24         
 bar_border <- "black"
 country_border <- "gray40"
 
 # 3. Create Graphic
 p <- ggplot() +
   # Map Base
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + 
   
   # --- BAR 1 (Social-oriented) ---
   geom_rect(data = map3.p1, aes(
     xmin = lng - bar_w,       
     xmax = lng,               
     ymin = lat,
     ymax = lat + (output.p1.cl1 * h_scale),
     fill = "Social-oriented SC sub-field"  # FIXED: Moved inside aes()
   ), 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # --- BAR 2 (Technology-oriented) ---
   geom_rect(data = map3.p1, aes(
     xmin = lng,               
     xmax = lng + bar_w,       
     ymin = lat,
     ymax = lat + (output.p1.cl2 * h_scale),
     fill = "Technology-oriented SC sub-field" # FIXED: Moved inside aes()
   ), 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # --- Manual Color Scale ---
   scale_fill_manual(values = c(
     "Social-oriented SC sub-field" = "green",
     "Technology-oriented SC sub-field" = "royalblue3"
   )) +
   
   # Theme Settings
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA),
     legend.position = "bottom",          
     legend.title = element_blank(),       
     legend.text = element_text(size = 10),
     legend.key = element_blank() # Recommended: removes grey boxes behind legend icons
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Export
 ggsave(
   filename = "output.p1.cl1cl2.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 #-------------------------------------------------------
 #OUTPUT PERIOD2 OUTPUT cl1 and cl2
 
 
 
 # 1. Load world map
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Settings
 bar_w <- 0.30          
 h_scale <- 0.08         
 bar_border <- "black"
 country_border <- "gray40"
 
 # 3. Create Graphic
 p <- ggplot() +
   # Map Base
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + 
   
   # --- BAR 1 (Social-oriented) ---
   geom_rect(data = map3.p2, aes(
     xmin = lng - bar_w,       
     xmax = lng,               
     ymin = lat,
     ymax = lat + (output.p2.cl1 * h_scale),
     fill = "Social-oriented SC sub-field"  # FIXED: Moved inside aes()
   ), 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # --- BAR 2 (Technology-oriented) ---
   geom_rect(data = map3.p2, aes(
     xmin = lng,               
     xmax = lng + bar_w,       
     ymin = lat,
     ymax = lat + (output.p2.cl2 * h_scale),
     fill = "Technology-oriented SC sub-field" # FIXED: Moved inside aes()
   ), 
   color = bar_border, 
   linewidth = 0.05) + 
   
   # --- Manual Color Scale ---
   scale_fill_manual(values = c(
     "Social-oriented SC sub-field" = "green",
     "Technology-oriented SC sub-field" = "royalblue3"
   )) +
   
   # Theme Settings
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA),
     legend.position = "bottom",          
     legend.title = element_blank(),       
     legend.text = element_text(size = 10),
     legend.key = element_blank() # Recommended: removes grey boxes behind legend icons
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Export
 ggsave(
   filename = "output.p2.cl1cl2.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )###########-----------------------
 
 #MEAN CNCI PERIOD 1 CL1 AND CL2
 
 # 1. Load world map
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Settings
 bar_w <- 0.30          # Width of a single bar
 h_scale <- 0.3         # Height multiplier
 country_border <- "gray40"
 
 # 3. Create Graphic
 p <- ggplot() +
   # Map Base
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + 
   
   # --- BAR 1 (Social-oriented) ---
   geom_rect(data = map3.p1, aes(
     xmin = lng - bar_w,       
     xmax = lng,               
     ymin = lat,
     ymax = lat + (mean.CNCI.p1.cl1 * h_scale),
     fill = "Social-oriented SC sub-field" # Corrected: inside aes()
   ), 
   color = "black", 
   linewidth = 0.05) + 
   
   # --- BAR 2 (Technology-oriented) ---
   geom_rect(data = map3.p1, aes(
     xmin = lng,               
     xmax = lng + bar_w,       
     ymin = lat,
     ymax = lat + (mean.CNCI.p1.cl2 * h_scale),
     fill = "Technology-oriented SC sub-field" # Corrected: inside aes()
   ), 
   color = "black", 
   linewidth = 0.05) + 
   
   # --- Manual Color Scale & Legend Labels ---
   scale_fill_manual(values = c(
     "Social-oriented SC sub-field" = "forestgreen",
     "Technology-oriented SC sub-field" = "tomato1"
   )) +
   
   # Theme Settings
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA),
     legend.position = "bottom",          # Legend at bottom
     legend.title = element_blank(),       # No title
     legend.text = element_text(size = 10),
     legend.key = element_blank()          # Removes gray background from legend icons
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Export
 ggsave(
   filename = "meanCNCI.p1.cl1cl2.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )
 
 #-------------------------------------------------------
 #MEAN CNCI PERIOD2 OUTPUT cl1 and cl2
 
 # 1. Load world map
 world <- ne_countries(scale = "medium", returnclass = "sf")
 
 # 2. Settings
 bar_w <- 0.20          # Width of a single bar
 h_scale <- 1.2         # Height multiplier
 country_border <- "gray40"
 
 # 3. Create Graphic
 p <- ggplot() +
   # Map Base
   geom_sf(data = world, 
           fill = "gray96", 
           color = country_border, 
           linewidth = 0.2) + 
   
   # --- BAR 1 (Social-oriented) ---
   geom_rect(data = map3.p2, aes(
     xmin = lng - bar_w,       
     xmax = lng,               
     ymin = lat,
     ymax = lat + (mean.CNCI.p2.cl1 * h_scale),
     fill = "Social-oriented SC sub-field" # Corrected: inside aes()
   ), 
   color = "black", 
   linewidth = 0.05) + 
   
   # --- BAR 2 (Technology-oriented) ---
   geom_rect(data = map3.p2, aes(
     xmin = lng,               
     xmax = lng + bar_w,       
     ymin = lat,
     ymax = lat + (mean.CNCI.p2.cl2 * h_scale),
     fill = "Technology-oriented SC sub-field" # Corrected: inside aes()
   ), 
   color = "black", 
   linewidth = 0.05) + 
   
   # --- Manual Color Scale & Legend Labels ---
   scale_fill_manual(values = c(
     "Social-oriented SC sub-field" = "forestgreen",
     "Technology-oriented SC sub-field" = "tomato1"
   )) +
   
   # Theme Settings
   theme_void() +
   theme(
     panel.background = element_rect(fill = "white", color = NA),
     plot.background = element_rect(fill = "white", color = NA),
     legend.position = "bottom",          # Legend at bottom
     legend.title = element_blank(),       # No title
     legend.text = element_text(size = 10),
     legend.key = element_blank()          # Removes gray background from legend icons
   ) +
   coord_sf(expand = FALSE)
 
 # 4. Export
 ggsave(
   filename = "meanCNCI.p2.cl1cl2.all.tiff", 
   plot = p, 
   device = "tiff", 
   dpi = 300, 
   width = 14, 
   height = 7, 
   units = "in", 
   compression = "lzw"
 )