用R做氣候圖

tags: 用R做, 學習筆記

氣候圖（climograph）是一種用來呈現某地長時間氣候的圖表，主要內容就是該地的氣溫與降水狀況，一般會呈現月均溫和各月雨量的分布。透過氣候圖，我們可以更直觀地了解一個地方的氣候型態。
除了單純的氣候圖之外，也有一種在植群生態學中使用的生態氣候圖（climate diagram），藉由特定的氣溫與降水量比例（莽原氣候

${10}^{\circ }C=30mm$、地中海型氣候

${10}^{\circ }C=20mm$），可以在氣候圖中呈現降水量與蒸散量的關係（若該月降水量超過 100mm 則為重濕），以及水的型態可否為植物所利用（是否低於

$0^{\circ }C$）等等。不過比較常看到的還是普通的氣候圖。

---

以下以中央氣象局香山觀測站（編號 C0D570，目前已撤銷）為例，用 R 的幾個套件從測站觀測資訊網把觀測資料抓下來，再繪製成氣候圖。
先來抓資料，把香山測站自設站到撤銷的所有完整年度資料抓下來：
（這部分的程式碼修改自阿好伯用來抓取風速風向資料的程式碼，感謝阿好伯。）

packages <- c("dplyr","jsonlite", "rvest", "magrittr", "lubridate","ggplot2")
invisible(lapply(packages, library, character.only = TRUE))

#載入資料，網址隨各站不同
url_start <- "https://e-service.cwb.gov.tw/HistoryDataQuery/YearDataController.do?command=viewMain&station=C0D570&stname=%25E9%25A6%2599%25E5%25B1%25B1&datepicker="

climatecrawler <- 
  function(x,y){
    begin_year <- x
    end_year <- y
    data1 <- data.frame()
    data2 <- data.frame()
    for(r in c(0:(end_year - begin_year))){
      url <- paste0(url_start,as.character(begin_year + r),"&altitude=15m") #要記得加上測站海拔
      #時間資料
      month <- url %>% 
        read_html() %>%
        html_nodes(xpath='//*[(@id = "MyTable")]//td[(((count(preceding-sibling::*) + 1) = 1) and parent::*)]') %>% 
        html_text(trim	= T) %>% 
        as.numeric()
      #均溫資料
      avgtemp <- url %>%
        read_html() %>%
        html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 8) and parent::*)]') %>%
        html_text(trim	= T) %>% 
        gsub("X","-9999", .) %>% 
        as.numeric()
      #高溫資料
      hightemp <- url %>%
        read_html() %>%
        html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 9) and parent::*)]') %>%
        html_text(trim	= T) %>% 
        gsub("X","-9999", .) %>% 
        as.numeric()
      #低溫資料
      lowtemp <- url %>%
        read_html() %>%
        html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 11) and parent::*)]') %>%
        html_text(trim	= T) %>% 
        gsub("X","-9999", .) %>% 
        as.numeric()
      #降水量
      precip <- url %>%
        read_html() %>%
        html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 19) and parent::*)]') %>%
        html_text(trim	= T) %>% 
        gsub("X","-9999", .) %>% 
        as.numeric()
      #日照時數
      sunshine <- url %>%
        read_html() %>%
        html_nodes(xpath='//td[(((count(preceding-sibling::*) + 1) = 30) and parent::*)]') %>%
        html_text(trim	= T) %>% 
        gsub("X","-9999", .) %>% 
        as.numeric()
      #合併資料
      data1 <- cbind(year = as.character(begin_year + r), month, avgtemp, hightemp, lowtemp, precip, sunshine)
      data2 <- rbind(data2, data1[-1,])
      data3 <- data2[order(as.numeric(data2$month)),]
    }
    return(data3)
  }

SSweather <- climatecrawler(2008, 2021)
write.csv(SSweather, "香山測站2008_2021觀測資料.csv")

接著，用抓到的資料作圖：

library(dplyr)
weather <- read.csv("香山測站2008_2021觀測資料.csv", header = T)
total <- as.data.frame(weather %>%
                         group_by(month) %>%
                         summarise(
                           temp = mean(avgtemp, na.rm = T),
                           hightemp = mean(hightemp, na.rm = T),
                           lowtemp = mean(lowtemp, na.rm = T),
                           precip = mean(precip, na.rm = T)
                         )
                       )

