空間分析　 實習一

助教 杜承軒 2021.03.06

建立繪製地圖的函數 Pollution_Map(arg1)；引數arg1是可自行設定的超越機率(e.g. 0.2)
1. 該函數會回傳該超越機率所對應的PSI值。
2. 以此數值為臨界值，繪製空氣污染地圖，超過該數值的測站，表示紅色，其餘為藍色。
3. 以此數值為臨界值，針對超過該數值的測站，按照測站類別(SiteType)，依照「一般測站、工業測站、交通測站」這三類，以box plot呈現PSI分布。

---

＊初步環境建置與讀取檔案

library(sf);library(tmap);library(pals);library(cartography);library(grid);library(ggplot2);library(GISTools)
setwd("D:/1092SA/Data")
EPA=st_read("EPA_STN1.shp", options="ENCODING=BIG5")
POP=st_read("Popn_TWN2.shp", options="ENCODING=BIG5")
windowsFonts(TOP=windowsFont("Topedia Sans TW Beta"))
windowsFonts(JF=windowsFont("jf金萱 半糖"))

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(POP,xlim=c(146500,351000),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.2)
 map = map + tm_shape(lowEPA) + tm_dots(col="blue", size= 0.2)
 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

---

Pollution_Map2=function(arg1){
 # (1)回傳超越機率對應的PSI值。
 ind=qnorm(arg1,59,17.4,F)
 
 # (2) 繪製空氣污染地圖。
 highEPA=EPA[EPA$PSI>ind,]
 lowEPA=EPA[!EPA$PSI>ind,]

 par(mfrow=c(1,2),family ="TOP")
 par(mar = c(1,0,1.5,0))
 plot(st_geometry(EPA),cex=0, main = "台灣空氣汙染地圖") #創一個範圍邊界
 plot(st_geometry(POP),col='green4',add=T,border='grey20')
 plot(st_geometry(highEPA),col='red',add=T,pch=20)
 plot(st_geometry(lowEPA),col='blue',add=T,pch=20)
 map.scale(140000, 2490000,50000, "公里",2,25)
 north.arrow(130000,2770000,7000,lab='N',col='black')
 legend("bottomleft", legend=c("高於臨界值", "低於臨界值"),col=c("red", "blue"),pch =16,bty='n')
 
 # (3) box plot呈現PSI分布。
 par(mar = c(2.5,2,1.5,1))
 attach(highEPA)
 boxplot(PSI[SiteType=="一般測站"],PSI[SiteType=="工業測站"],PSI[SiteType=="交通測站"],names=c("一般測站","工業測站","交通測站"),main="高於臨界值的PSI盒狀圖")
 detach(highEPA)
 par(mfrow=c(1,1))
 return(ind)
}

---

Pollution_Map2(0.3)

## [1] 68.12457

---

Pollution_Map2(0.5)

## [1] 59