Introduction

Empirical orthogonal function (EOF) analyses are normally used to study possible spatial patterns of climate variability and how they change with time.In EOF analysis, one also projects the original climate data on an orthogonal basis. However, this orthogonal basis is derived by computing the eigenvectors of a spatially weighted anomaly covariance matrix, and the corresponding eigenvalues provide a measure of the percent variance explained by each pattern. Therefore, EOFs of a space-time physical process can represent mutually orthogonal space patterns where the data variance is concentrated, with the first pattern being responsible for the largest part of the variance, the second for the largest part of the remaining variance, and so on.

Study area

Bay of Bengal,part of Indian Ocean has been selected on this study to EOF analysis. The selected latitude is 5 N to 25 N and longitude is 80 E to 100 E.

Data

In this study, we will use data from International Research Institute for Climate and Society(IRI) website. Data link

Reference
Reynolds, R.W., N.A. Rayner, T.M. Smith, D.C. Stokes, and W. Wang, 2002: An Improved In Situ and Satellite SST Analysis for Climate. J. Climate, 15, 1609-1625.

##required libraries
library(ncdf4)
library(fields)
library(RColorBrewer)
library(maps)

Code

##Install the necessary package for NetCDF

#Import the netCDF file
nc <- nc_open("NOAA NCEP EMC CMB GLOBAL Reyn_SmithOIv2 monthly ssta Data Files.nc")
print(nc,1)
## File NOAA NCEP EMC CMB GLOBAL Reyn_SmithOIv2 monthly ssta Data Files.nc (NC_FORMAT_CLASSIC):
## 
##      1 variables (excluding dimension variables):
##         float ssta[X,Y,T]   
##             pointwidth: 1
##             missing_value: NaN
##             units: degree_Celsius
##             expires: 1601856000
##             CE: 4.19999980926514
##             colormap: [null 12632256 8388608 [16711680 66] [16760576 43] [15453831 16] 13959039 [13959039 7] 64636 [65535 4] [36095 24] [255 31] [128 60] [128 7] 2763429]
##             colorscalename: sstacolorscale
##             scale_max: 4.19999980926514
##             CS: -4.19999980926514
##             scale_min: -4.19999980926514
##             ncolor: 254
##             maxncolor: 254
##             long_name: Sea Surface Temperature Anomaly
##             standard_name: sea_surface_temperature
##             CI: 0.5
## 
##      3 dimensions:
##         T  Size:466
##             standard_name: time
##             pointwidth: 1
##             calendar: 360
##             gridtype: 0
##             units: months since 1960-01-01
##         Y  Size:180
##             pointwidth: 1
##             standard_name: latitude
##             gridtype: 0
##             units: degree_north
##         X  Size:360
##             pointwidth: 1
##             standard_name: longitude
##             gridtype: 1
##             units: degree_east
#set coordinate variable: X,Y
lat <- ncvar_get(nc,"Y")
lon <- ncvar_get(nc, "X")
#extract sst anomalies:
anom <- ncvar_get(nc,"ssta")

##Set the regio for Bay of Bengal : lat(5:25),lon(80:100)
lat_rng <- c(5,25)
lon_rng <- c(80,100)

lat_indx <- which(lat >= lat_rng[1] & lat <= lat_rng[2])


lon_indx <-which(lon >= lon_rng[1] & lon <= lon_rng[2])

##extract the bay of bengal sst anomalies

anom_bob <- ncvar_get(nc,"ssta",start = c(lon_indx[1],lat_indx[1],1),count = c(length(lon_indx),length(lat_indx),-1))

dim(anom_bob)
## [1]  20  20 466
#define which location is ocean or land
s1 <- which(is.na(anom_bob[,,1]))
s2 <- which(!is.na(anom_bob[,,1]))
print(length(s2))
## [1] 274
##out of 7*10=70 grid cells,there are 17 cells on the land where no data are available.

#vectorize the SST anomalies
ssta <- matrix(0,nrow = dim(anom_bob)[3],ncol = length(s2))
for(i in 1:dim(anom_bob)[3])
  ssta[i,] <- anom_bob[,,i][-s1]

###Detect the dominant pattern
#extract the eofs of data
eof <- svd(ssta)$v

####

#define the location 
loc <- as.matrix(expand.grid(x=lon[lon_indx],y=lat[lat_indx]))[s2,]

coltab <- colorRampPalette(brewer.pal(9,"BrBG"))(2048)

#######plot the first eof
par(mar=c(4,4,2,3),oma=c(1,1,1,1))

quilt.plot(loc,eof[,1],nx=length(lon[lon_indx]),ny=length(lat[lat_indx]),xlab="longitude",ylab="latitude",main="1st EOF",col = coltab)

map(database = "world",fill = TRUE,col="gray",ylim = c(5,25),xlim = c(80,100),add = TRUE)

##plot the 2nd EOF
par(mar=c(4,4,2,3),oma=c(1,1,1,1))

quilt.plot(loc,eof[,2],nx=length(lon[lon_indx]),ny=length(lat[lat_indx]),xlab="longitude",ylab="latitude",main="2nd EOF",col = coltab)

map(database = "world",fill = TRUE,col="gray",ylim = c(5,25),xlim = c(80,100),add = TRUE)