以R語言爬取監測站歷史資料並以ggplot2繪製風玫瑰圖(風花圖,Wind Rose)_大寮測站為例

--- 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] ![](https://i.imgur.com/WRfddTo.jpg) :::danger 2024/01/20更新非常可惜的是以前的網站更新了QQ 沒辦法在像之前一樣更改網址爬取資料目前解決方式只好先分成兩步驟先使用Python的Selenium下載CSV檔再用ggplot2畫圖相關程式碼我就分享出來拋磚引玉看看有沒有大神能夠將整個程式碼在優化一下或是有更好用的爬蟲套件能分享一下不然光Selenium的安裝跟使用就耗費我好多時間XD Chrome的版本太新就是一個大問題chromedriver支援還沒跟上QQ 要想像之前架設一個Shiny能自動爬取跟畫圖看來是有困難QQ [使用Python selenium 爬取codis氣候資料服務系統資料並使用ggplot2繪製風花圖](/LW3Ys5qwQ4a4j6zrOdc8HA) ::: 今天來分享一個環工系比較可能用到的功能~~ 或是還有哪些科系會用到我就不清楚了下面網站可以查詢台灣的監測資料我們可以在上面查到許多台灣的降雨、風向、風速、日照等資料 ![](https://i.imgur.com/cpWgicg.jpg) 這次以公司附近的大寮測站為例之前也有一篇相關文章大家可以看看[使用R(RStudio Cloud)語言擷取觀測資料查詢系統氣象資料](/cAVB2a4XRsC4SvAFhmygNQ) ![](https://i.imgur.com/zLH9aIJ.jpg) ![](https://i.imgur.com/1W1mcQZ.png) ```r= 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/ 有空來調整看看 ::: ```r= #貼上網址刪除最後面的日期 #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("") } ``` ![](https://i.imgur.com/4h3OQhE.png) # Shiny APP ```r= # # 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()中再判斷選擇的是什麼而執行哪個函數 ::: ![](https://i.imgur.com/gnRUNBV.png) # RStudio Cloud 專案 https://rstudio.cloud/spaces/104923/project/2596118 # [Shinydashboard(儀錶板)手動添加數據繪製風花圖](/Kd2uUuldTrmfafDFt_GfCA) ![](https://i.imgur.com/KFdMCHb.png) 🌟全文可以至下方連結觀看或是補充全文分享至 https://www.facebook.com/LHB0222/ https://www.instagram.com/ahb0222/ 有疑問想討論的都歡迎於下方留言喜歡的幫我分享給所有的朋友 \o/ 有所錯誤歡迎指教 # [:page_with_curl: 全部文章列表](https://hackmd.io/@LHB-0222/AllWritings) ![](https://i.imgur.com/47HlvGH.png)

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.