---
disqus: ahb0222
GA : G-VF9ZT413CG
---
# 以R語言爬取監測站歷史資料並以ggplot2繪製風玫瑰圖(風花圖,Wind Rose)_大寮測站為例
> [color=#40f1ef][name=LHB阿好伯, 2021/02/15][:earth_africa:](https://www.facebook.com/LHB0222/)
###### tags: `R` `ggplot2`

[TOC]


:::danger
2024/01/20更新
非常可惜的是以前的網站更新了QQ
沒辦法在像之前一樣更改網址爬取資料
目前解決方式只好先分成兩步驟
先使用Python的Selenium下載CSV檔
再用ggplot2畫圖
相關程式碼我就分享出來
拋磚引玉看看有沒有大神能夠將整個程式碼在優化一下
或是有更好用的爬蟲套件能分享一下
不然光Selenium的安裝跟使用就耗費我好多時間XD
Chrome的版本太新就是一個大問題chromedriver支援還沒跟上QQ
要想像之前架設一個Shiny能自動爬取跟畫圖看來是有困難QQ
[使用Python selenium 爬取codis氣候資料服務系統資料並使用ggplot2繪製風花圖](/LW3Ys5qwQ4a4j6zrOdc8HA)
:::

今天來分享一個環工系比較可能用到的功能~~
或是還有哪些科系會用到我就不清楚了
下面網站可以查詢台灣的監測資料
我們可以在上面查到許多台灣的降雨、風向、風速、日照等資料

這次以公司附近的大寮測站為例
之前也有一篇相關文章大家可以看看[使用R(RStudio Cloud)語言擷取觀測資料查詢系統氣象資料](/cAVB2a4XRsC4SvAFhmygNQ)



packages <- c("jsonlite", "rvest", "magrittr", "lubridate","ggplot2")

installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# 載入所需套件
  invisible(lapply(packages, library, character.only = TRUE))

:::danger
~~複製網址後刪除最後的時間資料
https://e-service.cwb.gov.tw/HistoryDataQuery/D..... &datepicker= 2021-05-01
貼到程式碼中
url_start <- "https://e-service.cwb.gov.tw/HistoryDataQuery/.....&datepicker="
修改需要的時間範圍執行後續程式碼即可
start_date <- ymd("2021-01-01") #起始時間
end_date <- ymd("2021-03-31") #終止時間~~

網站又...又...更新了
https://opendata.cwa.gov.tw/dataset/observation/O-A0001-001
https://codis.cwa.gov.tw/
有空來調整看看

:::

#貼上網址刪除最後面的日期
#https://e-service.cwb.gov.tw/....&datepicker=2021-05-01
#https://e-service.cwb.gov.tw/....&datepicker=
url_start <- "https://e-service.cwb.gov.tw/HistoryDataQuery/DayDataController.do?command=viewMain&station=C0V730&stname=%25E5%25A4%25A7%25E5%25AF%25AE&datepicker="
start_date <- ymd("2021-01-01")  #起始時間
end_date <- ymd("2021-03-31")  #終止時間

{
data1 <- data.frame()
data2 <- data.frame()
for(d in c(0:(end_date - start_date))){
  url <- paste0(url_start,as.character(start_date + d),"#")
  #時間資料
  time_H <- url %>% 
    read_html() %>%
    html_nodes(xpath='//*[(@id = "MyTable")]//td[(((count(preceding-sibling::*) + 1) = 1) and parent::*)]') %>% 
    html_text(trim	= T) %>% 
    as.numeric()
  #風速資料
  Windspeed <- url %>%
    read_html() %>%
    html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 7) and parent::*)]') %>%
    html_text(trim	= T) %>% 
    gsub("X","-9999", .) %>% 
    as.numeric()
  #風向
  windDirection<- url %>%
    read_html() %>%
    html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 8) and parent::*)]') %>%
    html_text(trim	= T) %>% 
    gsub("X","-9999", .) %>% 
    as.numeric()
  data1 <- cbind(as.character(start_date + d), time_H, Windspeed,windDirection)  #合併資料
  data2 <- rbind(data2, data1)
}
data3 <- na.omit(data2)  #刪除包含NA值的列
data3$windDirection_N <- cut(as.numeric(data3$windDirection),
            breaks = c(0, 11.26, 33.76, 56.26, 78.76, 101.26, 127.76, 146.26, 168.76, 191.26, 213.76, 236.26, 258.76, 281.26, 303.76, 326.26, 348.75, 360),
            labels = c("N","NNE","NE","ENE","E","ESE","ES","SSE","S","SWS","SW","WSW","W","WNW","NW","NNW","N"),
           #   labels = c("北","北東北","東北","東東北","東","東東南","東南","南東南","南","南西南","西南","西西南","西","西西北","西北","西北西","北"),
            include.lowest = TRUE)
data3$Windspeed_N <- cut(as.numeric(data3$Windspeed),
                         breaks =c(0, 2.1, 4.1,  6.1, Inf),
                         labels = c("0~2", "2.1~4", "4.1~6", ">6.1"),
                         right = F)
