空間分析 作業一

園藝碩一 R10628306 趙旻誼

*讀檔與環境建置

library(sf)
library(tmap)
library(pals)
library(cartography)
library(grid)
library(ggplot2)
library(GISTools)

setwd("D:/110SA/data")
EPA=st_read("EPA_STN1.shp", options="ENCODING=BIG5")
Popn=st_read("Popn_TWN2.shp", options="ENCODING=BIG5")
windowsFonts(TOP=windowsFont("Times New Roman"))
windowsFonts(JF=windowsFont("凝書體 1.0 Regular"))

Pollution_Map=function(arg1){
  #(1)回傳超越機率對應的PSI值
  ind=qnorm(arg1,59,17.4,F)
  
  #(2)繪製空氣汙染地圖
  highEPA=EPA[EPA$PSI>ind,]
  lowEPA=EPA[!EPA$PSI>ind,]
  
  map = tm_shape(Popn,xlim=c(146500,3510000),ylim=c(2400000,2850000))+
    tm_polygons(col='green4',border.col = 'grey20',lwd=0.1)+
    tm_scale_bar(width = 0.3)+
    tm_compass(position = c(.05,.8))+
    tm_layout(frame=F,title="臺灣空氣汙染地圖",title.size = 1,title.position = c("center","top"),fontfamily="JF")
  map = map + tm_shape(highEPA) + tm_dots(col="red", size = 0.1)
  map = map + tm_shape(lowEPA) + tm_dots(col="blue", size = 0.1)
  map = map + tm_add_legend("symbol",labels=c("高於臨界值","低於臨界值"),col=c("red","blue"))
  
  #(3) box plot 呈現PSI分布
  highEPA=highEPA[highEPA$SiteType%in%c("一般測站","交通測站","工業測站"),]
  boxplot = ggplot(highEPA,aes(x=SiteType, y= PSI)) +
   geom_boxplot()+
   ggtitle("高於臨界值的PSI盒狀圖")+
   xlab("測站類別")+theme_minimal()+
   theme(plot.margin=unit(c(1,1,1,1),"cm"),plot.title = element_text(hjust = 0.5),text=element_text(family="JF"))

#繪圖
 grid.newpage()
  pushViewport(viewport(layout=grid.layout(1,2)))
 print(map, vp=viewport(layout.pos.col = 1))
 print(boxplot, vp=viewport(layout.pos.col = 2))
 
 return(ind)
}

Pollution_Map(0.3)

## [1] 68.12457

Pollution_Map(0.5)

## [1] 59