lme_extractions.Rmd

---
title: "OBIS extractions from Large Marine Ecosystem regions"
output:
  html_document:
    df_print: paged
---
# Marine Biodiversity Observation Network Pole to Pole of the Americas ([MBON Pole to Pole](https://marinebon.org/p2p/))

Written by E. Montes (enrique.montes@noaa.gov) on Auguts 28, 2020.


This code creates a map showing the boundaries of a selected Large Marine Ecosystems ([LME](http://lme.edc.uri.edu/index.php/lme-introduction)) and extracts records for selected taxa from the [Ocean Biodiversity Information System (OBIS)](https://obis.org/) using [robis](https://obis.org/manual/accessr/) tools. 

# Step 1
First, let's load required libraries
```{r}
library(robis)
library(rgdal) # for `ogrInfo()` and `readOGR()`
library(tools) # for `file_path_sans_ext()`
library(dplyr) # for `inner_join()`, `filter()`, `summarise()`, and the pipe operator (%>%)
library(ggplot2) # for `fortify()` and for plotting
library(sp) # for `point.in.polygon()` and `spDists()`
library(tidyr) # for `gather()`
library(readr) # for `write_tsv()`
library(leaflet)
library(lubridate)
```

# Step 2 - NOT NEEDED ANY MORE
Now let's provide the function fortify.shape(), which puts the shapefile data in the object class data.frame, so that it can be used by ggplot2, and extract portions of the data (from the fortified data.frame object) for a smaller domain
# ```{r}
# fortify.shape <- function(x){
#   x@data$id <- rownames(x@data)
#   x.f <- fortify(x, region = "id")
#   x.join <- inner_join(x.f, x@data, by = "id")
# }
# 
# subset.shape <- function(x, domain){
#   x.subset <- filter(x, long > domain[1] & 
#                        long < domain[2] & 
#                        lat > domain[3] & 
#                        lat < domain[4])
#   x.subset
# }
# ```

# Step 3
Let's read the shapefile "LMEs66.shp" containing all polygons, fortify the global data and then extract domain. See numbers [here](http://lme.edc.uri.edu/index.php/lme-introduction)
```{r}

setwd("/Users/enrique.montes/Documents/lme-extractions/data")
path.lme.coast <- ("~/Documents/lme-extractions/data")
fnam.lme.coast <- "LMEs66.shp"
dat.coast <- readOGR(dsn = path.lme.coast,
                     layer = file_path_sans_ext(fnam.lme.coast))

# fortify the global data and then extract domain
# dat.coast <- fortify.shape(dat.coast) # a 410951x8 dataframe # NOT NEEDED

# Specify the desired LME:
dat.sel_1 <- subset(dat.coast, LME_NUMBER == 15) # S. Brazil. 
dat.sel_2 <- subset(dat.coast, LME_NUMBER == 16) # E. Brazil
dat.sel_3 <- subset(dat.coast, LME_NUMBER == 17) # N. Brazil
dat.sel_4 <- subset(dat.coast, LME_NUMBER == 14) # Patagonia
dat.sel_5 <- subset(dat.coast, LME_NUMBER == 13) # Humboldt C.
dat.sel_6 <- subset(dat.coast, LME_NUMBER == 12) # Caribbean
dat.sel_7 <- subset(dat.coast, LME_NUMBER == 11) # P. Ctral A.
dat.sel_8 <- subset(dat.coast, LME_NUMBER == 5) # GoM
dat.sel_9 <- subset(dat.coast, LME_NUMBER == 4) # Gulf of California
dat.sel_10 <- subset(dat.coast, LME_NUMBER == 3) # CCS
dat.sel_11 <- subset(dat.coast, LME_NUMBER == 2) # G. Alaska
dat.sel_12 <- subset(dat.coast, LME_NUMBER == 1) # E. Bearing S.
dat.sel_13 <- subset(dat.coast, LME_NUMBER == 54) # Chukchi S.
dat.sel_14 <- subset(dat.coast, LME_NUMBER == 55) # Beauford S.
dat.sel_15 <- subset(dat.coast, LME_NUMBER == 66) # Canadian Arctic
dat.sel_16 <- subset(dat.coast, LME_NUMBER == 18) # Canadian E. Arctic
dat.sel_17 <- subset(dat.coast, LME_NUMBER == 9) # Labrador S.
dat.sel_18 <- subset(dat.coast, LME_NUMBER == 8) # Scotian Shelf
dat.sel_19 <- subset(dat.coast, LME_NUMBER == 7) # NE US
dat.sel_20 <- subset(dat.coast, LME_NUMBER == 6) # SE US
dat.sel_21 <- subset(dat.coast, LME_NUMBER == 10) # Hawaii
dat.sel_22 <- subset(dat.coast, LME_NUMBER == 63) # Hudson Bay Complex
dat.sel_23 <- subset(dat.coast, LME_NUMBER == 19) # Greenland Sea
dat.sel_24 <- subset(dat.coast, LME_NUMBER == 21) # Norwegian Sea
dat.sel_25 <- subset(dat.coast, LME_NUMBER == 20) # Barents Sea
dat.sel_26 <- subset(dat.coast, LME_NUMBER == 58) # Kara Sea
dat.sel_27 <- subset(dat.coast, LME_NUMBER == 57) # Laptev Sea
dat.sel_28 <- subset(dat.coast, LME_NUMBER == 56) # E. Sibarian Sea
```

# Step 4
Plotting the coastline and selected LME boundaries
```{r}
# Define lat/lon limits here
xlims <- c(-150, -25)
ylims <- c(-60, 60)

# Generate a base map with the coastline:
p0 <- ggplot() + theme(text = element_text(size=18)) + 
  geom_path(data = dat.coast, aes(x = long, y = lat, group = group), 
            color = "black", size = 0.25) + 
  coord_map(projection = "mercator") + 
  scale_x_continuous(limits = xlims, expand = c(0, 0)) + 
  scale_y_continuous(limits = ylims, expand = c(0, 0)) + 
  labs(list(title = "", x = "Longitude", y = "Latitude"))

p0

# highlight LME of interest
p.sel <- p0 +
  geom_path(data = dat.sel_6,
            aes(x = long, y = lat, group = group),
            colour = "chartreuse4", size = 1) +
  geom_path(data = dat.sel_8,
            aes(x = long, y = lat, group = group),
            colour = "coral", size = 0.75)
p.sel

```

# Step 5
This section extracts OBIS records of using the "occurrence" function or downloaded data, and plots time series or pie charts with distributions of large groups (annelids, molluscs, plants and echinoderms) in the upper 100m.
```{r}
# Set extraction params
# LME codes (from OBIS URL)
# N Brazil=40017; E Brazil=40016; S Brazil=40015; Patagonia=40014; Humboldt=40013; Caribbean=40012; 
# P Ctral A=40011; CCS=40003; GoA=40002; NE USA=40007; E Bearing=40001; Canada E Arctic=40018; GoM=40005; 
# Chukchi=40054; SE USA=40006; Labrador=40009; Scotian S=40008

area = 40005  
depth = 100
mol_code = 51
echi_code = 1806
anne_code = 882
plan_code = 3

## read data from OBIS
## Molluscs
mollusc.100_2 = occurrence(areaid = area, taxonid = mol_code, enddepth = depth) 

## Ehinoderms
echino.100_2 = occurrence(areaid = area, taxonid = echi_code, enddepth = depth)

## Annelida
anne.100_2 = occurrence(areaid = area, taxonid = anne_code, enddepth = depth)

## Platae
plant.100_2 = occurrence(areaid = area, taxonid = plan_code, enddepth = depth)

## Merge the data frames
#subset data (year and phylum)
sub_mollusc <- data.frame(mollusc.100_2$date_year, mollusc.100_2$phylum) 
# rename column headers
names(sub_mollusc)[names(sub_mollusc) == "mollusc.100_2.date_year"] <- "year"
names(sub_mollusc)[names(sub_mollusc) == "mollusc.100_2.phylum"] <- "phylum"

sub_echino <- data.frame(echino.100_2$date_year, echino.100_2$phylum)
names(sub_echino)[names(sub_echino) == "echino.100_2.date_year"] <- "year"
names(sub_echino)[names(sub_echino) == "echino.100_2.phylum"] <- "phylum"

sub_anne <- data.frame(anne.100_2$date_year, anne.100_2$phylum)
names(sub_anne)[names(sub_anne) == "anne.100_2.date_year"] <- "year"
names(sub_anne)[names(sub_anne) == "anne.100_2.phylum"] <- "phylum"

sub_plant <- data.frame(plant.100_2$date_year, plant.100_2$phylum)
names(sub_plant)[names(sub_plant) == "plant.100_2.date_year"] <- "year"
names(sub_plant)[names(sub_plant) == "plant.100_2.phylum"] <- "phylum"

total.100 <- bind_rows(sub_mollusc, sub_echino, sub_anne, sub_plant)

## Plot the data
ts_plot <- ggplot() +
  geom_histogram(data = total.100, aes(x = year, fill = phylum), binwidth = 2) + 
  scale_fill_brewer(palette = "Spectral") + 
  xlim(c(1960, 2025)) + 
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  theme(axis.text.x = element_text(size=14, angle=0), 
        axis.text.y = element_text(size=14, angle=0))
ts_plot

# ggsave(ts_plot, filename = "test.png", device = "png", width = 20, height = 10,  dpi=300)

## Plot pie chart 
# Group plants in a single group
all_tbl <- table(total.100$phylum)
all_df <- as.data.frame(all_tbl)
p_list = c("Chlorophyta", "Rhodophyta", "Tracheophyta")
rest_list = c("Annelida", "Echinodermata", "Mollusca")
p_idx = match(p_list, rownames(all_tbl))
rest_idx = match(rest_list, rownames(all_tbl))
p_sum = sum(all_df$Freq[p_idx], na.rm = TRUE)
freq_val <- c(all_df$Freq[rest_idx], p_sum)
group_id <- c(c(rest_list), "Plants")
f_tbl <- data.frame(group = group_id, freq = freq_val)
  
cols <- rainbow(nrow(f_tbl))
groups_pie <- pie(f_tbl$freq, labels = f_tbl$group, col = cols)
  
```


# Extract all records from an area and visualize results
```{r}
obis_df <- occurrence(areaid = area)
obis_df[complete.cases(obis_df[ , "date_year"]), ]

allts_plot <- ggplot() +
  geom_histogram(data = obis_df, aes(x = date_year, fill = kingdom), binwidth = 1) + 
  scale_fill_brewer(palette = "Spectral") + 
  # guides(fill="none") +
  xlim(c(1940, 2025)) +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  theme(axis.text.x = element_text(size=14, angle=0), 
        axis.text.y = element_text(size=14, angle=0))
allts_plot

last_df <- obis_df %>% 
  filter(date_year == 2021)
unique(last_df$datasetID)
```