使用R語言進行資料分析期末專案

資料來源

  1. EPA 測站
  2. CI Air box (民生公共物聯網資料平台)
  3. cems 固定排放源連續監測資料

目標

  1. 互動式即時地圖
  2. 目標汙染事件分析
  3. 嚴重汙染地區分析

概念

  1. function
  2. tidyverse
  3. leaflet . ggmap

install packages

中研院空氣盒子GPS分布

library(jsonlite)
library(dplyr)
GPS_url<-"https://raw.githubusercontent.com/LinkItONEDevGroup/LASS/master/LASS_DB/AirBox/school.json"
GPS_json<-fromJSON(GPS_url)
GPS_df<-do.call("rbind", GPS_json)
str(GPS_df)
## 'data.frame':    1942 obs. of  8 variables:
##  $ id  : chr  "28C2DDDD45AA" "28C2DDDD43ED" "28C2DDDD4408" "28C2DDDD43F8" ...
##  $ area: chr  "taipei" "taipei" "taipei" "taipei" ...
##  $ org : chr  "es" "es" "es" "es" ...
##  $ name: chr  "臺北市關渡國小" "臺北市金華國小" "臺北市民族國小" "臺北市國語實小" ...
##  $ en  : chr  "" "" "" "" ...
##  $ addr: chr  "[112]台北市北投區中央北路四段581號" "[106]臺北市大安區永康里愛國東路79巷11號" "[105]臺北市松山區民生東路4段97巷7號" "[100]臺北市中正區龍光里南海路58號" ...
##  $ lat : num  25.1 25 25.1 25 25 ...
##  $ lon : num  121 122 122 122 122 ...
library(leaflet)
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:ggmap':
## 
##     inset
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
m <- leaflet(GPS_df) %>%
  addTiles() %>%
  setView(lng = 120.239, lat = 22.992, zoom = 8)%>%
  addCircles()%>%
  setView(121.53435,25.01711, zoom = 13)
## Assuming "lon" and "lat" are longitude and latitude, respectively
m

環保署測站分布,即時值

library(leaflet)
library(magrittr)
library(jsonlite)
library(RCurl)
## Loading required package: bitops
## 
## Attaching package: 'RCurl'
## The following object is masked from 'package:tidyr':
## 
##     complete
library(tidyverse)

epa_instant_url<-"https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Datastreams?$expand=Thing,Observations($orderby=phenomenonTime%20desc;$top=1)&$filter=name%20eq%20%27PM2.5%27&$count=true"

epa_latest <-fromJSON(epa_instant_url)
epa.df <- as.data.frame(epa_latest$value)


x<-epa.df$Observations
x_tib <- x %>% 
    enframe() %>% 
    unnest()

#clear the data
epa.clear<-epa.df%>%
  select(Thing,Observations)%>%
  mutate(Name=Thing$name,
         Lat=Thing$properties$TWD97Lat,
         Lon=Thing$properties$TWD97Lon,
         Result=x_tib$result)

head(epa.clear)
##   Thing.name       Thing.description Thing.properties.AreaName
## 1       汐止 環保署空氣品質測站-汐止                北部空品區
## 2       萬華 環保署空氣品質測站-萬華                北部空品區
## 3       桃園 環保署空氣品質測站-桃園                北部空品區
## 4       古亭 環保署空氣品質測站-古亭                北部空品區
## 5       中山 環保署空氣品質測站-中山                北部空品區
## 6       新竹 環保署空氣品質測站-新竹                竹苗空品區
##   Thing.properties.SiteEngName Thing.properties.SiteName
## 1                        Xizhi                      汐止
## 2                       Wanhua                      萬華
## 3                      Taoyuan                      桃園
## 4                       Guting                      古亭
## 5                    Zhongshan                      中山
## 6                      Hsinchu                      新竹
##        Thing.properties.Address Thing.properties.Township
## 1 新北市汐止區樟樹一路137巷26號                    汐止區
## 2     臺北市萬華區中華路1段66號                    萬華區
## 3        桃園市桃園區莒光街15號                    桃園區
## 4  臺北市大安區羅斯福路3段201號                    大安區
## 5     臺北市中山區林森北路511號                    中山區
## 6          新竹市東區民族路33號                      東區
##   Thing.properties.TWD97Lat Thing.properties.SiteType
## 1                  25.06567                  一般測站
## 2                  25.04650                  一般測站
## 3                  24.99567                  一般測站
## 4                  25.02061                  一般測站
## 5                  25.06236                  一般測站
## 6                  24.80562                  一般測站
##   Thing.properties.City Thing.properties.TWD97Lon Thing.@iot.id
## 1                新北市                  121.6408            22
## 2                臺北市                  121.5080            63
## 3                桃園市                  121.3044            43
## 4                臺北市                  121.5296            16
## 5                臺北市                  121.5265            11
## 6                新竹市                  120.9721            56
##                                           Thing.@iot.selfLink
## 1 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(22)
## 2 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(63)
## 3 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(43)
## 4 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(16)
## 5 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(11)
## 6 https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Things(56)
##                                                                                                        Observations
## 1 2019-03-05T12:00:00.000Z, NA, 17, 6448521, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448521)
## 2 2019-03-05T12:00:00.000Z, NA, 19, 6448706, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448706)
## 3 2019-03-05T12:00:00.000Z, NA, 14, 6448770, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448770)
## 4 2019-03-05T12:00:00.000Z, NA, 13, 6448723, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448723)
## 5 2019-03-05T12:00:00.000Z, NA, 15, 6448689, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448689)
## 6 2019-03-05T12:00:00.000Z, NA, 11, 6448884, https://sta.ci.taiwan.gov.tw/STA_AirQuality/v1.0/Observations(6448884)
##   Name      Lat      Lon Result
## 1 汐止 25.06567 121.6408     17
## 2 萬華 25.04650 121.5080     19
## 3 桃園 24.99567 121.3044     14
## 4 古亭 25.02061 121.5296     13
## 5 中山 25.06236 121.5265     15
## 6 新竹 24.80562 120.9721     11

