library(tidyr) library(dplyr) library(ggplot2) library(sf) library(skimr) library(ggrepel) library(stringr)

mys_adm0 <- st_read("mys_adm0.gpkg") mys_adm1 <- st_read("mys_adm1.gpkg") mys_adm2 <- st_read("mys_adm2.gpkg") googmob <- read.csv("googmob.csv")

mys_adm1long <- mys_adm1 %>% gather(month, light, dec19:jul21) head(mys_adm1long)

palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High", "Low-Low", "Low-High", "High-Low") #map ggplot() + geom_sf(data = mys_adm2, aes(fill = factor(sat_cluster)), lwd=0) + geom_sf(data = mys_adm1, aes(fill=NA, alpha=0), show.legend = FALSE, lwd=0.05) + labs(title = "Night-time lights, 2019") + scale_fill_manual(name ="Median LISAs clusters", values = palette, labels = labelss) + theme_set(theme_bw()) + theme(plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #Measured in radiance #save ggsave("LISAs SAT.png")

palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582') labelss <- c("Not Significant", "High-High", "Low-Low", "Low-High", "High-Low") #map ggplot() + geom_sf(data = mys_adm2 , aes(fill = factor(i_cluster)), lwd=0) + geom_sf(data = mys_adm1, aes(fill=NA, alpha=0), show.legend = FALSE, lwd=0.05) + labs(title = "Gross Regional Domestic Product (GRDP), 2019") + scale_fill_manual(name ="Median LISAs clusters", values = palette, labels = labelss) + theme_set(theme_bw()) + theme( plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) #Measured in constant 2015 price #save ggsave("LISAs GRDP.png", width = 6.67 , height = 6.67 )

#reshape mys_adm1long2 <- mys_adm1 %>% gather(month, sat, sat19:sat20) %>% gather(month, gdp, gdp19:gdp20) head(mys_adm1long2) #r2 lm19 <- lm(mys_adm1$gdp19 ~ mys_adm1$sat19) summary(lm19)

#r2 lm20 <- lm(mys_adm1$gdp20 ~ mys_adm1$sat20) summary(lm20)

#graph ggplot(data= mys_adm1long2, aes(x = log(sat), y = log(gdp), col = factor(month))) + geom_point() + geom_smooth(method=lm, se=FALSE) + geom_label(aes(label = state)) + scale_color_manual(name ="Year (R²)", values = c("#41ae76","#005824"), labels = c("2019 (0.61)", "2020 (0.61)")) + labs(title = "Proportion of variation in GRDP predictable by lights", x ="Night-time lights intensity, Radiance", y = "Gross Regional Domestic Product (GRDP), VND") + theme_set(theme_bw()) + theme_update(plot.title = element_text(hjust = 0.5)) ggsave("rel.png")

#map ggplot() + geom_sf(data = mys_adm2, aes(fill = localR2), alpha=0.5, lwd=0) + geom_sf(data = mys_adm1,alpha= 0, lwd=0.2) + labs(title = "Proportion of variation in GRDP predictable by lights") + theme_set(theme_bw()) + scale_fill_gradient(name = "R²", low = "white", high = "darkgreen") + theme(plot.title=element_text(hjust=0.5, lineheight=1.2), axis.title.x = element_blank(), axis.title.y = element_blank()) + geom_label_repel(data = mys_adm1, aes(x = COORD_X, y = COORD_Y, label = state), size=3, alpha =0.8) #save ggsave("grdp_r2.png")

