In-class Exercise 2: EHSA

Published

November 25, 2023

Modified

November 30, 2023

Overview

Getting Started - Import packages

This function calls pacman to load sf, tidyverse, tmap, knitr packages;

  • tmap : For thematic mapping; powerful mapping package
  • sf : for geospatial data handling, but also geoprocessing: buffer, point-in-polygon count, etc
    • batch processing over GIS packages; can handle tibble format
  • sfdep : creates space-time cube, EHSA; replaces spdep
  • tidyverse : for non-spatial data handling; commonly used R package
  • plotly : makes R graphs interactive; zoom, onmouseover details
    • juse ggplotly(p) where p is a ggplot object
show code 
pacman::p_load(tmap, sf, sfdep, tidyverse, plotly, zoo, Kendall)

Loading the data

  • Hunan: geospatial dataset in ESRI shapefile format
    • use of st_read() to import assf data.frame
      • $geometry column is actually a list inside the df cell; that’s the power of the tibble dataframe
      • “features” of simple features refers to geometric features eg point line curve etc
    • note projection is WGS84; see `88
  • Hunan_GDPPC.csv: attribute format in csv format
    • unlike previous csv, this is a time-series data; columns are YEAR / COUNTY / GDPPC
    • !NOTE! year needs to be numerical and sequential; no datestamp no character not time object
  • !IMPORTANT! to retain geometry, you must left join to the sf dataframe (eg you can also hunan2012 right join hunan)
    • without sf dataframe, normal tibble dataframe will drop the geometry column
show code 
hunan <- st_read(dsn = "data/geospatial", 
                 layer = "Hunan")
Reading layer `Hunan' from data source 
  `C:\1darren\ISSS624\In-class_Ex\In-class_Ex2\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 88 features and 7 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 108.7831 ymin: 24.6342 xmax: 114.2544 ymax: 30.12812
Geodetic CRS:  WGS 84
show code 
hunan2012 <- read_csv("data/aspatial/Hunan_2012.csv")
hunan_GDPPC <- left_join(hunan,hunan2012)%>%
  select(1:4, 7, 15)
GDPPC <-read_csv("data/aspatial/Hunan_GDPPC.csv")

Create spacetime cube

  • Use the spacetime() function;
    • Specify space/time columns using .loc_col, .time_col _ Use is_spacetime_cube() to confirm operation performed successfully
    • GDPPC and GDPPC_st object look very similar, oclumns etc;
show code 
GDPPC_st <- spacetime(GDPPC, hunan, 
                      .loc_col = "County",
                      .time_col = "Year"
                      )
is_spacetime_cube(GDPPC_st)
[1] TRUE

Perform Gi* analysis

  • Gi* needed for
  • activate on each time-period, reconsider the geometry column for future calculations
    • each row for nb, wt has redundant info; however, this will be needed for Gi* calculations over time
    • more tedious, but
  • include_self into nb matrix to
show code 
GDPPC_nb <- GDPPC_st %>%
  activate("geometry") %>%
  mutate(nb = include_self(st_contiguity(geometry)), 
         wt = st_inverse_distance(nb, geometry, scale = 1, alpha = 1),
         .before = 1) %>%
  set_nbs("nb") %>%
  set_wts("wt")
  • use group_by(Year) to analyse year by year
  • use tidyr::unnest to expand output of local_gstar_perm function
show code 
gi_stars <- GDPPC_nb %>%
  group_by(Year) %>%
  mutate(gi_star = local_gstar_perm(
    GDPPC, nb, wt)) %>%
  tidyr::unnest(gi_star)

Perform Mann-Kendall

  • Prof you need to show us the cbg section
show code 
# cbg <- gi_stars %>%
  
  
# ggplot(data = cbg,
#        ) +
#   geom_line() +
#   theme_light()
# 
# 
# ehsa <- emerging_hotspot_analysis(
#   x = GDPPC_st,
#   .var = "GDPPC",
#   k = 1,
#   nsim = 99
# )

EHSA

  • ehsa has 88 rows; 10 years’ data condensed down by individual region

    • ehsa$classification shows classification of region into sporadic coldspot, oscillating hotspot etc

  • after plot, some have “no pattern”

  • note that NO PATTERN and NO VALUE is different;

    • no pattern has no pattern in time/space
    • no value is grey – not statistically significant, e.g.
show code 
ehsa <- emerging_hotspot_analysis(
  x = GDPPC_st,
  .var = "GDPPC",
  k = 1,
  nsim = 99
)
hunan_ehsa <- left_join(hunan_GDPPC,ehsa, by = join_by("NAME_2" == "location"))
ehsa_sig <- hunan_ehsa %>%
  filter(p_value < 0.05)

tmap_mode("plot")
tm_shape(hunan_ehsa) + 
  tm_polygons() +
  tm_borders(alpha = 0.5) + 
tm_shape(ehsa_sig) + 
  tm_fill("classification") + 
  tm_borders(alpha = 0.4)