環保署測站即時值呈現

#PM2.5  color function
PM_color <- function(x){
  if(x <= 10){
    return('#00e800') # green
  }
  else if(x > 10 && x <= 35 ){
    return('#CAFF70') # yellow green
  }
  else if(x > 35 && x <=50 ){
    return('#ffff00') # yello
  }
  else if(x > 50 && x <= 65 ){
    return('#ff7e00') # orange
  }
  else if(x > 65 && x <= 70 ){
    return('#ff0000') # red
  }
  else{
    return('#8f3f97') # purple
  }
}

#add  color column
epa.clear$Color <- sapply(epa.clear$Result, PM_color)

epa.clear<-epa.df%>%
  select(Thing,Observations)%>%
  mutate(Name=Thing$name,
         Lat=Thing$properties$TWD97Lat,
         Lon=Thing$properties$TWD97Lon,
         Result=x_tib$result,
         Color=sapply(epa.clear$Result, PM_color))

# draw EPA instant pm2.5 map
epa_map <- leaflet() %>%
  addProviderTiles(provider = "CartoDB") %>%
  addCircleMarkers(lng = epa.clear$Lon,
             lat=epa.clear$Lat,
             color =epa.clear$Color,
             popup = paste0(epa.clear$Name,
                                 "<br/>",
                                 'PM2.5即時濃度:',epa.clear$Result,''))%>%
   setView(lng=121.1,lat=23.8, zoom = 7)
epa_map

環保署歷史值

library(tidyverse)
library(leaflet)
epa_201904 <- read_csv("C:\\Rdata\\EPA_OD_201904.csv")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   SiteName = col_character(),
##   County = col_character(),
##   Pollutant = col_character(),
##   Status = col_character(),
##   PublishTime = col_character()
## )
## See spec(...) for full column specifications.
epa_time <- epa_201904[1:80,]
epa_station <- read_csv("C:\\Rdata\\epa_station.csv")
## Parsed with column specification:
## cols(
##   station_id = col_character(),
##   SiteName = col_character(),
##   SiteEngName = col_character(),
##   AreaName = col_character(),
##   County = col_character(),
##   Township = col_character(),
##   SiteAddress = col_character(),
##   TWD97Lon = col_double(),
##   TWD97Lat = col_double(),
##   SiteType = col_character()
## )
epa_all <-right_join(epa_time,epa_station,by.x = "SiteName", by.y = "SiteName")
## Joining, by = c("SiteName", "County")
epa_all$PM2.5_AVG[is.na(epa_all$PM2.5_AVG)]<-0 #把缺少的值转成0

epa.clear<-epa_all%>%
  mutate(Name=epa_all$SiteName,
         Lat=epa_all$TWD97Lat,
         Lon=epa_all$TWD97Lon,
         Result=epa_all$PM2.5_AVG,
         Color=sapply(epa_all$PM2.5_AVG, PM_color))

epa_map_time <- leaflet() %>%
  addProviderTiles(provider = "CartoDB") %>%
  addCircleMarkers(lng = epa.clear$TWD97Lon,
             lat=epa.clear$TWD97Lat,
             color =epa.clear$Color,
             popup = paste0(epa.clear$SiteName,
                                 "<br/>",
                                 '2019/4/1/00:00:00時PM2.5濃度:',epa.clear$Result,''))%>%
   setView(lng=121.53435,lat=25.01711, zoom = 10)
epa_map_time

問題

  1. ggmap google api
  2. airbox api
  3. timestamp

團隊內部

  1. 約時間
  2. 分工+互相學習