data3$Windspeed <- as.numeric(data3$Windspeed)

ggplot(data = data3, aes(x = windDirection_N, fill = Windspeed_N))+
  geom_bar( position = position_stack(reverse = TRUE)) + 
  theme_bw() +
  scale_fill_discrete(guide = guide_legend(reverse=TRUE), name = "Wind Speed (m/s)")+
  coord_polar(start = -0.3) + 
  xlab("") 
  }


# Shiny APP

#
# This is a Shiny web application. You can run the application by clicking
# the 'Run App' button above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#


packages <- c("jsonlite", "rvest", "magrittr", "lubridate","ggplot2","shiny")

#install.packages()
library("shiny")
library("jsonlite")
library("rvest")
library("magrittr")
library("lubridate")
library("ggplot2")
library("shiny")
WindData <- function(station_number,start_date,end_date){
  url_start <- paste0("https://e-service.cwb.gov.tw/HistoryDataQuery/DayDataController.do?command=viewMain&station=",station_number,"&stname=%25E5%25A4%25A7%25E5%25AF%25AE&datepicker=")
  start_date <- ymd(start_date)  #起始時間
  end_date <- ymd(end_date)  #終止時間
  data1 <- data.frame()
  data2 <- data.frame()
  for(d in c(0:(end_date - start_date))){
    url <- paste0(url_start,as.character(start_date + d),"#")
    #時間資料
    time_H <- url %>% 
      read_html() %>%
      html_nodes(xpath='//*[(@id = "MyTable")]//td[(((count(preceding-sibling::*) + 1) = 1) and parent::*)]') %>% 
      html_text(trim	= T) %>% 
      as.numeric()
    #風速資料
    Windspeed <- url %>%
      read_html() %>%
      html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 7) and parent::*)]') %>%
      html_text(trim	= T) %>% 
      gsub("X","-9999", .) %>% 
      as.numeric()
    #風向
    windDirection<- url %>%
      read_html() %>%
      html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 8) and parent::*)]') %>%
      html_text(trim	= T) %>% 
      gsub("X","-9999", .) %>% 
      as.numeric()
    data1 <- cbind(as.character(start_date + d), time_H, Windspeed,windDirection)  #合併資料
    data2 <- rbind(data2, data1)
  }
  windD <- factor(c("N","NNE","NE","ENE","E","ESE","ES","SSE","S","SWS","SW","WSW","W","WNW","NW","NNW","N"))
  data3 <- na.omit(data2)  #刪除包含NA值的列
  data3$windDirection_N <- cut(as.numeric(data3$windDirection),
                               breaks = c(0, 11.26, 33.76, 56.26, 78.76, 101.26, 127.76, 146.26, 168.76, 191.26, 213.76, 236.26, 258.76, 281.26, 303.76, 326.26, 348.75, 360),
                               windD,
                               include.lowest = TRUE)
  data3$Windspeed_N <- cut(as.numeric(data3$Windspeed),
                           breaks =c(0, 2.1, 4.1,  6.1, Inf),
                           labels = c("0~2", "2.1~4", "4.1~6", ">6.1"),
                           right = F)
  data3$Windspeed <- as.numeric(data3$Windspeed)
  p1 <- ggplot(data = data3, aes(x = windDirection_N, fill = Windspeed_N))+
    geom_bar( position = position_stack(reverse = TRUE)) + 
    theme_bw() + scale_x_discrete(drop = FALSE,
                                  labels = windD) +
    scale_fill_discrete(guide = guide_legend(reverse=TRUE), name = "Wind Speed (m/s)")+
    coord_polar(start = -0.2) + 
    xlab("") + ggtitle(paste(station_number,"(",start_date,"~",end_date,")"))
  p1
  
  
}
WindData2 <- function(station_number,start_date,end_date){
  url_start <- paste0("https://e-service.cwb.gov.tw/HistoryDataQuery/MonthDataController.do?command=viewMain&station=",station_number,"&stname=%25E5%25A4%25A7%25E5%25AF%25AE&datepicker=")
  start_date <- floor_date(ymd(start_date),unit = "month")  #起始時間
  end_date <- floor_date(ymd(end_date),unit = "month")  #終止時間
  int <- interval(start_date, end_date)
  m1 <- time_length(int, "month")
  data1 <- data.frame()
  data2 <- data.frame()
  for(d in c(0:m1)){
    m2 <- start_date + months(d)
    url <- paste0(url_start,as.character(paste0(year(m2),"-",ifelse(month(m2)>=10,
                                                                    month(m2),paste0(0,month(m2))))),"#")
    #時間資料
    time_d <- url %>% 
      read_html() %>%
      html_nodes(xpath='//*[(@id = "MyTable")]//td[(((count(preceding-sibling::*) + 1) = 1) and parent::*)]') %>% 
      html_text(trim	= T) %>% 
      as.numeric()
    #風速資料
    Windspeed <- url %>%
      read_html() %>%
      html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 17) and parent::*)]') %>%
      html_text(trim	= T) %>% 
      gsub("X","-9999", .) %>% 
      as.numeric()
    #風向
    windDirection<- url %>%
      read_html() %>%
      html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 18) and parent::*)]') %>%
      html_text(trim	= T) %>% 
      gsub("X","-9999", .) %>% 
      as.numeric()
    data1 <- cbind(as.character(start_date + months(d)), time_d, Windspeed,windDirection)  #合併資料
    data2 <- rbind(data2, data1)
  }
  data2
  windD <- factor(c("N","NNE","NE","ENE","E","ESE","ES","SSE","S","SWS","SW","WSW","W","WNW","NW","NNW","N"))
  data3 <- na.omit(data2)  #刪除包含NA值的列
  data3$windDirection_N <- cut(as.numeric(data3$windDirection),
                               breaks = c(0, 11.26, 33.76, 56.26, 78.76, 101.26, 127.76, 146.26, 168.76, 191.26, 213.76, 236.26, 258.76, 281.26, 303.76, 326.26, 348.75, 360),
                               windD,
                               include.lowest = TRUE)
  data3$Windspeed_N <- cut(as.numeric(data3$Windspeed),
                           breaks =c(0, 2.1, 4.1,  6.1, Inf),
                           labels = c("0~2", "2.1~4", "4.1~6", ">6.1"),
                           right = F)
  data3$Windspeed <- as.numeric(data3$Windspeed)
  p1 <- ggplot(data = data3, aes(x = windDirection_N, fill = Windspeed_N))+
    geom_bar( position = position_stack(reverse = TRUE)) + 
    theme_bw() + scale_x_discrete(drop = FALSE,
                                  labels = windD) +
    scale_fill_discrete(guide = guide_legend(reverse=TRUE), name = "Wind Speed (m/s)")+
    coord_polar(start = -0.2) + 
    xlab("") + ggtitle(paste(station_number,"(",paste0(year(start_date),"-",month(start_date)),
                             "~",paste0(year(end_date),"-",month(end_date)),")"))
  p1
  
  
}
# Define UI for application that draws a histogram
ui <- fluidPage(
  titlePanel("以R語言爬取監測站歷史資料並以ggplot2繪製風玫瑰圖(風花圖,Wind Rose)"),
  # Default value is the date in client's time zone
  
  
  # Application title
  
  # Sidebar with a slider input for number of bins 
  
  sidebarLayout(textInput ( "station_number" , " 監測站編號_預設大寮測站", value  =  "C0V730" ),
                sidebarPanel(dateRangeInput('dateRange',
                                            label = 'Date range input: yyyy-mm-dd',
                                            start = Sys.Date() -1, end = Sys.Date() -1),
                             submitButton(text = "更新", icon = NULL, width = NULL),
                             downloadButton("downloadData", "Download Plot")
                )), 
  helpText("測站編號查詢","https://e-service.cwb.gov.tw/HistoryDataQuery/"),
  radioButtons( "mode", "日/月資料",choices = c("日","月")),
  # Show a plot of the generated distribution
  mainPanel(
    plotOutput("WindRosePlote", width = "600px", height = "400px")
  )
)

