-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy patheez_extractions.Rmd
147 lines (120 loc) · 4.84 KB
/
eez_extractions.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
---
title: "OBIS extractions from Exclusive Economic Zones"
output:
html_document:
df_print: paged
pdf_document: default
---
# Marine Biodiversity Observation Network Pole to Pole of the Americas ([MBON Pole to Pole](https://marinebon.org/p2p/))
Written by E. Montes ([email protected]) on June 23, 2023.
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, warning=FALSE}
library(ggplot2)
library(rgdal) # for `ogrInfo()` and `readOGR()`
library(sf)
library(tools)
library(robis)
library(tidyverse)
```
# Step 2
Load shapefile
```{r, warning=FALSE}
# Download the EEZ shapefile for the Americas
setwd("/Users/enrique.montes/Documents/lme-extractions/EEZ_land_union_v3_202003")
path.eez.coast <- ("~/Documents/lme-extractions/EEZ_land_union_v3_202003")
fnam.eez.coast <- "EEZ_Land_v3_202030.shp"
eez <- readOGR(dsn = path.eez.coast,
layer = file_path_sans_ext(fnam.eez.coast))
# Read the shapefile
sel.eez<- subset(eez, SOVEREIGN1 == "Ecuador")
```
# Step 3
Plotting the coastline and selected EEZ boundaries
```{r, warning=FALSE}
# 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 = eez, 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 = sel.eez,
aes(x = long, y = lat, group = group),
colour = "coral", size = 0.75)
p.sel
```
# Step 4
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, warning=FALSE}
# Set extraction params
# EEZ codes (from OBIS URL)
# Ecuador = 58 (https://obis.org/area/58)
area = 58
depth = 100
mol_code = 51
echi_code = 1806
anne_code = 882
plan_code = 3
## read data from OBIS
## Molluscs
molluscs= occurrence(areaid = area, taxonid = mol_code, enddepth = depth)
print(molluscs$species)
## Ehinoderms
echinos = occurrence(areaid = area, taxonid = echi_code, enddepth = depth)
print(echinos$species)
## Annelida
annel = occurrence(areaid = area, taxonid = anne_code, enddepth = depth)
print(annel$species)
## Platae
plants = occurrence(areaid = area, taxonid = plan_code, enddepth = depth)
print(plants$species)
## Merge the data frames
sub_mollusc <- data.frame(molluscs$date_year, molluscs$phylum)
# rename column headers
names(sub_mollusc)[names(sub_mollusc) == "molluscs.date_year"] <- "year"
names(sub_mollusc)[names(sub_mollusc) == "molluscs.phylum"] <- "phylum"
sub_echino <- data.frame(echinos$date_year, echinos$phylum)
names(sub_echino)[names(sub_echino) == "echinos.date_year"] <- "year"
names(sub_echino)[names(sub_echino) == "echinos.phylum"] <- "phylum"
sub_annel <- data.frame(annel$date_year, annel$phylum)
names(sub_annel)[names(sub_annel) == "annel.date_year"] <- "year"
names(sub_annel)[names(sub_annel) == "annel.phylum"] <- "phylum"
sub_plant <- data.frame(plants$date_year, plants$phylum)
names(sub_plant)[names(sub_plant) == "plants.date_year"] <- "year"
names(sub_plant)[names(sub_plant) == "plants.phylum"] <- "phylum"
total <- bind_rows(sub_mollusc, sub_echino, sub_annel, sub_plant)
## Plot the data
ts_plot <- ggplot() +
geom_histogram(data = total, aes(x = year, fill = phylum), binwidth = 2) +
scale_fill_brewer(palette = "Spectral") +
xlim(c(1960, 2017)) +
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$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)
```