Time series analysis of data for Sea Surface Temperature(SST)

1.Introduction

The Bay of Bengal (BoB) is the north-eastern part of the Indian Ocean, bounded on the west and northwest by India, on the north by Bangladesh, and on the east by Myanmar and the Andaman and Nicobar Islands of India. Its southern limit is a line between Sangaman Kanda, Sri Lanka and the north westernmost point of Sumatra(Wikipedia).
The Sea Surface Temperature (SST) of the ocean is indicated by measurements taken at depths that range from 1 millimeter to 20 meters.Some measurements are made using shipboard instruments, but satellites now provide the majority of global SST data.Global climate change is largely dependent on the long term changes in the Sea Surface Temperature Anomaly (SSTA). Regional climate in the Bangladesh, Myanmar and a large part of Indian sub-continent primarily depends on the seasonal, intra-annual and decadal variability in SST in the Bay of Bengal. Sea surface temperature is directly linked to the formation of two important phenomena in the Bay of Bengal, tropical cyclones and mesoscale eddies. Understanding these phenomena in this region requires knowledge on the long term trends and temporal variability of the SST.
SST is often analysed using statistical techniques for dependent data (time series).The analysis of time series is an important and valuable approach adopted in several studies, for its ability to improve the spatial and temporal resolution of the major seasonal and inter-annual patterns in biological and oceanographic data (Vantrepotte and Mélin, 2010, Vantrepotte and Mélin, 2011). These studies provide indicators about long-term changes in natural conditions, such as climate change, which is why such indicators are advised or even mandated for coastal water monitoring programmes (e.g. Water Framework Directive, 2000/60/EC) (Vantrepotte and Mélin, 2010).

2.Study area

This study has been done over Bay of Bengal from 5°N latitude to 25°N latitude,80°E longitude to 100°E longitude.

#required libraries
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
theme_set(theme_bw())
#
world <- ne_countries(scale = "medium", returnclass = "sf")
#
ggplot(world)+geom_sf(fill="darkslategray",color="grey" )+coord_sf(xlim = c(80,100),ylim = c(5,25),expand = FALSE)+xlab("Longitude")+ylab("Latitude")+annotation_scale(location="bl",width_hint=0.5)+annotation_north_arrow(location="bl",which_north="true",pad_x = unit(0.55, "in"), pad_y = unit(0.4, "in"),style = north_arrow_fancy_orienteering)+annotate(geom = "text", x = 90, y = 18, label = "Bay of Bengal", fontface = "italic", color = "grey22", size = 5)+annotate(geom = "text", x = 90, y = 24, label = "Bangladesh", color = "light blue", size = 2.5)+annotate(geom = "text", x = 83, y = 22, label = "India", color = "light blue", size = 3)+annotate(geom = "text", x = 96, y = 21, label = "Myanmar",color = "light blue", size = 3)

fig.1. geographical map of the Study area.

3.Data and Methodology

To examine the linear trend and temporal variability of sea surface temperature in our area of interest, we will use NOAA Extended Reconstructed Sea Surface Temperature (ERSST), version 5 data from their website at https://psl.noaa.gov/ . The ERSST data cover the world ocean with a spatial coverage of 2∘ latitude × 2∘ longitude global grid (89×180). The SST data is available at NetCDF format from 1856-present. Then the SST data will be extracted in the region of 80∘E~100∘E, 5∘N~25∘N. The monthly mean data from 1920 to 2020 will be applied in this study.

3.1 Time series analysis

A Time series \({{Y_t}_(t=1) }= 1^N\) is a collection of observations indexed by time t. In this study, time series were extracted for the period from 1 January 1920 to 1 January 2020, covering a time horizon of 100 years (N = 1201 monthly mean observations). Total SST observation in the time span is 132110 and the missing value is 42035 (31.82%).

3.2.1 Time series modeling

Classical decomposition methods have been widely identified in the marine sciences literature (e.g. Loisel et al., 2014; Mélin et al., 2011; Vantrepotte and Mélin, 2010). With these methods, a time series is described using the following components: a trend component \((T_t)\), which consists of the underlying long-term direction over the time horizon; a seasonal component \((S_t)\), which is a repetitive pattern over time; and an irregular component \((I_t)\), which is the unexplained variation of the time series that is not attributed to trend or seasonality. An additive model \(Y_t = T_t + S_t + I_t\) is chosen to decompose the monthly time series into these components. The strong seasonality of some time series makes it difficult to measure their trend (DeLurgio, 1998). Therefore, identification of the seasonal pattern \((S_t)\), is an essential step so that the seasonal variation, can be eliminated without affecting the trend. This new time series is given by \(Y_t⁎ = Y_t − S_t = T_t + I_t\) , t = 1 , … , N.

3.2.2 Trend analysis

A Trend is a long term change in the mean level. The simplest trend is the familiar “linear trend + noise”,for which the observation at time t is a random variable \(X_t\) given by \(X_t=α+βt+ε_t\) where \(α\) and \(β\) are constants and \(ε_t\) denotes a random error term with zero mean.The mean level at time t is given by \(m_t=(α+βt)\); this is sometimes called the ‘trend term’.Other writers prefer to describe the slope \(β\) as trend; the trend is then the change in the mean level per unit time. The current study applies the usual approach composed of adjusting the seasonality from the time series and apply a linear model to the seasonally adjusted time series.

Required libraries

#required libraries
library(ncdf4)
library(lubridate)
library(ggplot2)
library(tidyverse)
library(ggthemes)
library(forecast)
library(reshape2)
library(zoo)
library(ggfortify)

