library(tidyverse) library(rgdal) library(rgeos) library(spgwr) library(RColorBrewer) library(ggrepel) library(estdaR)

#shapefile sh <- st_read("eap_subnational.gpkg") sh_c <- st_read("eap_country.gpkg") #skim dataset sh_df <- st_drop_geometry(sh) sh_df

lm20 <- lm(sh$MORAN_LAG ~ sh$MORAN_STD) summary(lm20)

#prep palette <- c('#C0C0C0', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant (49,712)", "High-High (3,810)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") #plot ggplot(sh, aes(x=sh$MORAN_STD, y=sh$MORAN_LAG)) + geom_point(aes(col=factor(LISA_CL))) + geom_smooth(formula=y ~ x, method="lm", se=FALSE) + geom_hline(yintercept=mean(sh$MORAN_LAG), lty=2) + geom_vline(xintercept=mean(sh$MORAN_STD), lty=2) + labs(title="Global Moran's I (0.845) Scatterplot", y="Average night-time lights intensity of neighbouring regions", x="Night-time lights intensity of location i", color="Clusters") + scale_color_manual(labels = labelss, values = palette) theme(plot.title=element_text(hjust=0.5, lineheight=1.2)) + theme_set(theme_bw()) #save ggsave("Global Moran's scatterplot.png")

#remove outliers sh_sub <- sh %>% filter(MORAN_STD < 1) %>% filter(MORAN_LAG < 1) #count number of outliers removed sh_sub_df <- st_drop_geometry(sh_sub) remaining <- sh_sub_df %>% filter (LISA_CL ==1) count(remaining)

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant (49,712)", "High-High (27)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") #plot ggplot(sh_sub, aes(x=sh_sub$MORAN_STD, y=sh_sub$MORAN_LAG)) + geom_point(aes(col=factor(LISA_CL))) + geom_smooth(formula=y ~ x, method="lm", se=FALSE) + geom_hline(yintercept=mean(sh$MORAN_LAG), lty=2) + geom_vline(xintercept=mean(sh$MORAN_STD), lty=2) + labs(title="(excld. 3,783 outliers)", y="Average night-time lights intensity of neighbouring regions", x="Night-time lights intensity of location i", color="Clusters") + scale_color_manual(labels = labelss, values = palette) theme(plot.title=element_text(hjust=0.5, lineheight=1.2)) + theme_set(theme_bw()) #save ggsave("Global Moran's scatterplot (excld outliers).png")

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High (3,810)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") ##not significant become white, borders become country border!! ##add numbers to legend #map ggplot() + geom_sf(data = sh, aes(fill = factor(LISA_CL)), lwd=0.2) + geom_sf(data = sh_c, fill = NA, aes(alpha=0), show.legend = FALSE) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(name ="Clusters", values = palette, labels = labelss) + #geom_text_repel(data = sh_sub, aes(x = X, y = Y, label = region), size=2, point.padding = NA) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save ggsave("LISAs.png") ##text title write ~ spatial distribution of night-time lights across east asia and pacific ##note: N=69,000

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High (3,810)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") #map ggplot() + geom_sf(data = sh, aes(fill = factor(LISA_CL)), lwd=0) + geom_sf(data = sh_c, fill = NA, aes(alpha=0), show.legend = FALSE) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(name ="Clusters", values = palette, labels = labelss) + #geom_text_repel(data = sh_sub, aes(x = X, y = Y, label = region), size=2, point.padding = NA) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save ggsave("LISAs.png") ##text title write ~ spatial distribution of night-time lights across east asia and pacific ##note: N=69,000

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High (3,810)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") ##not significant become white, borders become country border!! ##add numbers to legend #map ggplot() + geom_sf(data = sh, aes(fill = factor(LISA_CL)), color = NA) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(name ="Clusters", values = palette, labels = labelss) + #geom_text_repel(data = sh_sub, aes(x = X, y = Y, label = region), size=2, point.padding = NA) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) + geom_sf(data = sh_c, fill = NA, aes(alpha=0), show.legend = FALSE) #save ggsave("LISAs.png") ##text title write ~ spatial distribution of night-time lights across east asia and pacific ##note: N=69,000

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High (3,810)", "Low-Low (15,425)", "Low-High (33)", "High-Low (20)") ##not significant become white, borders become country border!! ##add numbers to legend #map ggplot() + geom_sf(data = sh, aes(fill = factor(LISA_CL)), color = NA) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(name ="Clusters", values = palette, labels = labelss) + #geom_text_repel(data = sh_sub, aes(x = X, y = Y, label = region), size=2, point.padding = NA) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save ggsave("LISAs.png") ##text title write ~ spatial distribution of night-time lights across east asia and pacific ##note: N=69,000

