Load required libraries
library(tidyverse)
library(rgdal)
library(rgeos)
library(spgwr)
library(RColorBrewer)
library(ggrepel)
library(estdaR)
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
✔ ggplot2 3.3.3 ✔ purrr 0.3.4
✔ tibble 3.1.0 ✔ dplyr 1.0.4
✔ tidyr 1.1.2 ✔ stringr 1.4.0
✔ readr 1.4.0 ✔ forcats 0.5.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
Loading required package: sp
rgdal: version: 1.5-23, (SVN revision 1121)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
Path to GDAL shared files: /usr/share/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
Path to PROJ shared files: /root/.local/share/proj:/usr/share/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.4-5
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
rgeos version: 0.5-5, (SVN revision 640)
GEOS runtime version: 3.9.0-CAPI-1.16.2
Linking to sp version: 1.4-5
Polygon checking: TRUE
Loading required package: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
NOTE: This package does not constitute approval of GWR
as a method of spatial analysis; see example(gwr)
Loading required package: spdep
Loading required package: sf
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Attaching package: ‘estdaR’
The following objects are masked from ‘package:spdep’:
geary, moran
Load data
#shapefile
sh <- st_read("eap_subnational.gpkg")
sh_c <- st_read("eap_country.gpkg")
#skim dataset
sh_df <- st_drop_geometry(sh)
sh_df
Reading layer `satellite' from data source `/work/eap_subnational.gpkg' using driver `GPKG'
Simple feature collection with 69000 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 73.5577 ymin: -55.11694 xmax: 188.6023 ymax: 53.56086
Geodetic CRS: WGS 84
Reading layer `A' from data source `/work/eap_country.gpkg' using driver `GPKG'
Simple feature collection with 25 features and 67 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 73.5577 ymin: -55.11694 xmax: 188.6023 ymax: 53.56086
Geodetic CRS: WGS 84
Preview global Moran's I
lm20 <- lm(sh$MORAN_LAG ~ sh$MORAN_STD)
summary(lm20)
Plot global Moran's I scatterplot
note: values are standardised
#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")
Warning message:
“Use of `sh$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Warning message:
“Use of `sh$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Saving 6.67 x 6.67 in image
Warning message:
“Use of `sh$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Warning message:
“Use of `sh$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Plot global Moran's I scatterplot (excld. outliers)
#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")
Warning message:
“Use of `sh_sub$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh_sub$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Warning message:
“Use of `sh_sub$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh_sub$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Saving 6.67 x 6.67 in image
Warning message:
“Use of `sh_sub$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh_sub$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Warning message:
“Use of `sh_sub$MORAN_STD` is discouraged. Use `MORAN_STD` instead.”
Warning message:
“Use of `sh_sub$MORAN_LAG` is discouraged. Use `MORAN_LAG` instead.”
Plot local Moran's I map
#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
Saving 6.67 x 6.67 in image
Execution error
KernelInterrupted: Execution interrupted by the Jupyter kernel.
#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
Saving 6.67 x 6.67 in image
Execution error
KernelInterrupted: Execution interrupted by the Jupyter kernel.
#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
Saving 6.67 x 6.67 in image
P-value
#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")
Saving 6.67 x 6.67 in image
Execution error
KernelInterrupted: Execution interrupted by the Jupyter kernel.
Plot local Moran's I map (median)
#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")
Saving 6.67 x 6.67 in image
Execution error
KernelInterrupted: Execution interrupted by the Jupyter kernel.
P-value
#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 directional local Moran's I
#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)
Warning message in RColorBrewer::brewer.pal(n, pal):
“n too large, allowed maximum for palette Set1 is 9
Returning the palette you asked for with that many colors
”
Warning message in RColorBrewer::brewer.pal(n, pal):
“n too large, allowed maximum for palette Set1 is 9
Returning the palette you asked for with that many colors
”
Warning message:
“Removed 51269 rows containing missing values (geom_point).”
Plot directional local Moran's I (arrow)
#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")
Saving 6.67 x 6.67 in image
Plot directional local Moran's I (rose)
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")
Warning message:
“position_stack requires non-overlapping x intervals”
Saving 6.67 x 6.67 in image
Warning message:
“position_stack requires non-overlapping x intervals”
P-value
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")
Saving 6.67 x 6.67 in image
Plot East Asia and Pacific Map
ggplot() +
geom_sf(data = sh, aes(fill = Country), color = NA, show.legend = FALSE) +
scale_fill_manual(values = col)