Lab 4：描述疾病擴散的時空趨勢

1 Lab 4：描述疾病擴散的時空趨勢

1.1 0. 套件載入

library(sf)
library(aspace)
library(tmap)
library(ggplot2)
library(dplyr)
library(RColorBrewer)

1.2 1. 讀取資料

setwd("D:/shp")

windowsFonts(TOP = windowsFont("Microsoft JhengHei"))

Dengue <- st_read("point_event.shp")

Reading layer `point_event' from data source `D:\shp\point_event.shp' using driver `ESRI Shapefile'
Simple feature collection with 5735 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 160412 ymin: 2427593 xmax: 328891 ymax: 2785873
CRS:           NA

TW <- st_read("Popn_TWN2.shp", options = "ENCODING=BIG5")

options:        ENCODING=BIG5 
Reading layer `Popn_TWN2' from data source `D:\shp\Popn_TWN2.shp' using driver `ESRI Shapefile'
Simple feature collection with 368 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -197572.5 ymin: 2295201 xmax: 606875.6 ymax: 2919551
Projected CRS: TWD97 / TM2 zone 121

st_crs(Dengue) <- st_crs(TW)

1.3 2. Task 1：探索不同時間尺度下的時間趨勢

Dengue$period <- ceiling((Dengue$WEEK + 1) / 8)

ggplot(Dengue) +
  geom_bar(aes(x = period)) +
  xlab("時段") +
  ylab("病例數") +
  scale_x_continuous(breaks = 1:7) +
  theme_minimal() +
  theme(text = element_text(family = "TOP"))

ggplot(Dengue) +
  geom_bar(aes(x = WEEK)) +
  xlab("週次") +
  ylab("病例數") +
  scale_x_continuous(breaks = seq(0, 55, 5)) +
  theme_minimal() +
  theme(text = element_text(family = "TOP"))

1.4 3. Task 2：探索不同時段下的空間趨勢

KH <- TW[TW$COUNTY == "高雄市", ]

IN <- st_contains(st_union(KH), Dengue, sparse = FALSE)
KH_Dengue <- Dengue[IN, ]

KD <- KH_Dengue %>%
  st_drop_geometry()

1.5 4. Standard Distance 標準距離

col <- RColorBrewer::brewer.pal(7, "Reds")

map <- tm_shape(
  KH,
  xlim = c(16, 20) * 10^4,
  ylim = c(249, 252) * 10^4
) +
  tm_borders()

for (i in 1:7) {
  
  KDi <- KD[KD$period == i, c(2, 3)]
  
  if (nrow(KDi) == 0) next
  
  sddatt <- calc_sdd(points = KDi)
  
  CENTRE <- c(
    sddatt$FORPLOTTING$CENTRE.x,
    sddatt$FORPLOTTING$CENTRE.y
  )
  
  CENTRE.sf <- CENTRE %>%
    st_point() %>%
    st_sfc() %>%
    st_sf()
  
  st_crs(CENTRE.sf) <- st_crs(KH)
  
  rad <- sddatt$ATTRIBUTES$SDD.radius
  SDD <- st_buffer(CENTRE.sf, rad)
  
  map <- map +
    tm_shape(CENTRE.sf) +
    tm_symbols(0.5, col[i]) +
    tm_shape(SDD) +
    tm_borders(col[i], 2)
}

map +
  tm_layout(frame = FALSE)

1.6 5. Standard Deviational Ellipse 標準差橢圓

col <- c("red", "orange", "gold", "green4", "deepskyblue", "blue", "purple")

map <- tm_shape(
  KH,
  xlim = c(16, 20) * 10^4,
  ylim = c(249, 252) * 10^4
) +
  tm_borders()

for (i in 1:7) {
  
  KDi <- KD[KD$period == i, c(2, 3)]
  
  if (nrow(KDi) == 0) next
  
  sdeatt <- calc_sde(points = KDi)
  
  CENTRE <- c(
    sdeatt[["FORPLOTTING"]][["CENTRE.x"]],
    sdeatt[["FORPLOTTING"]][["CENTRE.y"]]
  )
  
  CENTRE.sf <- CENTRE %>%
    st_point() %>%
    st_sfc() %>%
    st_sf()
  
  st_crs(CENTRE.sf) <- st_crs(KH)
  
  SDE.data <- data.frame(
    x = sdeatt$LOCATIONS$x,
    y = sdeatt$LOCATIONS$y
  )
  
  SDE.pt.sf <- st_as_sf(SDE.data, coords = c("x", "y"))
  st_crs(SDE.pt.sf) <- st_crs(KH)
  
  SDE.poly <- st_cast(st_combine(SDE.pt.sf), "POLYGON")
  SDE.sf <- st_sf(SDE.poly)
  
  map <- map +
    tm_shape(CENTRE.sf) +
    tm_symbols(0.5, col[i]) +
    tm_shape(SDE.sf) +
    tm_borders(col[i], 2)
}

map +
  tm_layout(frame = FALSE)