#prep !!!change to p-value specific palette <- c('#f2f2f2', '#b0dbf1', '#02198B') labelss <- c("Not Significant (49,7120)", "0.01 (12,028)", "0.001 (7,260)") #map ggplot() + geom_sf(data = sh, aes(fill = factor(P)), color = NA) + labs(title = "Local Moran's I") + scale_fill_manual(name ="P-value", values = palette, labels = labelss) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save ggsave("LISAs_p.png")

#prep palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High (7,450)", "Low-Low (8,255)", "Low-High (476)", "High-Low (19)") #map ggplot() + geom_sf(data = sh, aes(fill = factor(LISA_CL_M)), color = NA) + geom_sf(data = sh_c, fill = NA, aes(alpha=0), show.legend = FALSE, lwd=0.1) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(name ="Clusters", values = palette, labels = labelss) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save ggsave("Median LISAs.png")

#prep !!!change to p-value specific palette <- c('#f2f2f2', '#b0dbf1', '#02198B') labelss <- c("Not Significant (52,800)", "0.01 (8,570)", "0.001 (7,630)") #map ggplot() + geom_sf(data = sh, aes(fill = factor(P_M)), color = NA) + labs(title = "Local Moran's I (2018)") + scale_fill_manual(values = palette) + theme_set(theme_bw()) + theme(legend.title = element_blank(), plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #save #ggsave("Median LISAs p.png")

#compute spatial weights coord = data.frame(sh$X,sh$Y) knn <- knn2nb(knearneigh(as.matrix(coord), k = 6)) w <- nb2listw(knn, style = "W",zero.policy=TRUE) #standardise variables sat12r <- sh$Sat_2012 - mean(sh$Sat_2012) sat18r <- sh$Sat_2018 - mean(sh$Sat_2018) #compute directional LISAs t0 <- sat12r t1 <- sat18r Regime <- sh$Country dl <- d.LISA(t0, t1, w, Regime, k = 8, mean.rel = FALSE, nsim = 999, arrow = FALSE, only = NULL)

#plot a<-dl$lisamap b<-a$data b<-as.data.frame(b) col <- c('#e2e1a7','#e2e248','#77b99b','#8fc467','#8e8268','#dbeb83','#b6edb8','#deb057','#8c68d0','#59e695','#68ecdc','#dc4ad7','#e0704f','#cde4da','#db4982','#63ee57','#8d3de8','#8597d8','#d6c7dc','#dc78d1','#eabda3','#dea8e2','#75c0dc','#acec4d','#d3849c') ggplot(b, aes(x = x0, xend = x1, y = y0, yend = y1)) + geom_segment(aes(color = Regime), alpha=0.8, arrow = arrow(length = unit(0.1, "cm"))) + scale_color_manual(values = col) + labs(title = "Directional LISAs plot", color="Country", y="", x="") + theme_set(theme_bw())+ theme(plot.title=element_text(hjust=0.5, lineheight=1.2), axis.text.y = element_blank(), axis.text.x = element_blank()) ggsave("Directional LISAs plot.png")

a<- dl$rose b <- a$data angle <- b$angle angle =substr(angle,1,nchar(b$angle)-3) angle <-str_replace_all(angle, "-", "") angle <-as.numeric(angle) b$angle <- angle ticks.shdi <- c("57", "30", "12", "9", "73","42","9","12") ggplot(b, aes(x=angle, y=length, col=Regime)) + geom_col(width = 47) + coord_polar() + labs(title = "Rose Diagram", color="Country") + scale_color_manual(values = col) + theme_minimal() + theme( axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.y = element_blank(), axis.text.x = element_blank(), plot.subtitle = element_text(size = 12), plot.title = element_text(hjust = 0.5), legend.position = "none" )+ scale_x_discrete(breaks= ticks.shdi) ggsave("Directional LISAs rose.png")

a<- dl$p.value palette <- c('#f1eef6','#bdc9e1','#74a9cf','#0570b0') ggplot(a, aes(x=type, y=1, fill=p.value)) + geom_col(width = 1) + labs(title = " ", fill="P-value") + scale_fill_manual(values = palette ) + coord_polar() + theme_minimal() + theme( axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.y = element_blank(), axis.text.x = element_blank() ) ggsave("Directional LISAs p-value.png")

ggplot() + geom_sf(data = sh, aes(fill = Country), color = NA, show.legend = FALSE) + scale_fill_manual(values = col)