# Define server logic required to draw a histogram
server <- function(input, output) {
  
  output$WindRosePlote <- renderPlot(
    if(input$mode == "日"){
      WindData(input$station_number, input$dateRange[1], input$dateRange[2])
    } else{
      WindData2(input$station_number, input$dateRange[1], input$dateRange[2])
    }
  )
  output$downloadData <- downloadHandler(
    filename = function(){paste("input$downloadData",'.png',sep='')},
    content = function(file){
      ggsave(file,plot=WindData(input$station_number, input$dateRange[1], input$dateRange[2]))
    }
  )
  
}
# Run the application 
shinyApp(ui = ui, server = server)

## [測試網站](https://gtgrthrst.shinyapps.io/Wind_Rose/)
:::success
* 修改風花圖有些角度沒有數據時不顯示的問題
* 解決方式在原本ggplot2的程式碼加上scale_x_discrete(drop = FALSE, labels = windD)
* 增加"更新"按鍵減少原本更換時間就自動跑程式的問題
* 在shiny ui()中增加 submitButton(text = "更新", icon = NULL, width = NULL)
即可新增一個更新按鈕
* 增加日/月資料的選擇，若是需要一整季的資料可以選擇使用"月"資料減少等待時間
* 增加選擇按鈕radioButtons( "mode", "日/月資料",choices = c("日","月"))
寫兩個函數一個為使用日報表一個使用月報表資料
在server()中再判斷選擇的是什麼而執行哪個函數
:::

# RStudio Cloud 專案

https://rstudio.cloud/spaces/104923/project/2596118

# [Shinydashboard(儀錶板)手動添加數據繪製風花圖](/Kd2uUuldTrmfafDFt_GfCA)


🌟全文可以至下方連結觀看或是補充

全文分享至

https://www.facebook.com/LHB0222/

https://www.instagram.com/ahb0222/

有疑問想討論的都歡迎於下方留言

喜歡的幫我分享給所有的朋友 \o/

有所錯誤歡迎指教

# [:page_with_curl: 全部文章列表](https://hackmd.io/@LHB-0222/AllWritings)