Calculate departure from average ocean temperature for 1980-2019

PottmpOcean.R

## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
# Needed so the code runs later on.
datadir <- "pottmp_data"
fnames <- paste("pottmp",1980:2019,"nc", sep = ".")

## ----Download information, eval=FALSE------------------------------------
## # ONLY RUN THIS ONCE TO DOWNLOAD THE DATA
## # This downloads 39 files of appx 144 Mb each!
## # Create a directory to store the data.
## datadir <- "pottmp_data"
## if(!dir.exists(datadir)) dir.create(datadir)
## # Create all the filenames for downloading
## fnames <- paste("pottmp",1980:2019,"nc", sep = ".")
## 
## ftpdir <- "ftp://ftp.cdc.noaa.gov/Datasets/godas/"
## 
## for(i in fnames){
##   download.file(
##     url = paste0(ftpdir, i),
##     destfile = file.path(datadir, i),
##     mode = "wb" #download in binary mode or files can't be read!
##   )
##   Sys.sleep(30) # wait 30 seconds before downloading the next file.
## }
## 

## ----Get info from netCDF files------------------------------------------
library(ncdf4)
ncin <- nc_open(file.path(datadir, "pottmp.2019.nc"))

## ------------------------------------------------------------------------
pottmp2019 <- ncvar_get(ncin, "pottmp")
lon <- ncvar_get(ncin, "lon")
lat <- ncvar_get(ncin, "lat")
# Don't leave your file open!
nc_close(ncin) 
# Check dimensions
dim(pottmp2019)

## ------------------------------------------------------------------------
range(pottmp2019[,,1,2], na.rm = TRUE) 

## ------------------------------------------------------------------------
# calculate the mean of the 10 upper layers (5m - 105m)
tmpdata <- apply(pottmp2019[,,1:10,],
                 c(1,2,4),  # this way you take the mean of the 10 upper layers
                            # for every longitude(1), latitude(2) and month(4)
                 mean, na.rm = TRUE)
tmpdata <- tmpdata - 273.15

## ---- warning=FALSE------------------------------------------------------
# Quickly plot the data
thebreaks <- seq(-4,31, by = 1)
thecols <- rev(rainbow(length(thebreaks) - 1))

library(maps)
image(lon,lat, tmpdata[,,2], breaks = thebreaks, col = thecols)
map("world2", add = TRUE, fill = TRUE, col = "grey")

## ----Close the file again------------------------------------------------
pottmp1980_2019 <- lapply(fnames, function(x){
  ncin <- nc_open(file.path(datadir,x))
  on.exit(nc_close(ncin)) # prevents files staying open when an error occurs
  # obtain the data we want
  tmp <- ncvar_get(ncin, "pottmp",         # the variable of interest
                   count = c(-1,-1,10,-1)  # only read 10 upper layers of ocean
                   )
  # calculate grand mean by month and return
  apply(tmp, 4, mean, na.rm = TRUE) 
  
})
pottmp1980_2019 <- unlist(pottmp1980_2019) - 273.15

## ------------------------------------------------------------------------
# Calculate a numeric time variable for plotting
xtime <- rep(1980:2019, each = 12) + (0.5:11.5)/12
xtime <- head(xtime, -10)
# Put everything together in a dataset
pottmp <- data.frame(
  tmp = pottmp1980_2019,
  year = floor(xtime),
  month = factor(c(rep(month.abb,times = 39), "Jan","Feb"),
                 levels = month.abb)
)

## ------------------------------------------------------------------------
refmeans <- with(pottmp[pottmp$year %in% 1981:2010,], # Select years
                 {
                   tapply(tmp, month, mean, na.rm = TRUE) # mean by month
                 })
pottmp$anomaly <- pottmp$tmp - head(rep(refmeans,40),-10) # calculate difference

## ------------------------------------------------------------------------
plot(anomaly ~ xtime, data = pottmp, type = "n", 
     xlab = "", ylab = "Departure from 1981-2010 average",
     main = "Ocean temp in upper 105 meter (GODAS)")
abline(h = seq(-0.1,0.4, by = 0.1), col = "lightgrey")
lines(anomaly ~ xtime, data = pottmp, lwd = 2, col = "red")
points(anomaly~xtime, data = pottmp, pch = 19, cex = 0.5)

## ------------------------------------------------------------------------
library(ggplot2)

pottmp$date <- as.Date(paste(pottmp$year, pottmp$month, 15, sep = "-"), 
                       format = "%Y-%b-%d")

ggplot(pottmp, aes(x = date, y = anomaly)) +
  geom_line(lwd = 2, col = "red") + geom_point() +
  labs(x = "Date", y = "Departure from 1981-2010 average") +
  scale_x_date(date_breaks = "2 years", date_minor_breaks = "1 year",
               date_labels = "%b %Y") +
  theme(axis.text.x.bottom = element_text(angle = 90, vjust = 0.5)) +
  ggtitle("Ocean temp in upper 105 meter (GODAS)")