SIEW (2021) On the spatial concentration of economic activity in East Asia and Pacific: New evidence from satellite night-time lights
This is a computational notebook to replicate the results of the Directional LISAs for the above-mentioned paper, to execute, use this link: https://deepnote.com/project/Directional-LISAs-KZtSSOQGTQSqkio1uksneg/%2Fnotebook.ipynb ## please be aware that due to a large sample size, and programming in a language not native to Deepnote, it might take around 30 minutes of execution time ##
Install package
devtools::install_github("amvallone/estdaR")
Skipping install of 'estdaR' from a github remote, the SHA1 (02d72ccd) has not changed since last install.
Use `force = TRUE` to force installation
Load libraries
library(estdaR)
library(sf)
library(tidyverse)
library(spdep)
library(ggplot2)
library(RColorBrewer)
library(ggrepel)
library(skimr)
Loading required package: spdep
Loading required package: sp
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')`
Loading required package: sf
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Loading required package: ggplot2
Attaching package: ‘estdaR’
The following objects are masked from ‘package:spdep’:
geary, moran
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
✔ tibble 3.1.0 ✔ dplyr 1.0.4
✔ tidyr 1.1.2 ✔ stringr 1.4.0
✔ readr 1.4.0 ✔ forcats 0.5.1
✔ purrr 0.3.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
Load shapefile
sh <- st_read("satellite.gpkg")
shdf = read_csv("satellite.csv")
skim(shdf)
Reading layer `satellite' from data source `/work/satellite.gpkg' using driver `GPKG'
Simple feature collection with 69000 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 73.5577 ymin: -55.11694 xmax: 188.6023 ymax: 53.56086
Geodetic CRS: WGS 84
Warning message:
“Missing column names filled in: 'X1' [1]”
── Column specification ────────────────────────────────────────────────────────
cols(
X1 = col_double(),
id = col_double(),
Country = col_character(),
Region = col_character(),
Sat_2012 = col_double(),
Sat_2018 = col_double(),
Pop_2010 = col_double(),
Pop_2020 = col_double(),
X = col_double(),
Y = col_double()
)
── Data Summary ────────────────────────
Values
Name shdf
Number of rows 69000
Number of columns 10
_______________________
Column type frequency:
character 2
numeric 8
________________________
Group variables None
── Variable type: character ────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max empty n_unique whitespace
1 Country 0 1 4 16 0 25 0
2 Region 266 0.996 1 39 0 47346 0
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd p0 p25
1 X1 0 1 34500. 19919. 1 17251.
2 id 0 1 34500. 19919. 1 17251.
3 Sat_2012 0 1 2.90 9.40 -0.185 0.155
4 Sat_2018 0 1 4.13 14.9 0.209 0.413
5 Pop_2010 10 1.00 2440. 8938. 0 131.
6 Pop_2020 10 1.00 2620. 9319. 0 133.
7 X 0 1 119. 11.6 74.9 107.
8 Y 0 1 13.3 10.9 -54.6 10.1
p50 p75 p100 hist
1 34500. 51750. 69000 ▇▇▇▇▇
2 34500. 51750. 69000 ▇▇▇▇▇
3 0.279 0.906 304. ▇▁▁▁▁
4 0.644 1.80 843. ▇▁▁▁▁
5 347. 963. 280868. ▇▁▁▁▁
6 370. 1051. 447646. ▇▁▁▁▁
7 121. 124. 189. ▁▅▇▁▁
8 13.6 16.8 52.9 ▁▁▂▇▁
Compute spatial weights
coord = data.frame(sh$X,sh$Y)
knn <- knn2nb(knearneigh(as.matrix(coord), k = 8))
w <- nb2listw(knn, style = "W",zero.policy=TRUE)
Standardised variables
sat12r <- sh$Sat_2012 - mean(sh$Sat_2012)
sat18r <- sh$Sat_2018 - mean(sh$Sat_2018)
1) VIIRS satellite night-time lights radiance intensity
Compute directional LISAs (no arrows)
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).”
Compute directional LISAs (arrow plot)
a<-dl$lisamap
b<-a$data
b<-as.data.frame(b)
dl_plot <-
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"))) +
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(dl_plot, file="Directional LISAs plot.png")
Saving 6.67 x 6.67 in image
dl_plot
Rose diagram
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")
rose_shdi <-
ggplot(b, aes(x=angle, y=length, col=Regime)) +
geom_col(width = 47) +
coord_polar() +
labs(title = "Rose Diagram", color="Country") +
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(rose_shdi, file="Directional LISAs rose.png")
Saving 6.67 x 6.67 in image
Warning message:
“position_stack requires non-overlapping x intervals”
rose_shdi
Warning message:
“position_stack requires non-overlapping x intervals”
P-value
a<- dl$p.value
palette <- c('#f1eef6','#bdc9e1','#74a9cf','#0570b0')
p_shdi <-
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(p_shdi, file="Directional LISAs p-value.png")
Saving 6.67 x 6.67 in image
p_shdi