Before start

System requeriment: CDO

In our R session, we will need to - define the working directory: setwd("/home/Earth/nperez/MEDSCOPE/AemetEida/STEP1") - install and load R packages

library(s2dverification)
library(multiApply)

Workflow of aemetEida branch in CSTools

The functions PCA, NormalizeVar and SetWeight4LonLat should be loaded by running the following lines:

source("./R_PROYECTO/workspace/MEDSCOPE_v2.0/R/PCA.R")
source("./R_PROYECTO/workspace/MEDSCOPE_v2.0/R/NormalizeVar.R")
source("./R_PROYECTO/workspace/MEDSCOPE_v2.0/R/SetWeight4LonLat.R")

The first step is to define the region in which NAO will be computed.

coord <- list(lonmin = -80, lonmax = 40, latmin = 25, latmax = 85)

The output is:

> coord
$lonmin
[1] -80

$lonmax
[1] 40

$latmin
[1] 25

$latmax
[1] 85

Then, the Load function from s2verification can be used to read the monthly reanalysis of geopotential at 500 hPa for the coordinates previously defined

Note: This file has been preprocessed data with Script: $HOME/proyectos/MEDSCOPE/z500mm_ei.ksh

The parameters leadtimemin and leadtimemax have been defined taken into account that the origin of the data contained in the file .nc is year 1975 and the resolution of the data is monthly. The period 1993-2016 has been selected, however it includes from December 1993 to March 2017 the parameters could be calculated as leadtimemin = (1993 - 1975 + 1) * 12 leadtimemax = (2017 - 1975) * 12 + 3

reanal <- list(name = "erainterim", path = 
"./DATOS_PROYECTO/MEDSCOPE/MEDSCOPE_v2.0/FicherosDatos/Reanalisis/$VAR_NAME$500_GLOBAL0_1.00deg_ei_mm.nc")

data <- Load(var = "z", obs = list(reanal), sdates = c('19750101'), 
             output = 'lonlat', storefreq = 'monthly', lonmin = coord$lonmin, 
             lonmax = coord$lonmax, latmin = coord$latmin, latmax = coord$latmax, 
             leadtimemin = 228, leadtimemax = 507)

The content of the data can be reviewed visualizing the strecture:

> str(data$exp)
 NULL
> str(data$obs)
 num [1, 1, 1, 1:280, 1:61, 1:121] 49876 50381 49640 49268 51336 ...
 - attr(*, "dimensions"*)= chr [1:6] "dataset" "member" "sdate" "ftime" ...
> str(data$Dates$start)
 POSIXct[1:280], format: "1993-12-01" "1994-01-01" "1994-02-01" "1994-03-01" ...

The SeasonSelect function from ClimProjDiags package is used to select winter from December to February.

dataDJF <- SeasonSelect(data$obs, season = 'DJF', dates = data$Dates$start, timedim = 4, calendar = 'gregorian')

Funtion Apply allows to compute the climatology and the anomalies:

clim <- Apply(dataDJF$data, target_dims = list(c('ftime')), fun = mean, na.rm = TRUE)             
anom <- Apply(data = list(dataDJF$data, clim$output1), target_dims = list(c('lat', 'lon'), 
              c('lat', 'lon')), fun = function(x, y) {x - y})
anom <- aperm(anom$output1, c(3, 1, 2))/9.80665
names(dim(anom)) <- c("ftime", "lat", "lon")

The weights for latitudes and PCA can be computed:

aWeightLonLat <- SetWeight4LonLat(mWeight = 4 ,lat = data$lat, lon = data$lon)
lPCA <- PCA(anomVar = anom, lon = data$lon, lat = data$lat, neofs = 3, 
            wght = aWeightLonLat, normVariance = TRUE, 
            projwanom = FALSE)

The normalization of Principal Components can be perform by running: Note: the check on NormalizeVar should be removed to make it work.

names(dim(lPCA$PCs)) <- c("ftime", "PC")
PCstand <- Apply(data = list(lPCA$PCs, lPCA$PCs), target_dims = list(c(1), c(1)), 
       fun = NormalizeVar, output_dims = "ftime")

Plotting:

The geopotential data is loaded again for all the northern hemisphere.

dataNH <- Load(var = "z", obs = list(reanal), sdates = c('19750101'), 
    output = 'lonlat', storefreq = 'monthly', lonmin = 0, lonmax = 360, 
        latmin = 0, latmax = 90, leadtimemin = 228, leadtimemax = 507)

Again, the anomalies and the units correction is applied:

dataNHDJF <- SeasonSelect(dataNH$obs, season = 'DJF', dates = data$Dates$start, timedim = 4, calendar = 'gregorian')
climNHDJF <- Apply(dataNHDJF$data, target_dims = list(c('ftime')), fun = mean, na.rm = TRUE)            
anomNHDJF <- Apply(data = list(dataNHDJF$data, climNHDJF$output1), target_dims = list(c('lat', 'lon'),
              c('lat', 'lon')), fun = function(x, y) {x - y})
anomNHDJF <- aperm(anomNHDJF$output1, c(3, 1, 2))/9.80665
names(dim(anomNHDJF)) <- c("ftime", "lat", "lon")

Then, the covariance between the geopotential and the PC1 (which is one time series)

cova <- Apply(data = list(anomNHDJF, PCstand$output1), target_dims = list(c('ftime'), c('ftime')), fun = cov, use = "pairwise.complete.obs")

The output looks like:

> str(cova)
List of 1
 $ output1: num [1:91, 1:360, 1:3] 16.7 18 19.5 21.3 23.1 ...

Two kind of plots are possible with s2dverification functions:

PlotStereoMap(cova$output1[,,1], lon = dataNH$lon, lat = dataNH$lat, fileout ="CovPC1_stereo.png")
PlotEquiMap(cova$output1[,,1], lon = dataNH$lon, lat = dataNH$lat, filled.continents = FALSE, 
            fileout ="CovPC1_equip.png")

Stereo

Rectangular