Introduction

Sea surface temperature is the temperature of the top few millimeter of the ocean’s surface. An anomaly is a departure from average conditions.Some sea surface temperature anomalies are simply transient events,not part of a specific pattern or trend.Other anomalies are more meaningful. Strong,localized sea surface temperature anomalies may reveal that an ocean current,such as the Gulf stream current off the east coast of United States,has veered off its usual path for a time or is stronger or weaker than than usual.Sea surface temperature anomalies that persist over many years can be signal of regional or global climate change such as global warming.Sea surface temperature anomalies have practical as well as scientific applications. For example,in coastal areas,anomalous temperatures(either cool or warm) can favour one organism in an ecosystem over another,causing populations of one kind of bacteria,algae,or fish to thrive or decline.Warm sea surface temperaure anomalies can also warn coral reef ecosystem.
reference: NASA earth observatory

Data

We will obtain world wide Sea Surface Temperature anomalies data from NASA’s Goddard Institute for Space studies website and our data(1880-2020) link.

Short summary

After loading the required libraries and data, we will extract the variables(lon,lat,time) from our netCDF file and subset the lat,lon as we want to get data from Bay of Bengal in the North Indian Ocean.We select latitude from 16N to 23N and longtude from 88E to 92E.Then we will plot the timeseries of that area.

Code

#libraries
library(ncdf4)
library(fields)

#load the netCDF data from the working directory
nc<- nc_open("gistemp1200_GHCNv4_ERSSTv5.nc")

#print the nc file to get info of the data 
print(nc,1)
## File gistemp1200_GHCNv4_ERSSTv5.nc (NC_FORMAT_CLASSIC):
## 
##      2 variables (excluding dimension variables):
##         int time_bnds[nv,time]   
##         short tempanomaly[lon,lat,time]   
##             long_name: Surface temperature anomaly
##             units: K
##             scale_factor: 0.00999999977648258
##             cell_methods: time: mean
##             _FillValue: 32767
## 
##      4 dimensions:
##         lat  Size:90
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##         lon  Size:180
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##         time  Size:1687
##             long_name: time
##             units: days since 1800-01-01 00:00:00
##             bounds: time_bnds
##         nv  Size:2
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named nv BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## 
##     5 global attributes:
##         title: GISTEMP Surface Temperature Analysis
##         institution: NASA Goddard Institute for Space Studies
##         source: http://data.giss.nasa.gov/gistemp/
##         Conventions: CF-1.6
##         history: Created 2020-08-13 09:26:17 by SBBX_to_nc 2.0 - ILAND=1200, IOCEAN=NCDC/ER5, Base: 1951-1980

Quick data review

From the info above,our data file is NC FORMAT CLASSIC.It has two variables including the Surface temperature anomaly(short tempanomay) which has 3 dimension(lon,lat,time).Given Unit of surface temperature anomaly is Kelvin.

Code continuing..

#extract the variables
lat<- ncvar_get(nc,"lat") #latitude
lon<- ncvar_get(nc,"lon") #longitude
time<- ncvar_get(nc,"time") #time(1880-2020)
time<- as.Date(time,origin="1800-01-01",tz="UTC") #convert into date format
summary(time)
##         Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
## "1880-01-15" "1915-03-01" "1950-04-15" "1950-04-15" "1985-05-30" "2020-07-15"
#plot the first date of the data
anom<- ncvar_get(nc,"tempanomaly",start = c(1,1,1),count = c(-1,-1,1))
image.plot(lon,lat,anom,col = rev(rainbow(100)))
map(add = TRUE)
title(main = "SST anomalies of 1880-01-15")

#Now,subsetting lat-lon
lat_rng <- c(16,23)
lon_rng <- c(88,92)

lat_ind <- which(lat >= lat_rng[1] & lat <= lat_rng[2])
lon_ind <- which(lon >= lon_rng[1] & lon <= lon_rng[2])
#extract the timeseries
anom_ts <- ncvar_get(nc,"tempanomaly",start = c(lon_ind[1],lat_ind[1],1),count = c(length(lon_ind),length(lat_ind),-1))

anom_ts <- apply(anom_ts,3,mean,na.rm=TRUE)

plot(time,anom_ts,type="l",xlab="year",ylab="Temperature anomalies",main="Temp. anomalies of N. BoB from 1880 to 2020")

More customized plot Using ggplot

library(tidyverse)
library(ggplot2)#load ggplot          
library(lubridate)        

#now we want to make a dataframe              
tempseries <- data.frame(year=time,anomalies=anom_ts)
tempseries %>% ggplot(aes(x=time,y=anom_ts))+geom_line()+scale_x_date(limits =ymd(c("1880-01-15","2025-07-15")) ,breaks = seq(ymd("1880-1-15"),ymd("2020-07-15"),"20 years"),minor_breaks = "10 years",date_labels ="%Y")+labs(title = "Temperature anomalies of N.BoB from 1880 to 2020", x="year",y="Temperature anomalies" )+ scale_y_continuous(breaks = seq(-2,2,0.5))+theme_classic()

Now,we want to make a grid of 30 years interval data row by row using ggplot.We subset 30 years data and plot it by rows in a single plot.

lat_rng <- c(16,23) #bob latitude
lon_rng <- c(88,92) #bob longitude