# Compute Light percentage growth rate mys_adm0 <- mys_adm0 %>% mutate(feb20_gr = ((log(feb20/jan20)))*100) %>% mutate(mar20_gr = ((log(mar20 / feb20)))*100) %>% mutate(apr20_gr = ((log(apr20 / mar20)))*100) %>% mutate(may20_gr = ((log(may20 / apr20)))*100) %>% mutate(jun20_gr = ((log(jun20 / may20)))*100) %>% mutate(jul20_gr = ((log(jul20 / jun20)))*100) %>% mutate(aug20_gr = ((log(aug20 / jul20)))*100) %>% mutate(sep20_gr = ((log(sep20 / aug20)))*100) %>% mutate(oct20_gr = ((log(oct20 / sep20)))*100) %>% mutate(nov20_gr = ((log(nov20 / oct20)))*100) %>% mutate(dec20_gr = ((log(dec20 / nov20)))*100) %>% mutate(jan21_gr = ((log(jan21 / dec20)))*100) %>% mutate(feb21_gr = ((log(feb21 / jan21)))*100) %>% mutate(mar21_gr = ((log(mar21 / feb21)))*100) %>% mutate(apr21_gr = ((log(apr21 / mar21)))*100) %>% mutate(may21_gr = ((log(may21 / apr21)))*100) %>% mutate(jun21_gr = ((log(jun21 / may21)))*100) %>% mutate(jul21_gr = ((log(jul21 / jun21)))*100) mys_adm0l <- mys_adm0 %>% gather(month, light_gr, feb20_gr:jul21_gr) mys_adm0l$month <- strtrim(mys_adm0l$month,5) # Extract Mobility (Percentage change from baseline) googmob_adm0 <- googmob %>% filter(state=="Malaysia") # Compile GRDP (percentage growth rates) month <- c("feb20","mar20","apr20","may20","jun20","jul20","aug20","sep20","oct20","nov20", "dec20","jan21","feb21","mar21","apr21","may21","jun21","jul21") chg <- c(6,-5.6,-28.7,-19.7,-3,-2.8,-3.5,-1.8,-4.3,-4.0,-2.1,-3.6,-3.6,6.1,39.7,19.3,-3.6,-7.6) Malaysia <- c("Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia","Malaysia") monthly_gdp <- data.frame(Malaysia,chg,month) monthly_gdp <- reshape(monthly_gdp, idvar = "Malaysia", timevar = "month", direction = "wide") colnames(monthly_gdp) <- c("state","feb20","mar20","apr20","may20","jun20","jul20","aug20","sep20","oct20","nov20", "dec20","jan21","feb21","mar21","apr21","may21","jun21","jul21") monthly_gdp <- monthly_gdp %>% gather(month, g_gdp, feb20:jul21) # Plot ggplot() + geom_line(data = mys_adm0l, aes(month, light_gr, color = "2. Night-time lights", group=state)) + geom_line(data = monthly_gdp, aes(month, g_gdp, color = "1. GRDP", group=state)) + geom_line(data = googmob_adm0, aes(month, median, color = "3. Google mobility", group=state)) + scale_color_manual(name =" ", values = c("#41ae76","#ffc000","#649fff")) + theme_set(theme_bw()) + labs(title = "Comparing the dynamics of Gross Regional Domestic Product (GRDP), Night-time lights, and Google mobility in Peninsular Malaysia") + theme( plot.title=element_text(hjust=0.5, lineheight=1.2)) + scale_y_continuous(name ="Percentage change") + scale_x_discrete(name ="Time (month.year)", limits=c("feb20","mar20","apr20","may20","jun20","jul20","aug20","sep20","oct20","nov20","dec20","jan21","feb21","mar21","apr21","may21","jun21","jul21"), labels=c("feb20"="02.2020","mar20"="03.2020","apr20"="04.2020","may20"="05.2020","jun20"="06.2020","jul20"="07.2020","aug20"="08.2020","sep20"="09.2020","oct20"="10.2020","nov20"="11.2020","dec20"="12.2020","jan21"="01.2021","feb21"="02.2021","mar21"="03.2021","apr21"="04.2021","may21"="05.2021","jun21"="06.2021","jul21"="07.2021")) ggsave("national_lights_grdp_mobility.png",width = 15, height = 10)