NetCDF Data loading and variable extracting

####
##NOAA ERSST v5 data set of 2*2 grid
nc <- nc_open("sst.mnmean.nc")

#extracting variables(lon,lat,time)
lon <- ncvar_get(nc,"lon")
lat <- ncvar_get(nc,"lat")
time <- ncvar_get(nc,"time")

time <- as.Date(time,origin="1800-1-1 00:00:00",tz="UTC")

##subset the variables from the interested area and time  range
lon_lim <- c(80,100)
lat_lim <- c(5,25)
time_lim <- c(ymd("1920-01-01"),ymd("2020-01-01"))

lon_ind <- which(lon >= lon_lim[1] & lon <= lon_lim[2])
lat_ind <- which(lat >= lat_lim[1] & lat <= lat_lim[2])
time_ind <- which(time >= time_lim[1] & time <= time_lim[2]) 

##extract the sst variable
sst <- ncvar_get(nc,"sst",start = c(lon_ind[1],lat_ind[1],time_ind[1]),count = c(length(lon_ind),length(lat_ind),length(time_ind)))

##make a time series of sst
sst_ts <- apply(sst,3,mean,na.rm=TRUE)
#time-span
year <- time[time_ind]

##make a data-frame with time series and observation 
tempseries <- data.frame(year=year,sst=sst_ts)

Plotting the time series

##plot the time series 
tempseries %>% ggplot(aes(x=year,y=sst_ts))+geom_line(alpha=0.8)+scale_x_date(limits =ymd(c("1920-01-01","2020-12-01")) ,breaks = seq(ymd("1920-01-01"),ymd("2020-12-01"),"10 years"),minor_breaks = "1 years",date_labels ="%Y")+labs(title = "Monthly mean SST in the Bay of Bengal from 1920 to 2020", x="year",y="Temperature" )+ scale_y_continuous(breaks = seq(20,35,0.5))+theme_clean()

fig.2. Time series plot of SST (1920-2020).

Monthly mean,Annual mean and five years moving average

mov_avg <- tempseries %>% select(year,sst) %>% mutate(sst_1yr = rollmean(sst, k = 13, fill = NA, align = "right"),sst_5yr = rollmean(sst, k = 61, fill = NA, align = "right"))

mov_avg %>% gather(key="metrice",value = "value",sst:sst_5yr)%>% ggplot(aes(x=year,y=value,col=metrice))+geom_line()+scale_color_manual(values = c("bisque4","darkred","blue"),labels=c("Monthly mean","Annual mean","5 years moving Average"))+scale_x_date(limits =ymd(c("1920-01-01","2020-12-01")) ,breaks = seq(ymd("1920-01-01"),ymd("2020-12-01"),"10 years"),date_labels ="%Y")+ scale_y_continuous(breaks = seq(23,33,0.5))+labs( x="year",y="SST [°C]" )+theme_clean(base_size = 12,)+theme(legend.title = element_blank(),legend.position = c("top"),legend.direction = "horizontal")

Decomposing the components of the additive time series

#Now,lets make a time series object
ts <- ts(sst_ts,start = c(1920,1),end = c(2020,1),frequency = 12,class = "ts")

#decomposition of ts object into trend,seasonality and error by additive model
decomposed <- decompose(ts,type = "additive")

#plot the components of an additive time series
theme_set(theme_bw())
autoplot(decomposed)

fig.3. Additive classical decomposition of SST.

Plot the linear trend of seasonally adjusted time series

#remove seasonality from the ts
decomposed_trend <- ts - decomposed$seasonal

#make a dataframe of the trend 
set.seed(12345)
decomposed_trend_df <- data.frame(date =as.Date(as.yearmon(time(decomposed_trend))),sst=as.matrix(decomposed_trend))
decomposed_trend_df %>% ggplot(aes(x=date,y=sst))+geom_line(alpha=0.8)+geom_smooth(method= "lm",se=FALSE,col="red")+geom_text(x =1950 , y = 29.5, label = "y = 28.37 + 0.000027 x , r² = 0.486",hjust=1.5,colour="blue")+scale_x_date(limits =ymd(c("1920-01-01","2020-12-01")) ,breaks = seq(ymd("1920-01-01"),ymd("2020-12-01"),"10 years"),minor_breaks = "1 years",date_labels ="%Y")+labs(x="year",y="SST [°C]" )+ scale_y_continuous(breaks = seq(25,32,0.2))+theme_clean()

fig.4. Linear fit of the seasonally adjusted time series.

Temporal variability of Seasonality

ts_df <- data.frame(time= year,sst=sst_ts) 

ts_df <- ts_df %>% mutate(year=year(time)) %>% mutate(month= month(time,label=TRUE)) %>% mutate(season = case_when(month %in% c("Nov","Dec","Jan","Feb") ~ "Winter", month %in% c("Mar","Apr","May") ~ "Spring", month %in% c("Jun","Jul","Aug") ~ "Summer", month %in%  c("Sep","Oct") ~ "Fall", TRUE ~ NA_character_) ) %>% select(-time)  %>%  as_tibble()

ts_df$season <- as.factor(ts_df$season)

ts_df %>% ggplot(aes(x=month,y=sst,fill=season)) + geom_boxplot(position = position_dodge(width = 0.7))+scale_y_continuous(breaks = seq(20,35,0.5))+labs(x="month",y="SST [°C]")+scale_fill_manual(values = c("antiquewhite4","darkolivegreen4","chocolate4","cornflowerblue"))+theme_clean()+theme(legend.title = element_blank())

fig.5. The seasonal cycle by month of SST