使用R語言分析eBird資料庫來畫出鳥類分布地圖

2021 年 3 月 17 日

eBird是全球最大的賞鳥紀錄資料庫，世界各地的每一個賞鳥者，都可以透過網路或手機上的APP將賞鳥時所觀察到鳥類登錄到資料庫內。累積起來的資料量十分驚人，是研究鳥類分布或遷徙等十分重要的資料來源。也讓許多像小雨燕和黃山雀這樣的業餘賞鳥人也能為鳥類研究貢獻一己之力。

eBird同時也讓資料庫開放給全世界的賞鳥者和研究者，除了提供下載自己的賞鳥紀錄作分析參考，也可以向eBird申請其他的資料產品。

eBird所提供的資料產品大致有3種：（詳情可以參考eBird網站說明）

eBird 基礎資料集 (EBD)

eBird 基礎資料集(EBD)提供所有觀察紀錄的原始資料及其關聯資料(meta data)。EBD於每年的3, 6, 9, 12月的15日更新。EBD開放給所有使用者，可以透過線上填寫資料申請後直接下載。

eBird 參考資料集(ERD)

eBird 參考資料集(ERD)每年更新一次且聚焦在北美和南美洲鳥種

eBird 觀察資料集(EOD)

這個資料集統整了全球鳥類的分布和發生，EOD資料集也可以於全球生物多樣性資訊機構 (Global Biodiversity Information Facility, GBIF ）下載，可以想見資料量十分大。

今天我們就試著用台灣地區eBird基礎資料集(EBD)來畫出白頭翁和烏頭翁在台灣的分布範圍吧。

由於EBD資料集資料十分龐大，不建議使用Excel來開啟。要分析EBD的資料集，建議可以使用官方建議：R語言的auk package 。auk提供了許多語法協助分析，也利用UNIX的AWK指令預先將資料集進行篩選，找出我們要的資料再進行分析。

如果是在Windows上使用則需要先安裝Cygwin才能夠使用auk。

接下來我們到R-Studio。先安裝需要用到的packages並設定資料集的位置。

#安裝等一下要用到的packages
install.packages("auk")
install.packages("tidyverse")
install.packages("ggmap")

library(auk) 
library(tidyverse)
ebd_dir <- "/R/eBird/Taiwan/" #這裡請自行設定EBD資料集放置的資料夾
ebdfile <- file.path(ebd_dir, 'ebd_TW_prv_relNov-2020.txt') #填入檔案名稱

接下來把資料夾裡的資料先進行過濾，篩選出需要用到的資料，我們在此抓取白頭翁(Pycnonotus sinensis)和烏頭翁(Pycnonotus taivanus)的資料

ebd_filters <- auk_ebd(ebdfile) %>% 
  auk_species(c("Pycnonotus sinensis", "Pycnonotus taivanus")) %>%  #設定要篩選的物種
  auk_date(c("*-01-01", "*-12-31")) %>%  #設定時間的區間，在此設定是整年。也可以設定某個季節或某個月份
  auk_complete()

ebd_out_file <- "Bubuls.txt" #過濾後資料集儲存的檔名
ebd_filtered <- auk_filter(ebd_filters, file = ebd_out_file) #執行過濾

ebd <- read_ebd(ebd_out_file) #讀入過濾後資料集

#先大致看一下資料集內容
glimpse(ebd) 

#以下是顯示資料集內容, 共有22萬多筆資料
#Rows: 220,773
#Columns: 45
#$ checklist_id                 <chr> "S53406684", "S53408734", "S53403002", "S53415895", "S53410409", "S53412666", "...
#$ global_unique_identifier     <chr> "URN:CornellLabOfOrnithology:EBIRD:OBS718674790", "URN:CornellLabOfOrnithology:...
#$ last_edited_date             <chr> "2019-03-04 23:13:08", "2019-03-05 00:20:58", "2020-12-06 22:48:48.200685", "20...
#$ taxonomic_order              <dbl> 23533, 23533, 23533, 23533, 23533, 23533, 23533, 23533, 23533, 23533, 23533, 23...
#$ category                     <chr> "species", "species", "species", "species", "species", "species", "species", "s...
#$ common_name                  <chr> "Light-vented Bulbul", "Light-vented Bulbul", "Light-vented Bulbul", "Light-ven...
#$ scientific_name              <chr> "Pycnonotus sinensis", "Pycnonotus sinensis", "Pycnonotus sinensis", "Pycnonotu...

接下來我們要每一筆白頭翁和烏頭翁紀錄的座標，在地圖上標上不同顏色，這樣就可以了解白頭翁和烏頭翁兩種鳥類在台灣分布區域情況。
我們使用ggmap這個package來產生台灣地區的底圖。先設定好地圖四個邊界的經緯度，放大倍率，和地圖型態，在此我們選擇的是地形圖 (terrain)

library(ggmap)
#設定地圖四個邊界的經緯度
tw_borders <- c(bottom  = 21.5, 
                top     = 25.5,
                left    = 119,
                right   = 123)
tw_map <- get_stamenmap(tw_borders, zoom = 8, maptype = "terrain")  
ggmap(tw_map)

這樣一來，我們就可以產生一張我們需要的台灣地圖底圖。接下來，我們只要將每一筆資料點上座標就可以囉。

ggmap所畫的圖和ggplot2的語法物件都能相容，因此，再加上一個圖層並點上我們要的資料

ebd_select <- ebd %>% select(common_name, longitude, latitude) #篩選出鳥種名、經緯度座標
colnames(ebd_select)<-c('Species', 'lon', 'lat')

bubulmap<-ggmap(tw_map, base_layer = ggplot(aes(x=lon, y=lat), data=ebd_select)) #加上圖層

bubulmap +
  geom_point(aes(x=lon, y=lat, color=Species), data=ebd_select) #加上點

這就是最後完成的圖。粉紅色為白頭翁，藍色為烏頭翁。可以很清楚看到白頭翁分布在西半部（包含澎湖、小琉球等離島），而台灣特有種的烏頭翁分布在台灣東部（含綠島)。兩種的分界大致在南澳蘇花公路，南邊的交界在南迴公路的楓港附近，在南橫公路好像以最高點大關山隧道為界，西邊是白頭翁的地盤，而東側則是烏頭翁的天下。楚河漢界，是不是很有意思呢？

如果要修改圖上的說明、文字、座標等，都可以透過標準的ggplot2語法來做設定。看我們來把圖說移到綠島旁邊。

bubulmap +
  geom_point(aes(x=lon, y=lat, color=Species), data=ebd_select) +
  theme(legend.title=element_blank(),
        legend.position=c(0.85, 0.15),
        legend.background = element_blank(),
        legend.text=element_text(size=15, face="bold", color='white'))+
  scale_color_hue(labels = c("白頭翁", "烏頭翁"))

當然也可以換不同的底圖試試，像這樣顏色又更顯眼了呢。

tw_map_2 <- get_stamenmap(tw_borders, zoom = 8, maptype = "toner")
ggmap(tw_map_2)

bubulmap2<-ggmap(tw_map_2, base_layer = ggplot(aes(x=lon, y=lat), data=ebd_select))
bubulmap2 +
  geom_point(aes(x=lon, y=lat, color=Species), data=ebd_select) +
  theme(legend.title=element_blank(),
        legend.position=c(0.85, 0.15),
        legend.background = element_blank(),
        legend.text=element_text(size=15, face="bold", color='white'))+
  scale_color_hue(labels = c("白頭翁", "烏頭翁"))