tm_rng_1880_1915 <- ymd(c("1880-01-01","1915-01-01"))
tm_rng_1915_1950 <- ymd(c("1915-01-01","1950-01-01"))
tm_rng_1950_1985 <- ymd(c("1950-01-01","1985-01-01"))
tm_rng_1985_2020 <- ymd(c("1985-01-01","2020-01-01"))


lat_ind <-which(lat >= lat_rng[1] & lat <= lat_rng[2])
lon_ind <- which(lon >= lon_rng[1] & lon <= lon_rng[2])

tm_ind_1880_1915 <- which(time> tm_rng_1880_1915[1] & time< tm_rng_1880_1915[2])
tm_ind_1915_1950 <- which(time> tm_rng_1915_1950[1] & time< tm_rng_1915_1950[2])
tm_ind_1950_1985 <- which(time> tm_rng_1950_1985[1] & time< tm_rng_1950_1985[2])
tm_ind_1985_2020 <- which(time> tm_rng_1985_2020[1] & time< tm_rng_1985_2020[2])
  
y.1880_1915 <- time[tm_ind_1880_1915]
y.1915_1950 <- time[tm_ind_1915_1950]
y.1950_1985 <- time[tm_ind_1950_1985]
y.1985_2020 <- time[tm_ind_1985_2020]

a.1880_1915 <- ncvar_get(nc,"tempanomaly",start =c(lon_ind[1],lat_ind[1],tm_ind_1880_1915[1]),count = c(length(lon_ind),length(lat_ind),length(tm_ind_1880_1915)) )
a.1915_1950 <-ncvar_get(nc,"tempanomaly",start =c(lon_ind[1],lat_ind[1],tm_ind_1915_1950[1]),count = c(length(lon_ind),length(lat_ind),length(tm_ind_1915_1950)) )
a.1950_1985 <-ncvar_get(nc,"tempanomaly",start =c(lon_ind[1],lat_ind[1],tm_ind_1950_1985[1]),count = c(length(lon_ind),length(lat_ind),length(tm_ind_1950_1985)) )
a.1985_2020 <-ncvar_get(nc,"tempanomaly",start =c(lon_ind[1],lat_ind[1],tm_ind_1985_2020[1]),count = c(length(lon_ind),length(lat_ind),length(tm_ind_1985_2020)) )


t.1880_1915 <- apply(a.1880_1915,3,mean,na.rm=TRUE)
t.1915_1950 <- apply(a.1915_1950,3,mean,na.rm=TRUE)
t.1950_1985 <- apply(a.1950_1985,3,mean,na.rm=TRUE)
t.1985_2020 <- apply(a.1985_2020,3,mean,na.rm=TRUE)


###make a dataframe with all the subsets

anom_df <- data.frame(t.1880_1915=t.1880_1915,t.1915_1950=t.1915_1950,t.1950_1985=t.1950_1985,t.1985_2020=t.1985_2020,y.1880_1915=y.1880_1915,y.1915_1950=y.1915_1950,y.1950_1985=y.1950_1985,y.1985_2020=y.1985_2020)
library(cowplot)
library(gridExtra)

p1<- anom_df%>% ggplot(aes(x=y.1880_1915,y=t.1880_1915))+geom_line()+scale_x_date(limits =ymd(c("1880-01-15","1914-12-15")) ,breaks = seq(ymd("1880-1-15"),ymd("1914-12-15"),"5 years"),minor_breaks = "1 years",date_labels ="%Y")+xlab("")+ylab("") +scale_y_continuous(breaks = seq(-2,2,1))+theme_classic()
p2<- anom_df%>% ggplot(aes(x=y.1915_1950,y=t.1915_1950))+geom_line()+scale_x_date(limits =ymd(c("1915-01-15","1949-12-15")) ,breaks = seq(ymd("1915-01-15"),ymd("1949-12-15"),"5 years"),minor_breaks = "1 years",date_labels ="%Y")+xlab("")+ylab("") +scale_y_continuous(breaks = seq(-2,2,1))+theme_classic()
p3<- anom_df%>% ggplot(aes(x=y.1950_1985,y=t.1950_1985))+geom_line()+scale_x_date(limits =ymd(c("1950-01-15","1984-12-15")) ,breaks = seq(ymd("1950-01-15"),ymd("1984-12-15"),"5 years"),minor_breaks = "1 years",date_labels ="%Y")+xlab("")+ylab("") +scale_y_continuous(breaks = seq(-2,2,1))+theme_classic()
p4<- anom_df%>% ggplot(aes(x=y.1985_2020,y=t.1985_2020))+geom_line()+scale_x_date(limits =ymd(c("1985-01-15","2019-12-15")) ,breaks = seq(ymd("1985-01-15"),ymd("2019-12-15"),"5 years"),minor_breaks = "1 years",date_labels ="%Y")+xlab("")+ylab("") +scale_y_continuous(breaks = seq(-2,2,1))+theme_classic()

p_final<- plot_grid(p1,p2,p3,p4,label_size = 10,label_colour = "grey",label_fontfamily = "serif",ncol = 1,align = "vh",vjust = 1,scale = 1,rel_heights = c(1.5,1.5,1.5,1.5))

y.grobe <- textGrob("Temperature anomalies[k]",gp=gpar(col="black",fontsize=12),rot=90)
x.grobe <- textGrob("Year",gp=gpar(col="black",fontsize=12))

grid.arrange(arrangeGrob(p_final,left = y.grobe,bottom = x.grobe))