以iNatualist觀察紀錄製圖

安裝與加載必要的 R 套件

我們需要使用一些 R 套件來處理和繪圖，像是 tidyverse（處理資料）和 ggplot2（繪圖）。也可能需要 lubridate 來處理時間資料。

# install.packages(c("tidyverse", "lubridate", "httr", "jsonlite"))

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(lubridate)
library(httr)
library(jsonlite)

## 
## Attaching package: 'jsonlite'
## 
## The following object is masked from 'package:purrr':
## 
##     flatten

取得 taxon_id 和 place_id

我們可以使用 iNaturalist API 先查詢物種和地點的 ID。

查詢物種的 taxon_id 例如，查詢八哥屬和藍磯鶇的 taxon_id：

get_taxon_id <- function(species_name, rank) {
  base_taxa_url <- paste0("https://api.inaturalist.org/v1/")
  response <- RETRY(verb = "GET",
                    url = base_taxa_url,
                    path = "v1/taxa",
                    query = list(q = species_name,
                                 rank = rank,
                                 order = "desc")
  )
  data <- fromJSON(content(response, "text"), flatten = TRUE)
  taxon_id <- data$results$id  # 取得第一個匹配的taxon_id
  return(taxon_id)
}

# 取得八哥屬和藍磯鶇的taxon_id
taxon_id_monticola <- get_taxon_id("Monticola solitarius", "species") # 13141
taxon_id_acridotheres <- get_taxon_id("Acridotheres", "genus") # 14868

取得觀察資料

透過 taxon_id 和 place_id 取得觀察資料：

# Define place IDs of interest
places_id <- c(7887, 6737)  # Taiwan = 7887, Japan = 6737
place_id_taiwan <- places_id[1]
place_id_japan <- places_id[2]


# 取得藍磯鶇在台灣和日本的資料
## 藍磯鶇(研究等級)@台灣、日本
## https://www.inaturalist.org/observations?place_id=7887&quality_grade=research&subview=table&taxon_id=13141
## unzip("data/inat_monticola_20240922download/observations-479259.csv.zip", exdir = "data/inat_monticola_20240922download/")
data_monticola <- read_csv("data/inat_monticola_20240922download/observations-479259.csv")

## Rows: 2308 Columns: 39
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (23): observed_on_string, time_observed_at, time_zone, user_login, user...
## dbl  (10): id, user_id, num_identification_agreements, num_identification_di...
## lgl   (5): captive_cultivated, private_place_guess, private_latitude, privat...
## date  (1): observed_on
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## 八哥屬(研究等級)@台灣、日本
## quality_grade=research&identifications=any&place_id=7887,6737&taxon_id=13141&spam=false
## unzip("data/inat_acridotheres_20240922download/observations-479263.csv.zip", exdir = "data/inat_acridotheres_20240922download/")
data_acridotheres <- read_csv("data/inat_acridotheres_20240922download/observations-479263.csv")

## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 6993 Columns: 39
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (23): observed_on_string, time_observed_at, time_zone, user_login, user...
## dbl  (10): id, user_id, num_identification_agreements, num_identification_di...
## lgl   (5): captive_cultivated, private_place_guess, private_latitude, privat...
## date  (1): observed_on
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

合併與繪圖

利用 ggplot2 來繪製物種隨時間的分布圖。

library(rnaturalearth)

## Warning: package 'rnaturalearth' was built under R version 4.3.2

# extract eastern Asia region from world
world <- ne_countries(scale = "medium", continent = "asia")
# targetRegion2 <- world %>% filter(region_wb=="South Asia" | region_wb == "East Asia & Pacific" & subregion != "Australia and New Zealand")

targetRegion <- world %>% filter(admin %in% c("Japan", "Taiwan"))

# Summarize and prepare data for plotting
data_monticola_clean <- data_monticola %>%
  drop_na(time_observed_at) %>% 
  mutate(year = year(time_observed_at)) %>% 
  mutate(birdTaxa = "藍磯鶇")
