Trong bài này chúng ta sẽ cùng thực hành cách lấy dữ liệu từ GBIF để vẽ vùng phân bố cho chi Hoya tại Việt Nam

1 Load dữ liệu

Có nhiều cách để có thể lấy được data trực tiếp từ GBIF với các gói thư viện như spocc, rgbif, dismo. Ở đây chúng ta sẽ thực hiện với gói spocc.

Ta đặt điều kiện has_coords = T để chỉ tải các ghi nhận có thông tin về tọa độ

library(spocc)
hoya <- occ(query = "Hoya", from= 'gbif', has_coords = T)

Kết quả trả về là một dạng list, để lọc các ghi nhận của Việt Nam ta thực hiện như sau

hoya_vn <- subset(hoya$gbif$data$Hoya, country == "Viet Nam")
dim(hoya_vn)
## [1]  56 112

Kết quả cho thấy có 56 ghi nhận về Hoya tại Việt Nam ( chỉ tính các ghi nhận có tọa độ) với 112 trường thông tin. có thể xem với lệnh View(hoya_vn) hoặc str(hoya_vn)

Tới đây ta đã có được dữ liệu cần, tuy nhiên trước khi vẽ bản đồ ta cần qua bước làm sạch dữ liệu

2 Làm sạch dữ liệu

Để làm sạch dữ liệu ta có thể sử dụng nhiều gói thư viện khác nhau trong R. Tuy nhiên đối với dữ liệu về ghi nhận đa dạng sinh học được tải về từ GBIF như vầy thì chúng ta có gói thư viện dành riêng là scrubr

library(scrubr)
library(magrittr) # để sử dụng được %>%

2.1 Tọa độ

Làm sạch các dữ liệu về tọa độ với coord_impossible() để loại ghi nhận có tọa độ không thích hợp, coord_incomplete() cho tọa độ không hoàn chỉnh,coord_unlikely() cho trường hợp tọa độ không chắc

dframe(hoya_vn) %>%
  coord_impossible() %>%
  coord_incomplete() %>%
  coord_unlikely()
## # A tibble: 56 × 112
##            name longitude latitude issues  prov        key
## *         <chr>     <dbl>    <dbl>  <chr> <chr>      <int>
## 1          <NA>    103.98    22.05 gass84  gbif 1096847394
## 2          <NA>    103.98    22.03 gass84  gbif 1096847408
## 3          <NA>    105.72    20.25 gass84  gbif 1096847405
## 4          <NA>    105.67    20.29 gass84  gbif 1096847407
## 5  Hoya carnosa    104.74    19.06 gass84  gbif 1096847427
## 6          <NA>    103.98    22.05 gass84  gbif 1096847402
## 7          <NA>    105.72    20.25 gass84  gbif 1096847422
## 8          <NA>    105.72    20.25 gass84  gbif 1096847406
## 9          <NA>    105.64    21.47 gass84  gbif 1096847419
## 10         <NA>    105.65    21.45 gass84  gbif 1096847431
## # ... with 46 more rows, and 106 more variables: datasetKey <chr>,
## #   publishingOrgKey <chr>, publishingCountry <chr>, protocol <chr>,
## #   lastCrawled <chr>, lastParsed <chr>, crawlId <int>,
## #   basisOfRecord <chr>, individualCount <int>, taxonKey <int>,
## #   kingdomKey <int>, phylumKey <int>, classKey <int>, orderKey <int>,
## #   familyKey <int>, genusKey <int>, scientificName <chr>, kingdom <chr>,
## #   phylum <chr>, order <chr>, family <chr>, genus <chr>,
## #   genericName <chr>, specificEpithet <chr>, infraspecificEpithet <chr>,
## #   taxonRank <chr>, coordinatePrecision <dbl>,
## #   coordinateUncertaintyInMeters <dbl>, elevation <dbl>,
## #   elevationAccuracy <dbl>, stateProvince <chr>, year <int>, month <int>,
## #   day <int>, eventDate <date>, lastInterpreted <chr>, license <chr>,
## #   geodeticDatum <chr>, class <chr>, countryCode <chr>, country <chr>,
## #   eventID <chr>, identifier <chr>, catalogNumber <chr>,
## #   recordedBy <chr>, vernacularName <chr>, institutionCode <chr>,
## #   locality <chr>, county <chr>, occurrenceRemarks <chr>,
## #   collectionCode <chr>, gbifID <chr>, occurrenceID <chr>,
## #   continent <chr>, recordNumber <chr>, higherGeography <chr>,
## #   endDayOfYear <chr>, type <chr>, rights <chr>, startDayOfYear <chr>,
## #   higherClassification <chr>, verbatimEventDate <chr>,
## #   nomenclaturalCode <chr>, verbatimCoordinateSystem <chr>,
## #   datasetName <chr>, disposition <chr>, fieldNotes <chr>,
## #   ownerInstitutionCode <chr>, samplingProtocol <chr>,
## #   verbatimElevation <chr>, collectionID <chr>, eventTime <chr>,
## #   dateIdentified <chr>, modified <chr>, references <chr>,
## #   rightsHolder <chr>, verbatimLocality <chr>, taxonID <chr>,
## #   `http://unknown.org/occurrenceDetails` <chr>, identificationID <chr>,
## #   habitat <chr>, georeferenceProtocol <chr>, language <chr>,
## #   otherCatalogNumbers <chr>, institutionID <chr>, identifiedBy <chr>,
## #   fieldNumber <chr>, preparations <chr>, previousIdentifications <chr>,
## #   nomenclaturalStatus <chr>, bibliographicCitation <chr>,
## #   waterBody <chr>, island <chr>, accessRights <chr>,
## #   verbatimTaxonRank <chr>, municipality <chr>,
## #   identificationQualifier <chr>, created <chr>,
## #   associatedSequences <chr>, locationRemarks <chr>, ...

