visual_7_warmup

Spatiotemperal R09228001 YANG 楊宇翔

Data: EPA_daily.csv Comparing county-level monthly air pollution in time and space from Taiwan EPA monitoring data Original scale: X-Y point (space); daily (time) Upscaling to county-level (space) and monthly (time)
Plotting [1] Using stplot() for plotting county-level monthly air pollution in time and space [2] Creating boxplot for comparing air pollution among counties

load in file

rm(list=ls())
#install.packages("spacetime")
#install.packages("maps")
library(sp)
library(spacetime)
library(maps)
library(maptools)
library(dplyr)
library(sf)
library(reshape2) 
library(ggplot2)
library(lubridate)
setwd("~/Desktop/110-1/110-1 data visualization/w7_spatiotemperal/ST_data")

TW = st_read('taiwan.shp',options= "ENCODING=utf8",quiet=T) 
TW = subset(TW,!TW$COUNTYENG %in% c("Kinmen County","Lienchiang County","Penghu County"))
TW = TW %>% group_by(COUNTYENG) %>% summarise()

TW_s = as(TW, "Spatial")

TW_s$COUNTYENG

##  [1] "Changhua County" "Chiayi City"     "Chiayi County"   "Hsinchu City"   
##  [5] "Hsinchu County"  "Hualien County"  "Kaohsiung City"  "Keelung City"   
##  [9] "Miaoli County"   "Nantou County"   "New Taipei City" "Pingtung County"
## [13] "Taichung City"   "Tainan City"     "Taipei City"     "Taitung County" 
## [17] "Taoyuan City"    "Yilan County"    "Yunlin County"

setwd("~/Desktop/110-1/110-1 data visualization/w7_spatiotemperal/ST_data")
library(readxl)
data =read_excel("EPA_daily.xlsx")
time = as.POSIXct(data$Time, tz= "GMT")
data = melt(data,id=c("Index","Time","ID"))
EPA = read_sf("EPA_STN1.shp",options= "ENCODING=Big5")
EPA_sp = SpatialPoints(st_coordinates(EPA))
data_st = STFDF(EPA_sp,time,data[order(data[,2]),])
proj4string(data_st) = '+init=epsg:3826'


agr_stf = STF(TW_s,as.Date(c("2021-11-01", "2021-12-01", "2022-01-01","2022-02-01"))) 
proj4string(agr_stf) = '+init=epsg:3826'
arg = aggregate(data_st[,,"value"],agr_stf, mean, na.rm=TRUE)


for( i in 1:19){
  arg@sp@polygons[[i]]@ID = TW$COUNTYENG[i]
}
library(RColorBrewer)
stplot(arg,col.regions= brewer.pal(9, "GnBu"), cuts = 9, sp.layout = list(EPA_sp, pch=16, cex=0.4, col='black'),main = "spatial temperal map of CO concentration")

#使用melt()，把寬表（wide table）轉成長表（long table） # 刪除離島 # 加總到縣市

par(family)

## NULL

stplot(arg, mode = 'tp',ylab='Concentration',main='Concentraion of CO ')

stf_tw_day <- STF(TW_s, 
                  seq.Date(from = as.Date("2021/11/01", format = "%Y/%m/%d"), by = "day", length.out = 120)) 
proj4string(stf_tw_day)<- proj4string(data_st)
arg_2 <- aggregate(data_st[,,"value"], stf_tw_day, mean, na.rm=TRUE)

for( i in 1:19){
   arg_2@sp@polygons[[i]]@ID = TW$COUNTYENG[i]
}


sp_2 <- as(arg_2, "Spatial")
sp_2_data <- sp_2@data 

sp_2_data <- sp_2_data %>% 
  mutate(county = rownames(sp_2@data))

sp_2_data_melt <- sp_2_data %>%  melt(id = c("county"))  # 改為long table

ggplot(sp_2_data_melt, 
       aes(x = county, y = value, fill = county)) +
       geom_boxplot() +  
       xlab("County") + 
       ylab("Air pollution")+
  theme(text = element_text(family = "Heiti TC Light"))

The End