-
Notifications
You must be signed in to change notification settings - Fork 0
/
03-Oceanic_pH.Rmd
161 lines (132 loc) · 5.39 KB
/
03-Oceanic_pH.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
148
149
150
151
152
153
154
155
156
157
158
159
160
---
title: "Oceanic pH"
author: "[email protected] and [email protected]"
date: '2024-6-6'
output: md_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## __Indicator: Oceanic pH__
Oceanic pH is a measure of how greenhouse gas emissions have already impacted the ocean. This indicator demonstrates that oceanic pH has decreased significantly over the past several decades (i.e., the ocean has become more acidic). Increasing ocean acidification limits the ability of marine organisms to build shells and other calcareous structures. Recent research has shown that pelagic organisms such as pteropods and other prey for commercially valuable fish species are already being negatively impacted by increasing acidification (Feely et al. 2016). The full impact of ocean acidification on the pelagic food web is an area of active research (Fabry et al. 2008).
```{r}
### Load libraries
library(tidyverse)
library(lubridate)
library(here)
library(stringr)
library(nmfspalette)
library(ncdf4)
```
```{r}
### Load libraries for mapping
library(raster)
library(rasterVis)
library(mapdata)
library(maptools)
library(cmocean)
library(latticeExtra)
library(grid)
library(rerddap)
library(terra)
library(viridis)
```
```{r}
# Set report year (RptYr), to make things easier
RptYr <- 2020
# Set path to variable: Sea_Surface_Temperature
# This is where the data are and where the plots will go
Dir <- here("Oceanic_pH")
```
```{r}
### Load data
# Bounding box, from John Marra via email:
lon_range <- c(129.4088, 137.0541)
lat_range <- c(1.5214, 11.6587)
#pre-downloaded netCDF file for pH at 5.14 m within the Palau domain (1,5-11.75N,129.5-137W)
fn="ph_palau.nc"
ph_data <- nc_open(fn)
ph <- ncvar_get(ph_data, "ph")
time <- ncvar_get(ph_data, "time")
nc_close(ph_data)
```
```{r}
# Change time from hours to year-month-date format
time_df <- as.data.frame.table(as.POSIXct(time*3600 ,origin='1950-01-01 00:00'))[2]#['Freq']
colnames(time_df) <- c('time')
# Monthly spatial average
ph_df=as.data.frame(apply(ph, c(3), mean, na.rm = TRUE))
colnames(ph_df) <- c('pH')
ph_ts <- data.frame(time_df, ph_df)
```
```{r}
### Linear fit
n_obs <- seq(1, length(ph_ts$pH), 1)
ph_lm <- lm(ph_ts$pH~ n_obs)
# summary(ph_lm) shows that there's a significant decreasing
# trend over time for this indicator. There are some automated
# checks that can be added to make sure this is still the case
# in future years.
# Change over time
delta_ph <- ph_lm$fitted.values[length(n_obs)] - ph_lm$fitted.value[1]
```
```{r}
### Plot the time series
# Create axes limits to make things simpler
# These were determined through looking at quick rough plots and data limits
ph_xlim <- c(min(ymd_hms(ph_ts$time)), max(ymd_hms(ph_ts$time)))
ph_ylim <- c(8.01, 8.09)
# Access the NMFS color palette
oceans <- nmfs_palette("oceans")(3)
crustacean <- nmfs_palette("crustacean")(4)
# Plot
plot(ymd_hms(ph_ts$time), ph_ts$pH, type = "l", lwd = 2, col = oceans[2],
xlim = ph_xlim, ylim = ph_ylim, xlab = " ", ylab = "pH",
xaxt = "n", yaxt = "n", xaxs = "i", yaxs = "i")
par(new = TRUE)
plot(ymd_hms(ph_ts$time), ph_lm$fitted.values, type = "l", lwd = 2, col = crustacean[1],
xlim = ph_xlim, ylim = ph_ylim, xlab = " ", ylab = " ",
xaxt = "n", yaxt = "n", xaxs = "i", yaxs = "i")
axis((1), at = ymd_hms(ph_ts$time[seq(1, length(n_obs), 12)]), tck = 0.025, labels = year(make_date(seq(1993, RptYr, 1))))
axis((2), at = seq(8.01, 8.09, 0.01), tck = 0.025, las = 1)
axis((3), at = ymd_hms(ph_ts$time[seq(1, length(n_obs), 12)]), tck = 0.025, labels = FALSE)
axis((4), at = seq(8.01, 8.09, 0.01), tck = 0.025, labels = FALSE)
# _axt = "n" removes tick labels so that they can be customized later
# _axs = "i" removes whitespace beyond axes maxima
```
```{r}
# Read the file into R and make it to rasterstack
stack_ph = stack(fn)
# Convert raster data to dataframe for calculating climatology
df_temp = as.data.frame(rasterToPoints(stack_ph))
df_temp$z = rowMeans(df_temp[,3:dim(df_temp)[2]], na.rm = T)
# Convert dataframe to raster for mapping
rst = rasterFromXYZ(df_temp[,c("x", "y", "z")])
#create a rasterbrick with the rasterlayer
ph_clim <- brick(rst)
```
```{r}
### Mapping long-term climatology
# Get land information and make it into a spatial object
land <- maps::map('world', fill=TRUE, xlim=lon_range, ylim=lat_range, plot=FALSE)
ids <- sapply(strsplit(land$names, ":"), function(x) x[1])
bPols <- map2SpatialPolygons(land, IDs=ids, proj4string=CRS('+proj=longlat +datum=WGS84 +no_defs'))
# Add EEZ
llines.SpatVector <- function(x, ...) {
xy <- crds(x, list=TRUE)
names(xy) <- c("x", "y")
lattice::llines(xy, ...)
}
f <- "eez/eez.shp"
v <- vect(f)
lns <- as.lines(v)
# make map themes
mapTheme <- rasterTheme(region=cmocean('matter')(50))
# Make plot
#regular scale
levelplot(ph_clim , pretty=T, margin=F, par.setting=mapTheme, colorkey=list( height = .5, width = 1) ) + layer(sp.polygons(bPols)) + layer(llines(lns))
# add unit to colorbar
grid.text(expression(pH) , y=unit(0.6, "npc"),
x=unit(0.81, "npc"))
```
Oceanic pH is simulated by biogeochemical models using ocean and atmosphere reanalyses as forcings. It spans from 1993 to 2020 and is provided by Copernicus (https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_BGC_001_029/description). A highly significant (p<0.001) decreasing trend is found in the pH of Palau waters (5m) from 1993 to `r RptYr`, and pH declined by `r signif(abs(delta_ph),2)` over this period.