# Light mys_adm1 <- mys_adm1 %>% mutate(feb20_gr = ((log(feb20/jan20)))*100) %>% mutate(mar20_gr = ((log(mar20 / feb20)))*100) %>% mutate(apr20_gr = ((log(apr20 / mar20)))*100) %>% mutate(may20_gr = ((log(may20 / apr20)))*100) %>% mutate(jun20_gr = ((log(jun20 / may20)))*100) %>% mutate(jul20_gr = ((log(jul20 / jun20)))*100) %>% mutate(aug20_gr = ((log(aug20 / jul20)))*100) %>% mutate(sep20_gr = ((log(sep20 / aug20)))*100) %>% mutate(oct20_gr = ((log(oct20 / sep20)))*100) %>% mutate(nov20_gr = ((log(nov20 / oct20)))*100) %>% mutate(dec20_gr = ((log(dec20 / nov20)))*100) %>% mutate(jan21_gr = ((log(jan21 / dec20)))*100) %>% mutate(feb21_gr = ((log(feb21 / jan21)))*100) %>% mutate(mar21_gr = ((log(mar21 / feb21)))*100) %>% mutate(apr21_gr = ((log(apr21 / mar21)))*100) %>% mutate(may21_gr = ((log(may21 / apr21)))*100) %>% mutate(jun21_gr = ((log(jun21 / may21)))*100) %>% mutate(jul21_gr = ((log(jul21 / jun21)))*100) mys_adm1l <- mys_adm1 %>% gather(month, light_gr, feb20_gr:jul21_gr) mys_adm1l$month <- strtrim(mys_adm1l$month,5) # Extract Mobility (Percentage change from baseline) googmob_adm1 <- googmob %>% filter(!state=="Malaysia") #Plot ggplot() + geom_line(data = mys_adm1l, aes(month, light_gr, color = "1. Night-time lights", group=state)) + geom_line(data = googmob_adm1, aes(month, median, color = "2. Google mobility", group=state)) + labs(title = "Comparing the dynamics of Night-time lights and Google mobility across states in Peninsular Malaysia") + theme( plot.title=element_text(hjust=0.5, lineheight=1.2)) + scale_color_manual(name =" ", values = c("#ffc000","#649fff")) + scale_y_continuous("Percentage change") + scale_x_discrete(name ="Time (month.year)", limits=c("feb20","mar20","apr20","may20","jun20","jul20","aug20","sep20","oct20","nov20","dec20","jan21","feb21","mar21","apr21","may21","jun21","jul21"), labels=c("feb20"="02.2020","mar20"="03.2020","apr20"="04.2020","may20"="05.2020","jun20"="06.2020","jul20"="07.2020","aug20"="08.2020","sep20"="09.2020","oct20"="10.2020","nov20"="11.2020","dec20"="12.2020","jan21"="01.2021","feb21"="02.2021","mar21"="03.2021","apr21"="04.2021","may21"="05.2021","jun21"="06.2021","jul21"="07.2021")) ggsave("state_lights_mobility.png",width = 15, height = 10)

#Plot col <- c('#41222A','#EC410B','#C14C32','#80003A','#B30019','#696B7E','#FFC000','#50B4D8','#506432','#CC0001','#445A67', '#010066') ggplot() + geom_line(data = mys_adm1l, aes(month, light_gr, color = state, group=state)) + labs(title = "Dynamics of night-time lights across states in Peninsular Malaysia") + theme( plot.title=element_text(hjust=0.5, lineheight=1.2)) + scale_color_manual(name =" ", values = col) + scale_y_continuous("Percentage change in night-time lights") + scale_x_discrete(name ="Time (month.year)", limits=c("feb20","mar20","apr20","may20","jun20","jul20","aug20","sep20","oct20","nov20","dec20","jan21","feb21","mar21","apr21","may21","jun21","jul21"), labels=c("feb20"="02.2020","mar20"="03.2020","apr20"="04.2020","may20"="05.2020","jun20"="06.2020","jul20"="07.2020","aug20"="08.2020","sep20"="09.2020","oct20"="10.2020","nov20"="11.2020","dec20"="12.2020","jan21"="01.2021","feb21"="02.2021","mar21"="03.2021","apr21"="04.2021","may21"="05.2021","jun21"="06.2021","jul21"="07.2021")) ggsave("state_lights.png",width = 15, height = 10)

ggplot(data=mys_adm1) + geom_sf(aes(fill = state), lwd=0) + scale_fill_manual( values = col) + theme_set(theme_bw()) + theme(legend.position = "none", panel.grid = element_blank(), axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.background = element_blank()) + geom_label_repel(aes(x = COORD_X, y = COORD_Y, label = state), size=3, alpha =0.8) ggsave("mys_map.png")