data_acridotheres_clean <- data_acridotheres %>%
  drop_na(time_observed_at) %>% 
  mutate(year = year(time_observed_at)) %>% 
  mutate(birdTaxa = "八哥屬")
data_combine <- bind_rows(data_monticola_clean, data_acridotheres_clean) %>% 
  select(birdTaxa, year, longitude, latitude) %>% 
  mutate(yearInterval = cut_interval(year, length = 3))


# draw map by year
p <- ggplot() +
  geom_sf(data = targetRegion) +
  geom_point(data = data_combine, mapping = aes(x = longitude, y = latitude, group = birdTaxa, color = birdTaxa), shape = 15) +
  # geom_point(data = data_combine %>% filter(birdTaxa == "八哥屬"), mapping = aes(x = longitude, y = latitude), color = "#ffff0060") +
  facet_wrap(~yearInterval, nrow = 2) +
  scale_fill_brewer() +
  scale_color_manual(values = c("#ffff0040", "#0000ff60")) +
  theme_classic() +
  theme(axis.text.x = element_blank(),axis.text.y=element_blank(),
        axis.line = element_blank(),axis.ticks = element_blank())
p

台灣島及周邊島嶼地區繪圖

data_combine_taiwan <- bind_rows(data_monticola_clean, data_acridotheres_clean) %>% 
  select(birdTaxa, year, longitude, latitude) %>% 
  mutate(yearInterval = cut_interval(year, length = 2)) %>% 
  filter(longitude <= 122.5)

# draw map by year in Taiwan
targetRegion_taiwan <- world %>% filter(admin %in% c("Taiwan"))


p_taiwan <- ggplot() +
  geom_sf(data = targetRegion_taiwan) +
  geom_point(data = data_combine_taiwan, mapping = aes(x = longitude, y = latitude, group = birdTaxa, color = birdTaxa), shape = 15) +
  # geom_point(data = data_combine %>% filter(birdTaxa == "八哥屬"), mapping = aes(x = longitude, y = latitude), color = "#ffff0060") +
  facet_wrap(~yearInterval, nrow = 3) +
  scale_fill_brewer() +
  scale_color_manual(values = c("#ffff0020", "#0000ff60")) +
  theme_classic() +
  theme(axis.text.x = element_blank(),axis.text.y=element_blank(),
        axis.line = element_blank(),axis.ticks = element_blank())
p_taiwan

換成TBN資料

取得觀察資料 from TBN

於TBN網站以關鍵字”藍磯鶇”、“八哥屬”查找資料

dataTBN_monticola <- read_csv("data/tbn_monticola_20240922download/data.csv")

## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 43571 Columns: 43
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (30): 紀錄識別號, 原始紀錄物種, 觀測方式, 外部識別碼, 行政區, 記錄者, 數量單位, 類群, 物種UUID, TaiCOL編碼,...
## dbl   (9): 觀測年, 觀測月, 觀測日, 原始座標誤差(公尺), 緯度(十進位), 經度(十進位), 海拔, 個體數量, 數量...
## lgl   (2): 保育類等級, 敏感資料模糊化程度
## dttm  (1): 最後修改時間
## date  (1): 觀測日期
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

dataTBN_acridotheres <- read_csv("data/tbn_acridotheres_202420922download/data.csv")

## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 555806 Columns: 43
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (31): 紀錄識別號, 原始紀錄物種, 觀測方式, 外部識別碼, 行政區, 記錄者, 數量單位, 類群, 物種UUID, TaiCOL編碼,...
## dbl  (10): 觀測年, 觀測月, 觀測日, 原始座標誤差(公尺), 緯度(十進位), 經度(十進位), 海拔, 個體數量, 數量, 敏感資料模糊...
## dttm  (1): 最後修改時間
## date  (1): 觀測日期
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

合併與繪圖 from TBN