library(scales)
library(ggpubr)
library(ggthemes)
windowsFonts(GJ = windowsFont("Gen Jyuu Gothic P Medium"))

TempPrecp <- ggplot(total) +
  geom_bar(aes(y = precip *20/300, x = month), stat="identity", fill = "#19B4C5", colour = "#19B4C5") +
  geom_line(mapping = aes(y = temp, x = month), size = 1.5, color = "#F9BD21") + 
  geom_point(mapping = aes(y = temp, x = month), size = 4, color = "#F9BD21") + 
  #geom_ribbon(aes(x = month, y = temp,
  #                ymin = hightemp,
  #                ymax = lowtemp,
  #                fill = "#F9BD21"), alpha = 0.2) +
  scale_x_continuous(name = "月份", breaks = 1:12) +
  #scale_y_continuous(name = "月均溫 (℃) \n", breaks = seq(0, 30, by = 5),
  #                   sec.axis = sec_axis(~. *300/30, name = "降水量(mm) \n")) +
  scale_y_continuous(name = "月均溫 (℃) ", limit = c(0, 35), breaks = c(10, 20, 30),
                     sec.axis = sec_axis(~. *300/20, name = "降水量(mm) \n")) +
  theme_few() +
  ggtitle("香山氣象站氣候圖(2008–2021年)") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 32),
        plot.margin = unit(c(1,2.5,1,1),'cm'),
        legend.position="none") +
  font("xy.text", size = 18) 
ggpar(TempPrecp, font.family = "GJ", font.x = c(28, "azure4"), font.y = c(28, "azure4"))

最後輸出如下：

繪製出的氣候圖。

兩側的 Y 軸單位與數字不同，用來繪圖的溫度和雨量數字大小也不同，因此需要乘以最高月均溫與最大月降水量的比值。
此處繪圖順序是將降水量作為第 2 個 Y 軸，因此畫長條圖時要乘以

$(\frac{\mathrm{最}\mathrm{高}\mathrm{月}\mathrm{均}\mathrm{溫}}{\mathrm{最}\mathrm{大}\mathrm{月}\mathrm{降}\mathrm{水}\mathrm{量}})$，而在繪製軸線刻度時則相反，要將主軸線刻度（記為 ~.）乘以

$(\frac{\mathrm{最}\mathrm{大}\mathrm{月}\mathrm{降}\mathrm{水}\mathrm{量}}{\mathrm{最}\mathrm{高}\mathrm{月}\mathrm{均}\mathrm{溫}})$。

除了氣候圖之外，當然也可以只畫溫度或雨量圖：

total2 <- as.data.frame(weather %>%
                         group_by(year) %>%
                         summarise(
                           temp = mean(avgtemp, na.rm = T),
                           hightemp = mean(hightemp, na.rm = T),
                           lowtemp = mean(lowtemp, na.rm = T),
                           precip = sum(precip, na.rm = T)
                         )
)
tem <- ggline(total2, "year", "temp", size = 1, point.size = 1.5, color = "#7F1084") +
  labs(x = "\n 年度", y = "年均溫 (℃) \n", title = "香山氣象站 2008–2021年均溫變化") +
  scale_x_continuous(breaks = seq(2008, 2021, by = 1)) +
  scale_y_continuous(limits = c(22, 23.5), breaks = c(22.0, 22.25, 22.5, 22.75, 23.0, 23.25, 23.5)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 32),
        plot.margin = unit(c(1,2.5,1,1),'cm')) +
  font("xy.text", size = 18) +
  theme_hc()
ggpar(tem, font.family = "GJ", font.x = c(28, "azure4"), font.y = c(28, "azure4")) + 
  geom_hline(yintercept = 22.6583, lty = 5, color = "black")

繪製出的年均溫變化圖。

可以看到香山氣象站這 14 年來的的年均溫趨勢往上，真是讓人擔心啊。

而最開頭提到的生態氣候圖，可以參考嘉義大學林政道老師的程式碼。
🐕‍🦺 2022.10.06