Ngoài ra, gói này còn có thể làm sạch dữ liệu về thời gian ghi nhận.

2.2 Tên loài NA

Trong dữ liệu có nhiều trường hợp thiếu thông tin về tên loài của ghi nhận, ta tạm thời thay các trường hợp này bởi “Hoya sp.”

hoya_vn$name[which(is.na(hoya_vn$name)==TRUE)]<-"Hoya sp."

2.3 Sau khi làm sạch, dữ liệu ta có thể vẽ các biểu đồ theo nhu cầu của mình

Ví dụ: Ta có thể xem nhanh số lượng ghi nhận theo năm như sau:

hist(hoya_vn$year, col="cyan",
     xlab="Nam", ylab="So luong ghi nhan",
     main="So luong ghi nhan ve Hoya theo nam",
     xlim=c(1995,2015),ylim=c(0,12))

3 Bản đồ phân bố

Vẽ bản đồ ta cũng có nhiều cách, có thể thực hiện ví dụ một số trường hợp như sau:

3.1 Bản đồ động với Leaflet

Với dạng bản đồ này ta có thể phóng to, thu nhỏ, hay nhấn vào điểm để để xem thông tin

library(mapr)
map_leaflet(hoya_vn)

3.2 Bản đồ trên nền google map

library(ggplot2)
map_ggmap(hoya_vn, zoom=5, size=3) + ggtitle("Ban do ghi nhan Hoya tai Viet Nam (Data: GBIF)")

3.3 Vẽ nhanh trên nền phần đất liền

library(maps)
map('world','Vietnam', col="red")
points(hoya_vn$longitude, hoya_vn$latitude, col="blue", pch=20, cex=0.8)
map.axes()

3.4 Hay vẽ chi tiết hơn với ggplot2

library(ggplot2)
library(maptools)
vn_base <- readShapePoly("map/vietnam_base.shp")
#can chuyen sang data.frame
vn_base <- fortify(vn_base)

#plot
g <- ggplot()+ 
    geom_polygon(data=vn_base, aes(long, lat, group = group, fill = hole), color = alpha("darkred", 1/2), size = 0.7) + 
    scale_fill_manual(values = c("skyblue", "white"), guide=F) 
g + geom_point(data=hoya_vn, aes(x=longitude,y=latitude, shape=name, color=name),size=4)+
    scale_shape_manual(values=c(1:8), name="Species")+
    scale_color_manual(values=c(1:7,11), name="Species")

3.5 Bản đồ gộp số lượng ghi nhận trên leaflet

library(leaflet)
leaflet(data = hoya_vn) %>% addTiles() %>%
    addMarkers(~longitude, ~latitude,clusterOptions = markerClusterOptions())